Biogeosciences A dynamic model of wetland extent and peat accumulation : results for the Holocene

Substantial deposits of peat have accumulated since the last glacial. Since peat accumulation rates are rather low, this process was previously neglected in carbon cycle models. For assessments of the global carbon cycle on millennial or even longer timescales, though, the carbon storage in peat cannot be neglected any more. We have therefore developed a dynamic model of wetland extent and peat accumulation in order to assess the influence of peat accumulation on the global carbon cycle. The model is based on the dynamic global vegetation model LPJ and consists of a wetland module and routines describing the accumulation and decay of peat. The wetland module, based on the TOPMODEL approach, dynamically determines inundated area and water table, which change depending on climate. Not all temporarily inundated areas accumulate peat, though, but peat accumulates in permanently inundated areas with rather stable water table position. Peatland area therefore is highly uncertain, and we perform sensitivity experiments to cover the uncertainty range for peatland extent. The peat module describes oxic and anoxic decomposition of organic matter in the acrotelm, i.e., the part of the peat column above the permanent water table, as well as anoxic decomposition in the catotelm, the peat below the summer minimum water table. We apply the model to the period of the last 8000 years, during which the model accumulates 330 PgC as catotelm peat in the peatland areas north of 40° N, with an uncertainty range from 240 PgC to 490 PgC. This falls well within the range of published estimates for the total peat storage in high northern latitudes, considering the fact that these usually cover the total carbon accumulated, not just the last 8000 years we considered in our model experiments. In the model, peat primarily accumulates in Scandinavia and eastern Canada, though eastern Europe and north-western Russia also show substantial accumulation. Modelled wetland distribution is biased towards Eurasia, where inundated area is overestimated, while it is underestimated in North America. Latitudinal sums compare favourably to measurements, though, implying that total areas, as well as climatic conditions in these areas, are captured reasonably, though the exact positions of peatlands are not modelled well. Since modelling the initiation of peatland growth requires a knowldge of topography below peat deposits, the temporal development of peatlands is not modelled explicitly, therefore overestimating peatland extent during the earlier part of our experiments. Overall our results highlight the substantial amounts of carbon taken up by peatlands during the last 8000 years. This uptake would have substantial impacts on the global carbon cycle and therefore cannot be neglected.


Introduction
Estimates of the amount of carbon stored in boreal peatlands vary.Yu et al. (2010) estimate that peat deposits of about 547 PgC have accumulated in the northern high latitudes during the Holocene, though other estimates are substantially lower, e.g.273 PgC estimated by Turunen et al. (2002).Nonetheless it is clear that boreal peatlands store substantial amounts of carbon, which may be up to a fifth of the total global soil carbon estimated as 2344 PgC in the top three meters (Jobbagy and Jackson, 2000).On interannual timescales, the changes in peat storage are rather small, and peat accumulation has therefore been neglected in carbon cycle models so far.On millennial timescales, this is a substantial factor in the carbon cycle, though, which is why we have developed a dynamical model of wetland extend and peat accumulation as described in this paper.
In order to represent the carbon cycle forcing by peatlands on millennial timescales, previous authors used scenarios of peat accumulation (Wang et al., 2009;Kleinen et al., 2010).While the use of such scenarios is possible for the Holocene, where some measurement data to derive the scenarios exists, Published by Copernicus Publications on behalf of the European Geosciences Union.such data do not exist for previous interglacials, since the glaciation occurred in just those places, where present-day peat deposits are located.The investigation of carbon cycle dynamics in previous interglacials therefore requires the use of a dynamical model.
Models of peat accumulation have previously mainly been developed for single sites.Clymo (1984) developed a onedimensional model of peat accumulation and decay in a single peat column.Later Clymo et al. (1998) determined parameter values for this model by fitting to profiles from numerous peat bogs from Finland and Canada.Contrary to the Clymo (1984) model which focuses on the biochemical decomposition processes in the peat layers below the water table, Ingram (1982) has developed a model of peat from a hydrological perspective, describing the groundwater table in a domed peat deposit.Over time these models have included more and more processes, for example threedimensional bog growth, i.e., the lateral expansion of peatlands from the site of first initiation (Korhola et al., 1996), more sophisticated parameterisations of water table depth, or more plant functional types (PFTs) (Frolking et al., 2010).This line of development has so far culminated in the Frolking et al. (2010) model describing the development of annual peat layers and thereby resolving the accumulation and decomposition processes in considerable detail.
The other approach is to model peat accumulation and decomposition in the framework of dynamic global vegetation models (DGVM) or the land surface components of climate models.Such a global approach necessarily neglects some of the detail included in dedicated site models, due to computational constraints, but also due to the fact that detailed parameterisations often are dependent on site specific parameters.Examples of models following this global approach are Wania et al. (2009a,b) who developed an extension of the LPJ DGVM, accounting for organic soils with the rationale to derive methane emissions from wetlands.Their model relies on prescribed wetland areas and does not determine the extent of wetlands dynamically.Ringeval et al. (2010) similarly modeled wetland CH 4 emissions using maps of wetland extent.
Parameterisations for the determination of wetland extent have been developed and implemented in a number of cases.Kaplan (2002) used a digital elevation model (DEM) to determine areas, where wetlands could develop.Using an approach based on TOPMODEL (Beven and Kirkby, 1979), wetland parameterisations have been developed for the NCAR GCM (Niu et al., 2005), the ISBA land surface model ORCHIDEE (Habets and Saulnier, 2001;Decharme and Douville, 2006), and the MetOffice land surface scheme MOSES (Gedney and Cox, 2003;Gedney et al., 2004).In all of these cases an explicit accounting for the long term accumulation and decay of peat is missing, since these studies focus on methane emissions, not peat accumulation.
The present study therefore aims to combine these approaches, implementing a model for peat accumulation and decay, as well as a parameterisation for wetland extent in order to assess the carbon accumulation in peat over the Holocene.Peatland development is a process highly dependent on specific local conditions, and it would therefore be impossible for a model of relatively coarse resolution and global scale to exactly simulate the development of every single peatland.Our model has the much more moderate aim to assess the carbon uptake and storage by the land surface in peatland systems.The aim therefore is to model average peatlands as they would develop under the average topographic and climatic conditions in a grid location, while neglecting some of the heterogeneity of site conditions.

