Forests on drained agricultural peatland are potentially 1 large sources of greenhouse gases – insights from a full 2 rotation period simulation 3

18 The CoupModel was used to simulate a Norway Spruce forest on fertile drained peat over 60-19 years, from planting in 1951 until 2011, describing abiotic, biotic and greenhouse gas (GHG) 20 emissions (CO 2 and N 2 O). By calibrating the model against tree ring data we obtained a 21 ‘reference’ model by which we were able to describe the fluxes and controlling factors over 22 the 60 years. We discuss some conceptual issues relevant to improving the model in order to 23 better understand peat soil simulations. However, the present model was able to describe the 24 most important ecosystem dynamics such as the plant biomass development and GHG 25 emissions. The GHG fluxes are composed of two important quantities, the forest carbon (C) 26 uptake, 405 g C m -2 yr -1 and the decomposition of peat soil, 396 g C m -2 yr -1 . N 2 O emissions 27 contribute to the GHG emissions by 0.5 g N m -2 yr -1 , corresponding to 56.8 g C m -2 yr -1 . The 28 60-year-old Spruce forest has an accumulated biomass of 164 Mg C ha -1 . However,


Introduction
Peatlands contain around one third of the carbon (C) stored in global soils, which is equivalent to almost half that present in the atmosphere (FAO, 2012;IPCC, 2013).Undisturbed peatlands accumulate C as partially decayed vegetation, and the decay processes emit C in the form of carbon dioxide (CO 2 ) and methane (CH 4 ).Overall, the balance of photosynthesis and respiration in peatlands means that these systems act as C sinks, acting to mitigate climatic warming (e.g.Gorham, 1991).However, when peatlands are drained for intensified land use, i.e. agriculture or forestry, the stored peat starts to decompose aerobically.The accelerated soil decomposition emits large amounts of CO 2 , in contrast CH 4 emissions are greatly reduced, possibly even accounting for a net uptake of atmospheric CH 4 (Limpens et al., 2008).The decomposition also releases nitrogen, and another powerful GHG, nitrous oxide (N 2 O), could also be produced, primarily through microbial nitrification and denitrification processes (Firestone and Davidson, 1989).Globally, peatlands cover only 3 % of the Earth surface, of which in turn 10-20 % have been drained for agriculture or forestry, mainly in the boreal and tropical regions (FAO, 2012).However, these drained areas emit around 6 % of the global annual anthropogenic GHG emissions (IPCC, 2013).
To date, a number of studies have investigated the size of GHG fluxes from managed peatlands with different land uses, together with their interactions with environmental factors (Kasimir Klemedtsson et al., 1997;Von Arnold et al., 2005a, b;Alm et al., 2007;Beek et al., 2010;Lund et al., 2010;Lohila et al., 2011;Ojanen et al., 2013).Several factors have been found to influence the size of the emissions, including the groundwater level (GWL), land use intensity, climate zones, and soil fertility (Klemedtsson et al., 2005;Drösler et al., 2008;Leppelt et al., 2014).In general, nutrientrich fens with deep GWL are larger GHG sources than ombrotrophic bogs with shallow GWL, while intensive land use in tropical and/or temperate regions have much higher emissions than extensive land use in boreal regions (Byrne et al., 2004).Peatlands in Europe used as grassland, agricultural land, peat cuts, and abandoned peat are generally found to be net GHG sources (Byrne et al., 2004;Drösler et al., 2008).However, forested drained peatland can be everything from a source to a small GHG sink due to the growing forest, where the net primary production (NPP) of trees and understorey vegetation balances the soil emissions (Drösler et al., 2008;Klemedtsson et al., 2008;Hommeltenberg et al., 2014).Previous flux measurement studies have also shown contradictory results.Measurements from Scandinavia and Great Britain have shown the NPP to compensate for the soil CO 2 release, and thus the forests to act as net sinks (Hargreaves et al., 2003;Von Arnold et al., 2005a, b;Ojanen et al., 2013).Hommeltenberg et al. (2014) also reported an afforested drained bog in southern Germany to be a net GHG sink; however, if the 44-year history of the forest were included in the analysis, then the so-called "long-term carbon balance", showed the forest to be an overall GHG source.Von Arnold et al. (2005a) showed that accounting for N 2 O in the greenhouse budget calculation could shift drained birch peatlands from being minor GHG sinks into sources.This was also shown by Meyer et al. (2013) for a drained former agricultural peat soil with spruce forest, where soil N 2 O emissions, in terms of global warming potential (265 times of CO 2 in a 100-year perspective, IPCC, 2013), offset half the net ecosystem exchange (NEE).Large N 2 O emissions are most pronounced for fertile soils like former agricultural peatlands (Klemedtsson et al., 2005).So far most studies have only covered a few years at most.Consequently we still lack an understanding of the full GHG balance when viewed over the full forest rotation (Maljanen et al., 2010).
In the present study we aim to address this knowledge gap by exploring the GHG balance for a Norway spruce (Picea abies) forest on drained agricultural peatland (Skogaryd Research Site) over a full rotational time period.Since measurements are mostly short term, and because it is not possible to directly upscale the measured fluxes to the entire forest rotation period (Drösler et al., 2008;Hommeltenberg et al., 2014), we chose a modelling approach based on emission data over 5 years and data on forest growth rate over 45 years for a spruce forest on former agricultural peatland.