CLIMBER2-LPJ
In the present study we are using the model CLIMBER2-LPJ as described in Kleinen et al. (2010).CLIMBER2-LPJ consists of the earth system model of intermediate complexity (EMIC) CLIMBER2, coupled to the DGVM LPJ.This combination of models allows experiments on timescales of an interglacial due to the low computational cost of CLIMBER2, while accounting for the heterogeneity of land surface processes on the much more highly resolved grid of LPJ.
CLIMBER2 (Petoukhov et al., 2000;Ganopolski et al., 2001) is an EMIC consisting of a 2.5-dimensional statisticaldynamical atmosphere with a latitudinal resolution of 10°a nd a longitudinal resolution of roughly 51°, an ocean model resolving three zonally averaged ocean basins with a latitudinal resolution of 2.5°, a sea ice model, and a dynamic terrestrial vegetation model (Brovkin et al., 2002).In the present model experiments, the latter model is used only for determining biogeophysical responses to climate change, while biogeochemical effects, i.e. the corresponding carbon fluxes, are determined by LPJ.
To this EMIC we have coupled the DGVM LPJ (Sitch et al., 2003;Gerten et al., 2004) in order to investigate land surface processes at a resolution significantly higher than the resolution of CLIMBER2.We also implemented carbon isotope fractionation according to Scholze et al. (2003).
LPJ is run on an 0.5°× 0.5°grid and is called at the end of every model year simulated by CLIMBER2.Monthly anomalies from the climatology of the temperature, precipitation and cloudiness fields are passed to LPJ, where they added to background climate patterns based on the Climatic Research Unit CRU-TS climate data set (New et al., 2000).In order to retain some temporal variability in these climate fields, the anomalies are not added to the climatology of the CRU data set, but rather to the climate data for one year randomly drawn from the range 1901-1930.LPJ simulates the changes in carbon pools and the carbon flux F AL between atmosphere and land surface, which is passed back to CLIMBER2 and is employed to determine the atmospheric CO 2 concentration for the next model year.Biogeochemical feedbacks between atmosphere and land surface are thus determined by the combination of CLIMBER2 and LPJ, while biogeophysical effects are solely determined by the CLIMBER2 land surface model, which includes its own dynamical vegetation model.

Modelling peat accumulation
A natural peatland consists of three functionally distinct layers.At the top there is a live plant layer, where plants generate organic matter through photosynthesis.Below that is an upper layer of peat, which usually is less than 30 cm in height (Belyea and Baird, 2006), the so-called acrotelm, located above the permanent water table and therefore aerobic for at least part of the year.At the bottom is the peat located below the permanent water table, the so-called catotelm.This latter layer can be several meters in height, and significant amounts of carbon can therefore be stored in the peat column.
During peat formation, the process can be described as follows.The live plants at the surface generate organic matter through photosynthesis.Dead organic matter is added to a litter layer, from which it passes to the acrotelm very quickly.In the acrotelm the organic mater is decomposed, either aerobically or anaerobically, depending on the position of the water table, and once decomposition has proceeded far enough, the organic matter suffers a structural collapse, which significantly enhances density while shrinking pore volume.The water is squeezed out of the organic matter and the permanent water table rises slightly, adding some more organic material to the catotelm (Belyea and Baird, 2006).
In principle it is possible to model this process of peat formation explicitly, as Frolking et al. (2010) have shown.Modelling the change in density of the organic matter requires keeping track of annual cohorts of organic material in order to model how they pass through the peat column and to determine how density changes in each peat layer.This approach therefore requires substantial amounts of computer memory if implemented on a global or even hemispheric scale.We therefore decided against this approach but rather approximate the peat formation process by assuming a catotelm formation rate which is proportional to the amount of carbon in the acrotelm, modified by annual mean temperature, an approach very similar to how soil organic matter decomposition is modelled in LPJ.For rates of peat accumulation, as well as the distribution of peatlands and total peat storage, numerous estimates exist for boreal regions, while far fewer estimates exist for tropical peatlands.Since these estimates are essential for calibrating and validating the model, we currently limit our investigation to the regions north of 40°N.

Peat dynamics
If we compare the soil carbon dynamics in wetland and nonwetland soils, the main difference is that part of the soil column in wetland soils is below the water table, which leads to anaerobic decomposition of soil organic matter.LPJ contains a number of carbon pools, as shown in Fig. 1, on the left.There are live biomass pools for carbon (C) in leaves, wood and roots, here shown as a single pool C B for simplicity.Then there are pools for carbon stored in litter, here summarised as C L , and finally there are pools C S for carbon stored in soil.
For the peat version of LPJ, this pool structure needs to be extended by an additional belowground C pool containing the carbon in the catotelm that is decomposing under anoxic conditions all year round, shown in Fig. 1, right hand side, as C C .
In line with this structure of C pools, a vertical column of peat would have to be seen as shown in Fig. 1, on the right.At the surface there is a litter layer.Below that is the acrotelm, a soil layer where carbon is decomposed partly under oxic and partly under anoxic conditions, depending on the position of the water table.Finally, below the minimum water table position, there is a soil layer where anoxic decomposition occurs all year round, the catotelm.
These C pools can then be modelled as follows: Here, the F XY are the carbon fluxes between the C pools, and the R k are the respiration fluxes to the atmosphere.XY can have the meanings BL = biomass-litter, LA = litteracrotelm, and AC = acrotelm-catotelm, while k can be L = litter, A = acrotelm (with o for oxic and a for anoxic conditions), and C = catotelm.Leaching of dissolved organic carbon (DOC) is not considered explicitly, but rather assumed to be implicitly contained in the respiration fluxes.
In order to keep the peat version of LPJ as close to the original as possible, we keep flux formulations for existing carbon pools as they are in the original version, but the fluxes F AC and R C have to be added.The fluxes F XY are dependent on the carbon mass (or rather area density) C X in the originating C pool X, as well as a temperature (and moisture) dependent decay function, and the same holds for the respiration fluxes R X .In concrete terms this is: with the rate constants k i modified multiplicatively by a response function γ (T soil ,w Soil ) depending on soil temperature T soil , as well as soil moisture w Soil .Since we are considering wetland processes, we assume that moisture is not a limiting factor to decomposition and therefore don't consider a dependence on w Soil .The temperature dependence of decomposition is quite often modelled as an exponential function exp(ln(Q 10 )(T − T Ref )/T Ref ), and measured Q 10 factors vary widely.Scanlon and Moore (2000), for example, measured Q 10 values ranging from 1.0 to 7.7 for peat decomposition.For ecosystem respiration, on the other hand, recent results indicate that Q 10 is 1.4, despite the huge range of measured Q 10 determined in laboratory studies (Mahecha et al., 2010).LPJ generally uses the formulation by Lloyd and Taylor (1994) for temperature dependence, which gives a temperature dependence roughly similar to a Q 10 of 2. We also apply it in wetland systems.Finally, α determines the fraction of decomposed litter that is added to the soil, while the remainder is respired.This is not changed from the nonwetland version of LPJ, either.
The case of the acrotelm respiration rate R A is slightly more complicated, since acrotelm peat above the water table decomposes aerobically, while it decomposes anaerobically below.In Eq. 7 for the acrotelm respiration under oxic conditions, β is the fraction of the acrotelm above the water table, which decomposes at the rate k A , while the rest (Eq.8) decomposes anaerobically at the rate υk A , with υ the ratio of anaerobic to aerobic CO 2 production.We follow Wania et al. (2009b) in setting this to 0.35.Finally, catotelm formation is modelled similar to decomposition, with the formation flux F AC depending on the peat formation rate constant k P , while decomposition depends on a decomposition constant k C .Clymo et al. (1998) determined this rate constant by fitting peat core data to a similar decomposition model, and their value translates to k C = 3.35×10 −5 a −1 if corrected for mean annual temperature.The peat formation rate constant k P , as well as the acrotelm respiration rate k A , we determined by comparing model results to measured acrotelm mass (Malmer and Wallen, 1993) and peat formation rates (Yu et al., 2010).
The fraction β of the acrotelm above the water table is determined by comparing acrotelm height, calculated from acrotelm density ρ A , acrotelm carbon fraction c f,A and acrotelm mass density C A , and water table position w w .The latter is assumed to be relative to the acrotelm-catotelm interface, which is located at the 50 year mean summer minimum water table position.The acrotelm density ρ A was taken from Granberg et al. (1999), who gives a density for the surface peat and for the peat at the bottom of the acrotelm.We therefore assume the mean of these values to be the acrotelm density.All parameter values used are listed in Table 1.
With regard to biomass input into the peat column, we decided not to implement special wetland PFTs in the interest of keeping the model as close to the original as possible.Initial tests showed that the productivity of (modelled) mosses is very similar to grasses, making an additional PFT unnecessary.We did follow Wania et al. (2009b), though, in introducing their parameterisation for inundation stress in trees since tree growth is strongly inhibited in wetlands.