Site description
Data used for the present study were obtained from the Skogaryd research site (http://www.fieldsites.se/en/field-research-stations), located in southwest Sweden (58 • 23 N, 12 • 09 E), which is part of the Swedish Infrastructure for Ecosystem Science (SITES, www.fieldsites.se).The drained peat area at Skogaryd was previously a fen, classified as mesotrophic peat with a peat depth of more than 1 m, according to the soil classification scheme suggested by Karlsson (1989).It was initially drained by ditches in the 1870s and then used for agriculture until 1951.Norway spruce (P.abies) was then planted and the stand is now a mature mixed coniferous forest dominated by Norway spruce (95 % by stem volume), with a sparse presence of Scots pine (Pinus sylvestris) and Silver birch (Betula pubescens; Klemedtsson et al., 2010).The site has been intensively measured and monitored since 2006, providing abiotic and biotic data including CO 2 and N 2 O fluxes that could be used to validate the long-term model predictions.More detailed site description can be found in He et al. (2016), Klemedtsson et al. (2010), Meyer et al. (2013) and Ernfors et al. (2010).

Concept of drained peatland for forestry
When peatlands are drained for forestry or agriculture, resulting in a lower GWL, the aerobic soil volume increases (Fig. 1a).The previously water-logged peat soil then decomposes aerobically, losing soil C stock and also causing a lowering of the soil surface (surface subsidence ;Eggelsmann, 1976;Hooijer et al., 2012).During the first few decades after planting, the development of plant roots and leaf area cover increases the transpiration rate, so deepening the GWL (Fig. 1b).In other words, a growing forest will, in part, help to keep the soil drained.However, drainage becomes less efficient with time due to subsidence and ditches becoming filled with litter and moss, all of which can lead to an increased GWL (Fig. 1c), which is why ditch maintenance is performed regularly.After ditch maintenance the forest ecosystem restarts at the well-drained state (Fig. 1d), until the final clear-cutting when re-drainage has to be conducted.The entire cycle then starts again and can continue until all the peat is gone.

Brief introduction to the CoupModel
The CoupModel (coupled heat and mass transfer model for soil-plant-atmosphere systems) is an updated version of the previous SOIL and SOILN model (Jansson and Moon, 2001).The main model structure is a one-dimensional, layered soil depth profile, in which the water, heat, and C and N dynamics are simulated based on detailed descriptions of soil physical and biogeochemical processes.C and N dynamics are simulated both in the soil and in the plant, driven by the canopy-intercepted radiation, regulated by multiplicative response functions of air temperature, and plant availability of water and N. Two vegetation layers are simulated in the model, the spruce tree and the understorey layer (e.g.grasses and shrubs; He et al., 2016).The model is available at http://www.coupmodel.com/.A detailed description of the model, its parameterization and setup is given in He et al. (2016); here only the variables and parameters with different values are reported.

Model approach and design
The CoupModel conceptually divides the soil organic matter (SOM) into two pools called soil litter (fresh plant detritus) and humus, constituting a fast and a slow decomposing pool, respectively (Johnsson et al., 1987).When soil litter decays, carbon is either released as CO 2 , or added into a resistant fraction, the humus pool (Johnsson et al., 1987).
In this study, the soil humus pool was used to represent the old stored soil peat.Thus soil decomposition is composed of both peat decomposition (called humus decomposition in the model) and soil litter decomposition.Besides, CoupModel conceptualizes the soil profile into a number of soil layers, where the soil's physical structure (defined by the measured water retention characteristics) and the drainage depth (a parameter used for estimation of horizontal flow of water out of the site due to drainage) is assumed to be fixed over time (Fig. 1e and f), with the drainage depth set to 0.5 m as in He et al. (2016).Though the drainage depth is a very important parameter for the simulated GWL, a fixed drainage level is not to be confused with a fixed GWL as the latter is simulated (see Fig. 5f).The subsidence of the soil surface and any variation in drainage (Fig. 1a, b, c and d) during the plant development years (1951 to 2011) cannot explicitly be simulated.We thus make the following assumptions to simplify the system.
First, the soil layers are assumed the same over the 60 years simulated, and the soil physical characteristics in 1951 are assumed the same as measured in 2006.Whilst this assumption may not hold in detail, we consider any changes minor as (1) this site has been drained for many years (starting in 19th century), why physical soil compaction should not be important during the last 60 years, and (2) soil properties were not found to be the major GHG emission influencing factor (He et al., 2016).A range of drainag depths was used to quantify the model's sensitivity.The lower end of the range was chosen to be a drainage depth of 0.3 m, since this has been suggested to be the minimum requirement to sustain forest productivity on drained peatlands (Sarkkola et al., 2010;Ojanen et al., 2013).The higher drainage level of 0.8 m was set according to general forest management practices, also taking into consideration the maximum simulated soil depth of 1 m.
Second, in order to define the initial soil C content in 1951, we use the soil C measurements made at Skogaryd in 2007, back-calculated to 1951 by assuming an annual peat loss of 260 g C m −2 yr −1 from 1951 to 2007.This annual loss was taken from the recent IPCC wetland supplement (IPCC, 2014), where it represents the emission factor for forest on drained nutrient-rich peatlands in the temperate region.The model's sensitivity to this initial condition was assessed by varying IPCC emission factors (EFs) between 200 and 330 g C m −2 yr −1 when calculating total soil C in 1951.In addition, an extremely large initial soil C is also used in the sensitivity analysis which was back-calculated using the highest peat decomposition rate of 630 g C m −2 yr −1 (Meyer et al., 2013) measured at Skogaryd during 2008.The back calculated total soil C is assumed uniformly distributed in the soil profile of 1 m depth, based on the measured data in 2007 (He et al., 2016).
Third, the soil C / N ratio in 1951 is assumed to be the same as measured in 2006, and the N deposition rate was also assumed to be constant as in He et al. (2016) during the entire simulated period.The model's sensitivity to this was tested by varying the initial soil C / N ratio between 20 and 45, the latter being a value measured at a nearby un-drained peatland near Skogaryd.
Fourth, the model only simulates the C and N dynamics in the uppermost 1 m depth of soil.
The model was initially run with the calibrated single parameter representation using the same mean parameter values as used by He et al. (2016).However, each calibrated parameter has a range of possible values, its so-called posterior distribution, which we varied in order to fit the model results to the 45-year (1966 to 2011) tree-ring-derived biomass data and extended abiotic data (2006 to 2011).We call the model parameterized to fit those data the "vegetation fitted" model, used for sensitivity analysis by varying the drainage depth, initial soil C, as well as the initial soil C / N ratio.

Tree ring sampling and data processing
The previous calibration of the CoupModel mainly focused on the soil processes while plant development was less emphasized (He et al., 2016).In order to calibrate the model results of the plant biomass development, we acquired incremental core samples from the spruce trees in Skogaryd during spring 2013, to estimate forest biomass.In total, 25 samples were obtained from randomly chosen trees.The cores were taken at breast height (1.3 m above ground).The annual growth rings in the tree cores were cross-dated according to standard dendrochronological methods (Stokes and Smiley, 1968) to assign an exact calendar year of formation to each ring.Tree ring width data were obtained by analysis of scanned images of carefully surfaced cores using the software CooRecorder (cybis.se).The annual variation in height growth was modelled with the Korf's function using cumulative radial growth during the previous years, calibrated by extensive inventory data, collected in 2010 (Meyer et al., 2013).Since the inventory data lacked information concerning trees with a diameter smaller than 10 cm, and because the sample depth of trees decreases back in time, the forest biomass calculations were only considered to be valid from 1966 (a date when all trees had a diameter above 10 cm and the sample replication was complete).The forest biomass was calculated for stem, living branches, dead branches, stumps and roots including fine roots, following the allometric equations (Marklund, 1988) for spruce in Minkkinen et al. (2001) and Meyer et al. (2013), using the inputs of measured annually resolved radial growth and modelled annual longitudinal growth.The total biomass of the tree stands was calculated as a sum of the average biomass of the individual trees, where the planting density was assumed to be 3000 trees ha −1 , which was a typical planting density during the 1950s in Sweden (Drossler et al., 2013).A thinning was conducted by the land owner in 1979 when the number of trees was reduced to ca. 1000 trees ha −1 , according to the survey data presented in Meyer et al. (2013).Using these tree ring biomass data, the thinning management was estimated to have removed 72 % of the spruce biomass.The forest thinning practices were assumed and made according to general Swedish forest management guidelines (Svensson et al., 2008).In addition, a heavy storm hit Skogaryd forest in 2010 and blew down 10 % of the tree biomass.The fallen trees were removed from the experimental site after the storm event.Therefore an additional harvest was included in the CoupModel to simulate this removal of storm-fallen biomass.

Data for model forcing
To drive the model, we used daily mean meteorological data (1961 to 2011) from the Swedish Meteorological and Hydrological Institute (SMHI) Såtenäs station (58 • 44 N, 12 • 71 E), (www.smhi.se)situated approximately 60 km east of Skogaryd.Precipitation, air temperature, wind speed and relative humidity data from Såtenäs were strongly correlated (R 2 > 0.8) with those from Skogaryd from 2006 to 2011, and were of similar magnitude.Another driving variable needed in CoupModel is the global short-wave radiation.As these data are not available from Såtenäs station, they were deduced by the model from the potential global radiation and atmospheric turbidity, using the measured total cloud-cover fraction (for more details see http://www.coupmodel.com).Since meteorological data were only available from 1961, the meteorological data from 1961 to 1971 were duplicated to represent the climate between 1951 and 1961.

GHG budget compilation
For a total GHG budget of the system we include harvest removal and products.We assume that the biomass removed by thinning management in 1979 and the storm harvest in 2010 was mainly used for paper production, as is common practice in Sweden (Swedish Forest Agency, 2005).We therefore use the emission factors suggested in the IPCC guidelines (IPCC, 2006), in which paper is assumed to decay exponentially with a half-life of 2 years.

Plant and soil development from 1951 to 2011
The simulated tree biomass dynamics during the 60 years agrees well with the estimated tree biomass from radial growth observations beginning in 1966.After an initial phase of slow growth during the establishment of the spruce trees' leaf area, growth increased almost linearly (Fig. 2d).The spruce trees gradually increased their leaf (needles) cover until a closed canopy formed in the 1980s with a maximum leaf area index (LAI) of around 6, which was similar to field measurements (Fig. 2b).The simulated annual average spruce tree growth over the whole period is 413 g C m −2 yr −1 with the maximum growth rate of 848 g C m −2 yr −1 in 1974 (Fig. 2c).However, the "vegetation fitted" model showed a slow establishment of the spruce in the first decade due to a modelled competition from grasses and other field vegetation, thus underestimating the spruce growth before 1970, mainly caused by lack of information on initial stages.The LAI and the NPP of spruce generally follow the dynamics of the plant's ability to intercept radiation (Fig. 2a); however, the model slightly overestimates annual spruce tree growth from the 1970s to the 1990s, and underestimates it from 1996 until 2011 (Fig. 2c).Furthermore, the large increase of simulated plant growth observed in 2006 was not observed in the tree ring data.The total tree biomass in 2011 is modelled to be 16.0 kg C m −2 , which is very similar to the biomass estimated from the tree ring data, 16.2 kg C m −2 (Fig. 2d).
The thinning conducted in 1979 removed 6.8 kg C m −2 plant biomass, and the storm in 2010 caused an additional removal of 1.8 kg C m −2 ; these quantities were used for indirect emission calculations (Fig. 2d).The modelled amounts of leaf and root biomass in 2007 also match estimations using allometric equations reported by Meyer et al. (2013).The modelled and estimated values for leaf biomass were 0.95 and 1.06 kg C m −2 , respectively, and the values for total roots, both coarse roots, (> 2 mm) and fine roots, (< 2 mm) were 2.9 and 3.0 kg C m −2 , respectively.The modelled value for spruce stem biomass was 12.8 kg C m −2 , which was higher than the estimated 11.2 kg C m −2 .This discrepancy may be explained by the estimated total spruce tree biomass by Meyer et al. (2013) being smaller than that estimated from tree ring data.The maximum biomass of understorey vegetation was simulated to be around 2 kg C m −2 10 years after planting, but it decreased gradually thereafter (Fig. 2e).
Table 1 shows the soil C budget of each modelled soil layer (down to 1 m) in 1951 and 2011.The soil C content at the uppermost 5 cm layer increases due to the addition of plant litterfall (Fig. 3), where the modelled C content in the first metre of soil is shown to match the observed data.Except the deepest layer, the other soil layers all lose soil C where losses decrease by depth.This is due to a soil water content increase, where decomposition is zero in the saturated soil (like the 90-100 cm layer; Table 1).Over the whole of the simulated 60 years, the accumulated soil litter decomposition almost equaled that of the soil peat (treated as humus in the model), where ca.80 % of the litter is respired and the rest added into the resistant soil C fraction, the soil peat (called humus formation in the figure).Over the 60 years, the soil litter was close to balance as the accumulated plant litterfall almost equal to the accumulated soil litter decomposition and humus formation (Fig. 3).Thus the total losses of soil C are mostly from decomposition of historical soil peat.

Comparing vegetation fitted model output with observational data from 2006 to 2011
The simulation beginning in 1951 using the "vegetation fitted" model showed a good fit with data collected during 2006 until 2011 of GWL, total net radiation and soil temperature data.The linear correlations (R 2 ) between the simulated and measured data were all above 0.8 with the mean errors close to zero (Fig. 4).Discrepancies were found in May 2010, when the measured GWL peaked (high GWL) which by the model was underestimated (Fig. 4c), and during summers and autumns when the model overestimated both radiation and soil temperature (Fig. 4a, b).Besides showing reasonable description of abiotic factors, the model results were also similar to observed data between 2007 and 2008 on NEE flux, both in terms of seasonal pattern and magnitude (Fig. 4d).However, the simulations seem to slightly underestimate the CO 2 uptake during summertime and overestimate the respiration flux in the autumn (Fig. 4).The model performance for N 2 O emissions during 2006 to 2011 was generally similar as in the previous calibration study (He et al., 2016), where the annual emission size was reasonably simulated but the model had some difficulties in capturing every measured emission peak.

Annual NEE and N 2 O from 1951 to 2011
The annual 60-year NPP for the spruce forest, including biomass and litter, was on average 673 g C m −2 with less than 100 g C m −2 during the first 10 years after planting, and with a value that fluctuates around 1000 g C m −2 yr −1 over the last 40 years (Fig. 5b).Peat respiration (decomposition) shows a slight decreasing trend during the simulated period, with an annual average of 399 g C m −2 (Fig. 5c).The decreasing trend may be explained by a lower amount of soil peat left in the surface (Table 1 and Fig. 3) and an increasing GWL (Fig. 5f) where inter-annual variations are mainly regulated by the weather (Fig. 5a).NPP and peat decomposition are the two major components of NEE, in which the system showed itself to be both a sink and a source during the first 19 years (1951 to 1970), but thereafter to be a continuous CO 2 sink, except for 1980 and 2002 (Fig. 5d).The thinning management in 1979 had a large impact on the NEE which changed the system to that of a source of 820 g C m −2 yr −1 for the following year.After 1981, the forest ecosystem was a continuous sink of CO 2 with an average NEE of 217 g C m −2 yr −1 except for being a minor source of 82 g C m −2 yr −1 for 2002 (Fig. 5d).
The model predicts low N 2 O emissions for the initial 10 years.Instead it predicts most of the N 2 O to be emitted from 1966 to 1988, a period when the understorey vegetation becomes sparse and the spruce trees are still small.Three years (1969, 1973 and 1977) of extreme emissions were due to very high precipitation events in summer.The last 30 years the emissions were lower again, due to a closed canopy and high N uptake by trees and a more shallow GWL.Over the 60 years, the simulated annual N 2 O emission varied from a minimum of 0.01 to a maximum 7 g N m −2 yr −1 , with an average of 0.7 g N m −2 yr −1 (Fig. 5e).

Overall GHG balance from 1951 to 2011
Over the full 60-year time period the forest trees acted as a C sink and the soil as a source, of fairly similar size (Fig. 6).This could be viewed as a relocation of C from the soil to the trees, since our model predicts the total soil C loss to be 75 kg CO 2 m −2 over the 60 years, while total plant biomass (including spruce forest and understorey vegetation) sequesters 58 kg CO 2 m −2 .The accumulated NEE shows the young forest ecosystem to be a net CO 2 source, and it is not until 1990, 39 years after the forestation, that the ecosystem uptake balances previous cumulative emissions and it reaches zero CO 2 emission before becoming an overall carbon sink.If including the N 2 O emissions during the 60year rotation period, taking the most commonly used 100year time horizon global warming potential from the IPCC (1 g N 2 O = 265 g CO 2 eq, IPCC, 2013), the source strength of the forest ecosystem increases and the system switches to an overall small GHG source.
However, if including the fate of the biomass removed as thinnings, usually used for paper production, resulting in indirect CO 2 emissions from consumed paper makes this extended system (from the production site to the fate of the products) a large GHG source of 38 kg CO 2 m −2 by the end of the simulation (Fig. 6).Soon, the whole forest will be harvested releasing most of the captured carbon into the atmosphere, 16 kg C m −2 (Fig. 2d), and if everything were released from these soils there would be 96.9 kg CO 2 m −2 released over a period of 60 years. 1,2,3Back-calculated initial soil C using the reported range of IPCC EF's 200; 330 and 630 g C m −2 yr −1 respectively. 4Positive change of NEE means the forest ecosystem sequesters more atmospheric CO 2 than the "vegetation fitted" model; negative change means sequestering less atmospheric CO 2 or a possible source to the atmosphere.

Model sensitivity
Accumulated plant biomass is most sensitive to a higher soil C / N ratio or a shallower drainage depth (Table 2).The peat decomposition is instead more sensitive than the accumulated plant biomass to larger initial soil C or increasing drainage depth (Table 2).Also, magnitudes of NEE and N 2 O fluxes are very sensitive to these variations, the NEE becoming a CO 2 source at larger initial soil C, since peat decomposition rate becomes larger than the accumulated plant biomass.The model sensitivity also shows higher N 2 O emissions under shallower rather than deeper drainage (Table 2).When these various factors were combined, the peat decomposition varied by −38 to +33 %, being largest when the combination was deep drainage with the largest initial soil C, and a low initial soil C / N ratio.The accumulated biomass varied between −69 and +6 %, being smallest when the combination was shallow drainage with a low initial soil C and a large soil C / N ratio.However, the overall total GHG emissions, including the thinning and storm harvested biomass and its associated CO 2 losses, the emissions increased by 11 to 57 % (Table 2), suggesting that the total GHG balance was still a source to the atmosphere.

Comparison of our simulated results with observational and published data
The GHG balance over a rotational period for forestry on drained peatland is mainly determined by two large fluxes viz.plant C assimilation and peat decomposition.We therefore first discuss the validity of these two variables by com-paring our simulated results with values published in the literature.

Plant growth
Our simulated spruce growth at 413 g C m −2 yr −1 was higher than the normal growth rate of 162 to 270 g C m −2 yr −1 in southwest Sweden, but lower than the potential growth rate of 472 to 607 g C m −2 yr −1 under experimentally optimal nutrient conditions (Bergh et al., 2005).This high growth rate can be explained by the fertile soil at the Skogaryd site, which was a drained fen before it was used for agriculture, and then forestry.The high rate of nitrate leaching, estimated at 4.3 g N m −2 yr −1 also suggests that nutrients are not likely to be limiting.That the forest growth at this site is close to maximum has also been demonstrated in a modelling study by Tarvainen et al. (2013), who showed that if canopy N content was increased by 30 %, canopy C uptake would only increase by 2-4 % and none of the 37 nutrients tested would directly limit photosynthesis.The very small increase of plant growth (+6 %) in our model sensitivity analysis (Table 2), obtained when more deeply drained soil plus a larger initial soil C and a lower C / N ratio assumed, can also be explained by the already high fertility at the site, so any extra nutrient availability would have a negligible impact.Our simulated understorey vegetation was small during most of the simulated years; however, it dominated the organic matter dynamics and GHG fluxes in the first 2 decades after plantation, a finding similar to that of Laiho et al. (2003).

Soil CO 2 and N 2 O fluxes
Our simulated average peat decomposition rate of 399 g C m −2 yr −1 during the period 1951 to 2011 is lower than the value measured in 2008, which was 630 g C m −2 yr −1 (Meyer et al., 2013).However, this high peat decomposition rate could be attributed to an inter-annual weather variation, which is corroborated by the high plant growth measured in 2008, 830 (±390) g C m −2 yr −1 .Our simulated N 2 O emission, 0.52 (±0.1) g N m −2 yr −1 during 2007 to 2009 is similar to the observed data collected these years, 0.71 (±0.59) g N m −2 yr −1 and measurements 2006 to 2011, 0.38 (±0.12) g N m −2 yr −1 (Holz et al., 2015).Only during these years, our predicted level of emissions was 0.50 (±0.12) g N m −2 yr −1 .Our simulated CO 2 and N 2 O fluxes are therefore generally comparable with the measured data.
Our simulated peat decomposition and N 2 O emissions are generally comparable in size with measured flux data from afforested drained peatland published in the literature (Table 3).However, when compared with the IPCC EF's for temperate drained nutrient-rich forest soil, which are given as 260 (200 to 330) g C m −2 yr −1 for CO 2 and 0.28 (−0.06 to 0.61) g N m −2 yr −1 for N 2 O (IPCC, 2014), our simulated values were found to be larger.This could be explained by the higher soil fertility at the Skogaryd site and also a deeper GWL (mean of 0.52 m during the simulated 60 years), compared to what pertained at those sites used for constructing the IPCC EF's.That the GWL is of crucial importance for emission levels for drained peat soils has also been shown by Couwenberg et al. (2011) and Leppelt et al. (2014).This could justify our assumption that our somewhat high estimates were due to deep and long-lasting drainage.
The unexpectedly low simulated N 2 O emission in the first years after planting was mainly caused by a high N uptake by the understorey vegetation, probably dominated by grasses, making less N available for nitrification and denitrification.From 1960 the understorey vegetation decreases and in 1966 the spruce trees are more dominant but still rather small, thus the model shows total N plant uptake to be smaller than in both the earlier and later periods and more inorganic N available in the soil (Supplement Fig. S1).During this period, the GWL was also deep (Fig. 5), resulting in high peat mineralisation rates and a net release of inorganic N.Both vegetation and GWL could thus explain the high N 2 O emission during the period 1966 to 1988.Besides, a high soil N availability also accounts for the "vegetation fitted" model consistently overestimating spruce growth during this period (Fig. 2c).However, our model also predicts some extreme annual N 2 O emissions, i.e. in 1969O emissions, i.e. in , 1973O emissions, i.e. in and 1977 (Fig. (Fig. 5e), modelled to occur simultaneously with large precipitation events lasting several days in summer.Since the model shows N 2 O to be produced close to the GWL at 50-70 cm depth and the model simply assumes the N gases to be directly emitted into the atmosphere after production thus no further reduction into N 2 is simulated.If a vertical gas transport process and further processing during transport was considered as performed by Stolk et al. (2011), this could further have converted N 2 O into N 2 thus simulating lower emissions.This is corroborated by gas profile measurements in Skogaryd where the soil N 2 O gas concentration increases with soil depth (He et al., 2016).If we remove the extreme annual emissions (1969, 1973 and 1977) in our calculations, the annual average N 2 O emission would change from 0.7 to 0.5 g N m −2 (thus 30 % lower).The accumulated CO 2 and N 2 O fluxes (NEE + N 2 O in Fig. 6) would in 2012 be 1000 g CO 2 equ m −2 instead of 7000 g CO 2 equ m −2 .However the forest ecosystem would still be a GHG source to the atmosphere.

Challenges of modelling long-term dynamics of an organic soil
Overall our modelling application indicates, given a few assumptions, that the CoupModel is generally able to simulate the decadal-scale dynamics of the drained organic soils used for forestry.However, our modelling exercise also reveals that there are some issues which still need to be more explicitly accounted for when simulating organic soils and which require further model development.These are the nature of the soil organic matter and physical changes of a peat soil.

A need for explicitly specifying the nature of soil organic matter
A multi-pool approach was developed for modelling SOM dynamics from mineral soils and has been shown to work well for forest mineral soils, e.g.(Svensson et al., 2008;Wu et al., 2013).However, for organic soils, because there is no explicit peat pool in the model, we have had to assume the peat to comprise an unknown mixture of the fast and the slow pool.In the present study we have assumed the initial values of SOM as only representative of the slow pool.The decomposition coefficients for the fast and slow pool were obtained by calibrating the model coefficient against the measured fluxes as we did in our previous study (He et al., 2016).However in this long-term simulation there is a continuous addition of spruce litter leading to resistant soil organic matter and a change in substrate quality over the simulation period for the slow pool.Although most existing models do not explicitly specify the nature of the organic matter (Smith et al., 1997), they can still simulate the total organic matter dynamics fairly well over a relatively short period.Metzger et al. (2015) found that the CoupModel could capture major C fluxes and the ecosystem dynamics when applied to five European treeless peatlands, where they pointed out that the total C flux was mainly determined by the decomposition coefficients of the total SOM.Continuous addition of organic matter into the slow pool from litter decomposition must also change the decomposition coefficient for the slow pool over time.However, this is seldom accounted for.In order to understand the long-term dynamics of organic matter, which might differ in origin and components, a more precise consideration of the changes of soil organic matter characteristics would be helpful.

Modelling physical changes of peat soil
For mineral soils in which the physical structure of the soil does not normally change over time, the CoupModel works well by assuming the soil layer profile to be fixed over time (Jansson, 2012;Jansson and Karlberg, 2011).However, this is not the case for organic soils where the soil structure is mainly built by soil organic matter, which gradually disappears through decomposition.Thus the soil's physical char-acteristics change over time, e.g. the pore structure, which could change the soil hydraulic conductivity and preferential flows (Kechavarzi et al., 2010).Moreover, decomposition makes the top soil disappear during a forest rotation, resulting in surface subsidence (Minkkinen and Laine, 1998;Leifeld et al., 2011;Hooijer et al., 2012).This causes the GWL to come closer to the soil surface, which in the normal case requires further drainage or ditch management.This process has not so far been implemented in the CoupModel, which currently is not able to account for surface subsidence, mainly due to lack of feedback coupling between the soil's biological and physical properties in the model.The model physical subroutine simulates the water and heat flow and then links this to the biochemical processes by response functions of water moisture and soil temperature.While there is no feedback to the soil physical processes arising from organic matter decomposition or other changes of the soil.All these processes remain a major challenge when applying the CoupModel to the long-term dynamics of a forest ecosystem on drained peatland.To quantify the uncertainty from surface subsidence, in the present study, the system was simplified by assuming a fixed drainage depth, whereas a range of values was used to quantify the model's sensitivity.The variation of the drainage depth had a considerable impact on the soil peat decomposition, as shown by the model sensitivity analysis (Table 2), which in turn highlights the need, when developing future models, to explicitly account for these processes when performing long-term simulations.

Initial soil C, N and soil C / N ratio
A major difficulty in the simulation was the unknown initial soil conditions.We chose to use the EF's 260 (200 to 330) g C m −2 yr −1 for CO 2 from the IPCC wetland supplement (IPCC, 2014), which compiles up-to-date observational data from similar sites under temperate climate conditions.Another alternative could be to use the subsidence rate to calculate the soil C losses, which has been applied in other published studies e.g.(Leifeld et al., 2011;Hommeltenberg et al., 2014).By taking the measured subsidence, 0.22 m (ranging from −0.15 to 1.03 m) during ca.60-year post-drainage period for Finnish drained afforested fens (Minkkinen and Laine, 1998), analogizing the measured total soil C in the upper 0.5 m in 2007, which was 55.3 kg C m −2 (Meyer et al., 2013), the estimated soil losses during the 60-year period would be 24.3 kg C m −2 , which is equivalent to a loss of 405 g C m −2 yr −1 , close to current modelling estimates, 399 g C m −2 yr −1 .Increased initial soil C in our sensitivity analysis shows both peat decomposition and plant growth to increase (Table 2).Compared to the "vegetation fitted" model, the combination of a small initial soil C, a large soil C / N ratio, and a shallow drainage, gives a larger reduction in plant growth than in peat decomposition, which is why the overall emissions of GHG increase.

GHG balance for the forest ecosystem
Our modelling indicates forest on drained agricultural peatland to be a strong net CO 2 source for the first 39 years of the forest rotation which changes into a CO 2 sink thereafter due to strong tree growth (Fig. 6).This means that, despite soil decomposition being high, the high growth rate of forest over 60 years compensates for most C losses.Meyer et al. (2013) also showed the forest ecosystem in Skogaryd to be an overall GHG sink (410 g CO 2 eq ha −1 m −2 ) in 2008, a year when the plant growth rate was at its maximum, thus offsetting the high rate of peat decomposition.Our findings are also generally in line with the few previous field investigations conducted on afforested drained agricultural peatlands where Mäkiranta et al. (2007) and Lohila et al. (2007) found a 30-year old Scots pine forest on drained agricultural bog to be, overall, a small source of CO 2 (50 g C m −2 yr −1 ), which was explained by a small leaf area index (varying between 0.7 and 2 during the observational period).Another study by Hommeltenberg et al. (2014) reported an afforested drained bog in Germany, previously used for agriculture, to emit 500 g C m −2 yr −1 .By combining eddy covariance measurements and biometric estimation, they concluded it to be a major CO 2 source, emitting a total of 13.4 kg C m −2 over the last 44 years.However, their short-term measurements (2010 to 2012) also indicated that forest growth offsets peat decomposition, a result similar to our study.
Growing forests on drained peat is done at the cost of the soil peat, which has generally accumulated slowly during the last millennia (the last 4000 years in Skogaryd).When the forest growth has been larger than the soil loss, the system has been interpreted as being an overall sink (Meyer et al., 2013;Hommeltenberg et al., 2014).However, the soil loss and the forest gain can be viewed as a "relocation" of the peat carbon into timber carbon.The simulated NEE (Fig. 6) shows that the system remains a sink for 2 decades but growth rate probably declines over time, as shown in the simulated period from 2011 to 2031.After 2031, to continue keeping forest on these lands, a lower growth (Luyssaert et al., 2008) will further decline the sink strength while the peat soil will continue to decompose as long as it is kept aerated by a living transpiring forest.Sudden fires would also be a risk releasing the forest biomass C.However, the forest in Skogaryd is not a nature reserve but a managed forest already mature for harvesting, commonly done at 80 years of age in southern Sweden.The harvested wood products over a forest rotation are used for both timber and paper, about 40 and 60 % (Sweden CRF table 4.Gs2 for year 2013, submitted to the UNFCC 2015) having a half-life of 30 and 2 years respectively (IPCC, 2006).Thus the carbon will soon be released as CO 2 .The best alternative would be the use timber for wooden buildings which otherwise should have been built by using concrete (Gustavsson et al., 2006).The displacement of concrete by wood could according to a meta-analysis by Sathre and O'Connor (2010) avoid emissions by 2.1 times the C content of the timber.However, even then, most buildings do not last more than a century and only a few buildings are functional for longer periods.Thus, most harvested biomass will soon be burnt releasing the stored C.These indirect emissions following the consumption of wood would shift the system from an overall small sink into a large GHG source (Fig. 6).Another alternative use of the biomass could be as biochar in agricultural soils (Ojanen et al., 2013), which potentially could shift the system into an overall GHG sink.However, this is equivalent to releasing C from peat and storing it in agricultural soils, and it is not clear for how long the char carbon persists.Additionally, there are some other direct and indirect GHG sources that become apparent during the full forest rotation period which we have not accounted for, such as methane emissions in drainage ditches and loss of dissolved organic C or particulate organic C.However, these contributions to the overall GHG balance are in general of minor importance and thus not likely to alter the overall picture (Meyer et al., 2013;Hommeltenberg et al., 2014).In summary, the overall message is that a forest rotation on fertile drained peat soil has a long-term GHG cost, never reaching a balance, and thus the wood products produced on peat soil cannot be regarded as renewable products.
In Sweden, forests on drained peatland cover 1.7 Mha (Maljanen et al., 2010;Von Arnold et al., 2005a), of which 0.4 Mha has high fertility, comparable to the soil in the present study.According to our simulations, these forests emit around 1.74 kg CO 2 eq m −2 yr −1 (peat decomposition and N 2 O emissions).Thus these fertile drained peat soils in Sweden emit 7 Mtonnes CO 2 eq annually, which is equivalent to 12 % of the emissions coming from all other sectors in Sweden when excluding LULUCF (land use, landuse change and forestry).From a climate change perspective, forested drained peatlands should be highlighted for actions, especially following forest clear-cut.Instead of digging the ditches deeper for replanting a new forest, making the soil wetter would reduce the soil decomposition, as shown by our sensitivity analysis and other studies (e.g.Karki et al., 2014).However, these measures need support from policy makers since landowners often only recognize revenues from forest production, not the cost of GHG emissions.

Conclusions
Our simulation study shows that the GHG fluxes in a forested drained peatland are composed of two important quantities: C uptake by forest growth, and C losses from the soil.By fitting the CoupModel to the spruce growth, up-scaled from radial tree-growth observations, we obtained a "vegetation fitted" model by which we were able to describe the C and N fluxes over 60 years.We show that the forest C growth is tightly coupled to soil C losses, and if the forest is harvested and used, there will only be losses over time.

Figure 1 .
Figure 1.Conceptual representation of the dynamics of plants and peat soil development over a forest rotation period.The upper panels (a, b, c, d) represent the conceived reality and (e) and (f) represent the CoupModel conceptualization.For all the figures, spruce tree and understorey vegetation, e.g.grasses are considered but for clarity, understorey vegetation is only shown in (a)."C 2007" in (f) represents the measured total soil C in the upper 0.5 m of the soil profile in 2007, and "C 1951" is the total soil C in the upper 1 m of the soil profile, as back-calculated from the equation: 2 × "C 2007" + (2007-1951) × IPCC EF's.Any variation of climate during the forest development in this conceptual figure is not considered.

Figure 2 .
Figure 2. (a) Simulated (black line) spruce absorbed radiation; (b) simulated and measured (red hollow circle) leaf area index; (c) annual spruce tree growth rate; (d) total spruce tree biomass; (e) spruce tree biomass for different components.In (e) the solid red symbols show the calculated plant biomass of leaf biomass, root and stem biomass using the allometric function given by Meyer et al. (2013).

Figure 3 .
Figure 3. Simulated development of major soil C pools in the first metre of soil, from 1951 to 2011.The red circle shows the measured total soil C in 2007 (±95 % confidence intervals) by Meyer et al. (2013).

Figure 4 .
Figure 4. (a) Simulated (black line) and measured (red hollow circle) total net radiation; (b) soil surface temperature (0-5 cm depth; c) GWL; (d) NEE.Measured data used to create these plots are 5day averages, except for NEE where daily averages have been used.

Figure 5 .
Figure 5.For the period 1951 to 2011: (a) annual precipitation (mm yr −1 ) and air temperature ( 0 C); (b) the simulated annual NPP of spruce trees (g C m −2 yr −1 ); (c) simulated annual peat decomposition rate (g C m −2 yr −1 ); (d) simulated annual NEE (g C m −2 yr −1 ); (e) simulated annual N 2 O emissions (g N m −2 yr −1 ); (f) simulated annual GWL (m).The dashed reference line separates the duplicated 1951 to 1961 and real climate 1961 to 2011.The source or sink is based on the atmospheric perspective, e.g. the soil emissions are sources, and plant uptakes are sinks.

Figure 6 .
Figure 6.Simulated total GHG balance for the forest ecosystem from 1951 to 2011 and extended to 2031.The simulated results of 2011 to 2031 are obtained by running the "vegetation fitted model" with meteorological data from 1991 to 2011 extended to represent the climate of 2011 to 2031.It should be noted that the GHG balance presented in this figure assumes no final harvest.

Table 1 .
Soil C content in the soil profile during 1951 to 2011 estimated by the vegetation fitted model, kg C m −2 .
* Negative change means an increase of soil C.

Table 2 .
Model sensitivity: change compared with "vegetation fitted" model during 1951 to 2011.

Table 3 .
A comparison of soil peat CO 2 and N 2 O emissions in the present study with values published in the literature.−2 yr −1 ) (g N m −2 yr −1 ) Calculated by assuming 50 % of measured soil respiration to have originated from root-based activity.

www.biogeosciences.net/13/2305/2016/ Biogeosciences, 13, 2305-2318, 2016 range
The model sensitivity analysis conducted provides evidence that a wide of drainage depths, site fertilities and initial soil C contents lead to similar overall results.Further model developments are however needed to better simulate the drained peat soil over a forest rotation period.

Supplement related to this article is available online at doi:10.5194/bg-13-2305-2016-supplement.
Author contributions.Hongxing He, Åsa Kasimir, Per-Erik Jansson and Leif Klemedtsson planned and initialized the study.Hongxing He conducted the data analysis and modelling under supervision from Åsa Kasimir, Magnus Svensson and Per-Erik Jansson.Jesper Björklund and Lasse Tarvainen helped Hongxing He with the tree ring data collection and analysis.Hongxing He and Åsa Kasimir wrote the paper with all authors commenting and participating in the interpretation of the results and contributing to the discussions.