Dynamic wetland model
While the resolution of LPJ at 0.5°already is rather high in comparison to typical resolutions of climate model land surface schemes, this still translates to a grid cell size of 50 km × 30 km at 60°N.Since most wetlands are of smaller extent than this, an approach is required that determines the grid cell fraction covered by wetlands.Since it is our aim to apply the model to previous interglacials, as well as times earlier in the Holocene, using a simple wetland map to determine grid cell wetland area based on present day observations is insufficient.Instead, a scheme to determine wetland extent dynamically is required.
For this purpose we have implemented an approach based on the TOPMODEL hydrological framework (Beven and Kirkby, 1979).TOPMODEL is a conceptual rainfall-runoff model that is designed to work at the scale of large watersheds using the statistics of topography, instead of requiring detailed topographic information.It is based on the compound topographic index (CTI) χ i = ln(α i /tan(β i )) with α i a dimensionless index for the area draining through point i

Biogeosciences
with χ i the local CTI index in point i, χ the grid cell mean CTI index, and f a parameter describing the exponential decline of transmissivity with depth for each soil type.Using this equation we can determine the grid cell fraction that is inundated, i.e., with a water table at or above the surface, as well as the mean water table height in the inundated fraction.The inundated area consists of all points within the grid cell that have a local water table depth z i ≥ 0, but since a local water table that is well above the surface implies either running water, i.e., a river, or a lake, we also set a maximum CTI value χ max , similar to Gedney and Cox (2003), which depends on soil type but otherwise is constant in space and time.We therefore assume the grid cell area with 0 ≤ z i ≤ z max , z max = z(χ max ) to be the grid cell wetland area.Finally, the wetland water table position w w is determined from Eq. 10, using the mean CTI index of the grid cell wetland fraction.Values for the parameters f and χ max were determined separately for each soil class in the LPJ soil map, with the exception of organic soils.These were assigned the values for medium coarse soils in order to prevent prescribing wetland location.Parameter values are listed in Table 2.
In order to take into account the modification of soil infiltration properties by permafrost, we follow Fan and Miguez-Macho (2011) who modify f multiplicatively by a function k depending on January temperature T Jan .In our case, a suitable modification seems to be: In order to determine the grid cell mean water table, we slightly modified the Stieglitz et al. (1997) approach to a formulation appropriate for LPJ: with z b the bottom of the soil column, z the height of the soil column, w the soil column average soil moisture, and w thres the minimum soil moisture for a water table to form.This modification of the original approach became necessary, since soil moisture in LPJ is a variable determining the plant available water as a fraction of field capacity, i.e., w = 0 at the wilting point and w = 1 at field capacity.Similarly, since the LPJ soil column is very shallow having only 2 layers and a total soil column height of 1.5 m, we are simply using the soil column average soil moisture w instead of the layer soil moistures.
Sensitivity experiments showed that a minimum soil moisture w thres = 0.1 yielded good results.In addition, an initial comparison with Lehner and Döll (2004) wetland area showed that wetland area in grid cells with a mean CTI index χ ≤ 5.5 is negligible.Wetlands are therefore only determined for grid cells with χ > 5.5.
The scheme described above determines the monthly grid cell fraction that is inundated.For peat to develop, areas that are just inundated for a few days during the course of the year are not relevant, but rather the areas that are inundated permanently, or at least during the growing season.Modelled soil water dynamics during the winter season, on the other hand, cannot be trusted in high latitude areas since LPJ does not consider permafrost, leading to too low water table positions during periods when the ground is frozen.We therefore use the inundated area and water  extent.The scheme to determine wetland extent described in this section is a dynamic scheme, i.e., wetland extent can change as the climatic conditions and therefore the soil moisture change.In order to limit the interannual fluctuations in peatland extent, we are using a 50 year running mean peatland extent in all calculations.As peatland extent changes, carbon pools have to be updated.For these transfers of carbon between the wetland and the non-wetland part of the grid cell, we are making two assumptions.First we assume the peat deposit to have a parabolic overall shape, i.e., carbon storage from the edge to the centre of the peat deposit follows a parabola, and second we assume that a wetland that shrunk previously and then expands again expands into the same area it covered previously.
In case wetland extent shrinks, some of the carbon stored in the catotelm pool C C will have to be passed into the soil pool of the non-wetland part of the grid cell.The fraction of C C to be passed is determined from the proportional change in wetland size assuming that the peat deposit has an overall parabolic shape with peat removed from the outer edge.The carbon is then transferred into the soil carbon pool of the non-wetland part of the grid cell, but a record is kept of previous wetland extent and amount of carbon transferred from the catotelm pool in case the wetland should grow again.In case of a growth in wetland extent, carbon is proportionally transferred from the soil carbon pool to the oxic wetland carbon pool.In case a wetland expands into areas previously covered by a wetland, carbon is also added to the anoxic pool C C .

Model experiments
We have performed two sets of model experiments.One is a model experiment under constant preindustrial boundary conditions, while the other experiment is a transient fully coupled model run with interactive CO 2 for the last 8000 years, from 8 ka BP to preindustrial.Both experiments were initialised from a 2000 year spinup under the corresponding boundary conditions.For each experiment, we performed three different model runs at slightly changed parameter set-tings in order to determine an uncertainty range for the model result.Since we estimate the peatland extent to be the most important uncertain parameter in the model, we vary the peatland area by setting it to the summer (JJAS) mean, minimum and maximum wetland extent as elaborated in Sect.3.1.We only consider peat accumulation in the high northern latitudes, i.e., north of 40°N, tropical peatlands are not considered at the moment.

Glacial extent during the early Holocene
When assessing the carbon uptake by peatlands during the Holocene, one further factor needs to be considered.At 8 ka BP, the start of our experiments, some remnants of the Laurentide ice sheet still existed.Furthermore, some areas in North America and Scandinavia had been depressed below sea level by the overlying ice mass during the glacial period.Our model, on the other hand, does not contain a dynamic model of the ice sheets at sufficiently high resolution, which would allow a consistent assessment of the effect of these changes.In order to obtain a first estimate of the consequences of the ice sheet presence, we used a scenario of ice sheet development and postglacial rebound determined with the ICE-5g model (Peltier, 2004).As a final step of the analysis, we use this scenario to determine the areas unsuited to peat formation and remove these areas from the assessment.

Wetland extent
For the global extent of peatlands, existing estimates suffer from shortcomings.Generally two types of estimates can be distinguished: Estimates based on peatland inventories and estimates based on remote sensing data.Estimates based on inventories suffer from different categorisations employed by different agencies.The map of global peatland areas shown by Yu et al. (2010), for example, marks the entire area of Sweden as a peatland, since it is reported as such in the national data the map is based on, which clearly is an overestimate.Remote sensing data, on the other hand, has the advantage of globally consistent coverage, but the disadvantage that so far no refined remote sensing datasets of peatlands exist (Krankina et al., 2008).Since the model primarily determines inundated area, we use datasets of wetland area as a proxy for peatland area.Here we are relying on two data sets.We used the global lakes and wetlands database (GLWD) (Lehner and Döll, 2004), which shows the annual maximum extent of wetlands, based on a combination of maps and remote sensing data.In addition, we applied the data set of remotely sensed surface water extent by Prigent et al. (2007) showing the monthly surface water extent from January 1993 to December 2000, which agrees reasonably well with GLWD for the maximum extent (Papa et al., 2010).For wetlands in the high latitudes it is very likely, though, that snow cover will mask wetlands during the snow season, making surface water extent measurements by remote sensing impossible during this time.In addition, passive satellite sensors rarely are able to sense inundation below a tree canopy, leading to an underestimate in forested areas.
In Fig. 2 we show zonal sums of wetland (inundated area) extent.Since snow cover makes the remote sensing of inundated area impossible, we limit the analysis to the extended summer season June to September (JJAS).For the total summer maximum inundated area north of 40°N, modelled extent is 3 % larger than the extent in Prigent et al. (2007), and 10 % larger than GLWD.The total maximum inundated area therefore is captured quite reasonably.The summer mean area is overestimated by 5 %, while the summer minimum area is overestimated by 30 %.If we consider the latitudinal variation of summer maximum extent in Fig. 2a, the model captures the latitudinal distribution reasonably well.For the summer mean extent in Fig. 2b, the model slightly underestimates the area around 55°N, while it overestimates the area north of 65°N.These biases are more pronounced for the summer minimum area shown in Fig. 2c.
Investigating the spatial distribution of wetlands more closely, the map in Fig. 3a shows the summer mean wetland extent as a grid cell fraction.The general locations of wetlands agree reasonably well with the map of peatland areas by Yu et al. (2010).In a direct comparison to the wetland extent by Prigent et al. (2007), shown in Fig. 3b, the model underestimates wetland extent in North America, while it over- estimates it in most parts of Eurasia, with the exception of southern Finland.For the major peatland areas like the Hudson's Bay lowlands and western Siberia, the model tends to slightly underestimate wetland extent in these areas, while it overestimates wetland extent elsewhere.While the latter could in principle be changed by adjusting the parameters f and χ max , the underestimation of wetland extent in regions with extensive areas of relatively flat terrain is a shortcoming of the TOPMODEL approach.It follows from Eq. 10 that the TOPMODEL approach redistributes the available water within a grid cell.Therefore the maximum grid cell fraction that can become a wetland is limited.In the underlying LPJ hydrology grid cell water content is limited to field capacity.The grid cell mean water table therefore never is above the surface, and even if the grid cell water table is near the surface, some fraction of the grid cell will always have a water table below the surface.
So far we have discussed the wetland extent, or rather the extent of inundated area.What fraction of this wetland extent becomes a peatland, on the other hand, is highly uncertain.For peat to form it is required that anoxic conditions prevail in the catotelm all year round, while oxic conditions may occur in the acrotelm for part of the year.The summer minimum wetland extent therefore almost certainly underestimates the peatland area, while the summer maximum extent usually occurs right after snowmelt when the soil is still partially frozen and prevents the infiltration of water.The summer maximum extent therefore almost certainly overestimates the peatland area.In order to consider this uncertainty in peatland extent, we perform a set of three sensitivity experiments.In these experiments, peatland area is set to the summer minimum, mean and maximum wetland extent.
Over the course of the Holocene, climate changed significantly.Due to the decrease in high latitude insolation through orbital changes, temperatures and therefore evapotranspiration decreased, while precipitation mostly increased in the mid to high latitudes.This would lead to wetter soils, a higher grid cell mean water table, and also a larger wetland extent.Model results for the extent of inundated area at 8 ka BP are shown in Fig. 4a, with changes during the last 8 ka shown in Fig. 4b.During the last 8 ka, wetland area increased slightly in northern Siberia and eastern Canada, while it slightly shrank in the Canadian Arctic.At the southern edge of the domain, inundated area shrank.As shown in Fig. 4c, some areas were located below the Laurentide ice sheet or below sea level at 8 ka BP.If ice sheets are taken into account, these areas would therefore show an increase in inundated area by 100 %.
Published estimates of peatland area change, for example MacDonald et al. (2006) and Yu et al. (2010), are based on basal dates of peatland initiation.These show a much larger increase of roughly 50 % in peatland area during the time considered.This discrepancy occurs because we estimate something different.The published estimates assess when peatlands were initiated and how this affected the total peatland area.Our figure, on the other hand, assesses how the present-day peatlands would change under the changed hydrological conditions of previous times.This implicitly assumes that topography is unchanged, an assumption not valid in cases where no peatland was present previously, since the accumulation of peat changes topography by filling in depressions and generally raising the land surface.In order to estimate the same measure as the reported changes in peatland area, it would be required to know the topography below the present day peatlands in order to explicitly model peatland initiation and subsequent topography change, which is presently impossible since topography below the peat layers in unknown for large areas.

Peat accumulation: The acrotelm
In the wetland areas the model accumulates carbon as peat, i.e., organic matter that decomposes very slowly since decomposition takes place under anaerobic conditions.Contrary to data on catotelm carbon, very few studies exist on the acrotelm.We are aware of only a single study comparing sites at different locations with a common methodology.In this study Malmer and Wallen (1993) investigate the acrotelm at 12 sites in Canada, Scotland, and Scandinavia.Malmer and Wallen measured carbon to nitrogen ratios in the peat cores.The ratio first decreases quickly as one goes deeper in the peat, due to the fact that carbon decomposes, while nitrogen is conserved.Further down in the core, it decreases more slowly, and Malmer and Wallen identified the intersect of the two C/N trends as the "decay decrease level", i.e., the transition from acrotelm to catotelm.In their study, they determine decomposition rates for the acrotelm and report both the height of the acrotelm, as well as the amount of organic matter contained therein.This dataset shows widely varying acrotelm heights and organic matter contents for very similar climatic conditions, indicating that the exact acrotelm thickness is very much dependent on local conditions the model cannot capture since grid spacing is 0.5°× 0.5°.Assessing the general acrotelm properties, the model develops an acrotelm layer with a mean thickness of 0.39 m, with a 5 % percentile of 0.09 m and a 95 % percentile of 0.66 m.This is within the generally assumed range (Charman, 2002), though slightly biased towards the thicker acrotelms reported for fens.

Peat accumulation: the catotelm
For the peat accumulated in the catotelm, numerous studies exist.They range in scope from studies of single sites, as in Yu (2006), for example, over regional summaries for Siberia (Kremenetski et al., 2003;Beilman et al., 2009) and North America (Gorham et al., 2003) to a recently published global summary of peat accumulation rates by Yu et al. (2010).While the single site publications usually are rather difficult to compare due to different conventions used, as well as different measures reported, the regional and global summaries employ standard methodologies and can therefore be compared rather well.While Yu et al. (2010) report catotelm accumulation rates, Gorham et al. (2003); Kremenetski et al. (2003) and Beilman et al. (2009) report the basal date and basal depth of the peat accumulated.From this we determined the long term apparent rate of carbon accumulation (LORCA) by dividing the basal depth by the basal date in years BP and converting the height accumulation rates into a carbon accumulation rate by using the C fraction and density of catotelm peat.If multiple measurements existed for a single LPJ grid cell, we compared to the average.This comparison is shown in Fig. 5, where the error bars for the model accumulation rate are determined from the three sensitivity experiments with differing peatland extents.Overall, the spread of measured values is substantially larger than for the model results, reflecting the influence of local conditions on peat accumulation, and in some grid points the disagreement between model and measurement can be quite large.The values generally scatter around the equal accumulation line, though, with mean accumulation slightly larger in the model than for the measurements.While the uncertainty range for the modelled accumulation rate is small for low ac- cumulation rates, it becomes substantially larger for higher rates.In these regions, the spread of measured values also tends to be larger.
The modelled catotelm accumulation therefore seems to capture the general pattern of measured peat accumulation, with measured and modelled accumulation rates generally scattered around the equal accumulation line, though it underestimates the large spread of accumulation measurements due to local conditions.

Carbon uptake
The carbon uptake by peatlands determines the relevance of peatlands for the global carbon cycle.With respect to the overall area north of 40°N, Fig. 6 shows the change in grid cell peat carbon density with respect to the grid cell, i.e., taking into account both the peatland area and the peat accumulation.
Modelled increases in carbon density over the course of the experiment, i.e., between 8 ka BP and present day, are quite variable (Fig. 6).Carbon density increases are largest in Scandinavia, where values reach 50 kgC m −2 .Eastern Canada and north-western Russia have lower values at 35 to 45 kgC m −2 .Most other areas, on the other hand, show substantially lower accumulations.In terms of peat height change, shown in Fig. 7, peat layers have increased by about 6 m in eastern Europe, as well as in North America south of the Hudson's Bay and on the British Isles.With peat height changes between 4 and 6 m, peat growth in Scandinavia is slightly less pronounced, which is similar along the North Sea coast.Peat accumulation in the Canadian Arctic, as well as Eastern Siberia, is rather small, though not quite negligible, but peat accumulation in the Asian wetland areas at the southern boundary of the study domain, in Kazakhstan and the surrounding areas, is very small.In principle, these changes in peat height could be compared to measurements of peat depth, which are routinely obtained.We refrain from doing so for two reasons, though.On the one hand, there seems to be a measurement and reporting bias, in that only the deepest part of the peatland is sampled (Korhola et al., 2010), therefore only showing maximum peat From the point of view of Holocene carbon cycle dynamics, the final important question is how much carbon is actually taken up and stored in peatlands.Here our model shows a substantial spread of results, depending on how the peatland area is determined with respect to the inundated area in the different sensitivity experiments.As shown in Fig. 8, the model accumulates 327.2 PgC over the 8 ka considered in our experiments if we assume that peatland area is equal to the mean inundated area.Using the minimum area, a total of 241.8 PgC is accumulated, while accumulation is 486.5 PgC for the maximum inundated area.Carbon accumulation in peatlands therefore is roughly 330 PgC for the peatland areas north of 40°N, with an uncertainty range from 240 PgC to 490 PgC.
If we take into account that some areas were still covered by an ice sheet at 8 ka BP or depressed below sea level by the previously existing ice sheet, as discussed in Section 2.6, this number is reduced to 317 PgC for the mean case, a difference of 10 PgC in comparison to the case where we neglect these changes.

Discussion
A dynamic model of wetland extent and peat accumulation, such as the one we have constructed, is extremely difficult, if not impossible to validate completely.Data on wetland extent and location is available from some national agencies, but these national inventories quite often are compiled using different measures and methodology and therefore are not comparable.Existing syntheses like the Lehner and Döll (2004) data have attempted to bridge these different method- ologies in order to determine global estimates of wetland area, but uncertainties remain large.Temporal changes in inundated area can be assessed using remote sensing data (Prigent et al., 2007), which compares favourably with Lehner and Döll (2004) for the maximum extent (Papa et al., 2010), but remote sensing data is uncertain as well.Any kind of ground cover, be it trees or snow, makes remote sensing of wetland extent impossible.A further complication is that these data sets describe the present day situation, and many wetlands in Europe as well as the more densely populated parts of North America have been drained in order to convert them to agricultural use.Our model, on the other hand, aims at determining the natural extent of wetlands, and anthropogenic disturbances are not taken into account.The approach we have chosen to determine wetland area is relatively simple and very much dependent on the hydrology of the underlying model.The latter has been evaluated positively with regard to streamflow (Gerten et al., 2004), but the rather shallow two-layer soil column in LPJ in combination with the limitation of soil moisture to field capacity and the lack of permafrost dynamics give reason for doubt with regard to soil water table dynamics.This is impossible to evaluate, though, since measurements of soil moisture exist only in very few places.Finally, the TOPMODEL approach itself is limiting as well, since a redistribution of the water within the grid cell, in combination with a grid cell mean water table that is always below the surface, implicitly limits the maximum wetland extent per grid cell.Peatland area very likely is somewhere in the range of maximum and minimum summer inundated area, more refined estimates require higher resolution measurements of peatland extent for calibration.In order to assess this uncertainty, we perform sensitivity experiments using the summer mean, maximum and minimum extent of inundated area to derive an uncertainty range for peatland area.
We are not aware of previously published attempts at determining the extent of peatlands in a dynamic way.Previous studies either focused on peat accumulation at single sites (Clymo, 1984;Frolking et al., 2001;Borren and Bleuten, 2006;Frolking and Roulet, 2007;Frolking et al., 2010;St-Hilaire et al., 2010), which does not require an estimate of wetland area, or on methane emissions (Kaplan, 2002;Gedney et al., 2004;van Huissteden et al., 2006;Bohn et al., 2007;Wania et al., 2009aWania et al., ,b, 2010;;Ringeval et al., 2010) from wetlands, which requires an estimate of inundated area, but not of the permanent wetland area.These approaches therefore either are map-based, or rely on a TOPMODEL approach similar to ours to determine the inundated area but do not specify permanently inundated areas, which are relevant for peat formation.Due to all these factors, the extent of peatlands our model determines appears to be the largest uncertainty.We take this uncertainty into account by performing three sensitivity experiments, of which we assume the experiment using summer mean inundated area as peatland extent to best capture the preindustrial peatland extent and the experiments using summer maximum and minimum areas spanning the uncertainty range.The modeled distribution of inundated area compares rather well to the Prigent et al. (2007) estimate, if just the latitudinal distribution of inundated area is considered.The full spatial distribution, on the other hand, is not captured as well, wetland area is underestimated in North America and overestimated in Eurasia.Since climate varies with latitude as a first approximation, we assume that the good fit to the latitudinal distribution also leads to a reasonable representation of peat accumulation, though the exact distribution and area of peatlands is not captured by the model.
Despite the simplicity of the approach we have chosen to model peat dynamics, which neglects potentially important factors like soil pH, exact species composition of aboveground vegetation, litter composition, etc., comparison to the little measurement data that exists appears quite favourable.Acrotelm measurements (Malmer and Wallen, 1993) vary widely, but our modelled acrotelm heights seem to fall within the ranges usually reported (Charman, 2002).Similarly, catotelm accumulation rates are in a similar range as peat core measurements (Kremenetski et al., 2003;Beilman et al., 2009;Yu et al., 2010), though our model strongly underestimates the variability in local accumulation rates.
The model estimate of total carbon accumulation in peatlands (242-486 PgC accumulated since 8 ka BP, with a best estimate of 327 PgC) is within the range of global estimates (Turunen et al., 2002;Yu et al., 2010).These latter estimates are based on point measurements which are interpreted as area averages and are upscaled by an estimate of peatland area, though, and therefore are somewhat uncertain as well.
The modelled carbon accumulation very likely is a slight overestimate during the earlier part of our experiments, since a number of factors are neglected.The peat accumulation shown in Fig. 8 is the peat accumulation over the last 8000 years, a timescale we chose since most of the continental ice sheets were melted at the time the model was initialised for.A significant number of peatlands started growing earlier than this, though, especially in areas that were not covered by ice sheets during the last glacial.In some other places like coastal areas of Scandinavia and the Hudson's Bay lowlands in Canada, peatland initiation took place later, though, usually because those areas were depressed below sea level by the overlying glacial ice sheet and postglacial rebound took some time to elevate these areas above sea level.We have estimated the effect of these latter factors by postprocessing the model output with a glacier and topography mask obtained from the ICE-5g model (Peltier, 2004), which reduces total carbon accumulation by 10 PgC for the mean extent case.The total increase of peatland areas, though, the model does not capture since we do not consider the effects of peat accumulation on topography, mainly because topography below peat deposits, which would be required, is unknown.Since published estimates of peatland area change approach 50 % for the last 8000 years, our model likely overestimates total accumulation.

Conclusions
We have extended the coupled climate carbon cycle model CLIMBER2-LPJ by a module determining permanent wetland extent and peat accumulation north of 40°N.Wetland area reflects the latitudinal distribution of wetlands well, while the longitudinal distribution is biased towards Eurasia, where area tends to be overestimated, while area is underestimated in North America.Acrotelm magnitude and catotelm peat accumulation show a smaller spread than measured values, reflecting the fact that the model neglects some of the heterogeneity in local conditions, but the overall agreement with the little measurement data that exists is good enough to give us some confidence in model results.We initialised the model for conditions at 8000 years before present, and determined the evolution of climate, wetland extent and peat accumulation until the present, assuming preindustrial conditions at present day.
Over the course of the 8000 years, wetland areas increased in Siberia and eastern Canada.Wetland areas at the southern end of the study domain either shrunk, as did wetlands in the Canadian Arctic, or did not change.The change in wetland extent reflects the changes in climate over the the course of the Holocene, a decrease in summer temperature and an increase in precipitation.It is substantially lower, though, than estimates in the literature based on peatland basal dates (MacDonald et al., 2006;Yu et al., 2010) due to the fact that peatland initiation and changes in topography due to peat growth are not considered.
During this time our model accumulates 327 PgC as catotelm peat in the areas above 40°N, with an uncertainty range of 242-486 PgC stemming from the uncertainty of peatland extent.While we have some confidence in these results, there certainly is scope for further improvement.Globally consistent data sets of peatland extent, inundation and topography would improve model calibration and validation.For the model implementation, more sophisticated hydrology models that explicitly consider permafrost processes and ground water dynamics are a requirement for an improved representation of wetland extent.

Fig. 4 .
Fig. 4. Summer (JJAS) mean wetland extent at 8 ka BP.Upper panel (a): extent as grid cell fraction.Middle panel (b): Difference to extent at 0 ka.Lower panel (c): As (a), but areas covered by ice sheet or below sea level masked (white).

Fig. 5 .
Fig. 5. Model catotelm peat accumulation rates compared to measured values.For model values, error bars are shown that span the uncertainty range in catotelm formation in the sensitivity experiments.Measurements are compiled from Yu et al. (2010), reporting accumulation rates, and Gorham et al. (2003); Kremenetski et al. (2003) and Beilman et al. (2009), where we converted basal date and peat height into accumulation rate.Averages are shown for grid cells containing multiple measurements.Mean values are shown in red.

Fig. 6 .
Fig. 6.Peat carbon density increase, relative to grid cell (not wetland fraction), accumulated over the last 8 ka.

Fig. 8 .
Fig. 8.Total carbon accumulation as catotelm peat in wetlands for cases mean, max and min.'ice' designates the accumulation, if ice sheets as well as postglacial rebound are considered.

Table 1 .
Sivapalan et al. (1987)t module.i ) the local slope at that point.The CTI can be derived from digital elevation models and near global datasets are readily available, for example the HYDRO1k data set in a resolution of 1 km(USGS, 1996).FollowingSivapalan et al. (1987), we are not using the CTI values themselves, but we rather approximate the distribution of CTI values with in a grid cell by fitting a gamma distribution to them.The parameters of this distribution can be derived from grid cell CTI mean, standard deviation and skewness.While this approach may be less precise in a few grid cells, it greatly reduces the required input data.The central equation of TOPMODEL determines the local water table z i in point i in relation to the grid cell mean water table z:

Table 2 .
Topmodel parameters f and χ max for LPJ soil classes.