Impact of an abrupt cooling event on interglacial methane emissions in northern peatlands

Rapid changes in atmospheric methane (CH 4), temperature and precipitation are documented by Greenland ice core data both for glacial times (the so called DansgaardOeschger (D-O) events) as well as for a cooling event in the early Holocene (the 8.2 kyr event). The onsets of D-O warm events are paralleled by abrupt increases in CH 4 by up to 250 ppb in a few decades. Vice versa, the 8.2 kyr event is accompanied by an intermittent decrease in CH 4 of about 80 ppb over 150 yr. The abrupt CH 4 changes are thought to mainly originate from source emission variations in tropical and boreal wet ecosystems, but complex process oriented bottom-up model estimates of the changes in these ecosystems during rapid climate changes are still missing. Here we present simulations of CH 4 emissions from northern peatlands with the LPJ-Bern dynamic global vegetation model. The model represents CH 4 production and oxidation in soils and transport by ebullition, through plant aerenchyma, and by diffusion. Parameters are tuned to represent site emission data as well as inversion-based estimates of northern wetland emissions. The model is forced with climate input data from freshwater hosing experiments using the NCAR CSM1.4 climate model to simulate an abrupt cooling event. A concentration reduction of ∼ 10 ppb is simulated per degree K change of mean northern hemispheric surface temperature in peatlands. Peatland emissions are equally sensitive to both changes in temperature and in precipitation. If simulated changes are taken as an analogy to the 8.2 kyr event, boreal peatland emissions alone could only explain 23 % of the 80 ppb decline in atmospheric methane concentration. This points to a significant contribution to source changes from low latitude and tropical wetlands to this event.


Introduction
CH 4 is a greenhouse gas contributing to the ongoing global warming. CH 4 concentrations have increased from their preindustrial value of ∼ 700 ppb to approximately 1800 ppb in 2008  due to anthropogenic CH 4 emissions. Air enclosures in polar ice cores allow for the reconstruction of methane variations in the past and show that CH 4 changed in concert with Northern Hemisphere temperature during the glacial/interglacial transitions  as well as during rapid climate variations. Examples of abrupt CH 4 variations recorded in Greenland ice cores are the 8.2 kyr event in the early Holocene ( Fig. 1) (Blunier et al., 1995;Chappellaz et al., 1997;Spahni et al., 2003;Kobashi et al., 2007), the Dansgaard-Oeschger events (D-Oevents) in the glacial period (Chappellaz et al., 1993;Brook et al., 2000;Blunier and Brook, 2001;Flückiger et al., 2004;Huber et al., 2006;Baumgartner et al., 2012) or the Younger Dryas (YD) cooling event (Chappellaz et al., 1993;Baumgartner et al., 2012). During the last glacial period Greenland temperature and CH 4 are found to covary within ∼ 30 yr, which is the limit of ice core data resolution and the width of the age distribution of air enclosed in bubbles in that ice (Huber et al., 2006). This suggests that CH 4 sources responded synchronously to the rapid climate changes. The nominal sensitivity of methane concentration variations on changes in Greenland temperature derived from ice cores is about ∼ 22 ppb K −1 for the 8.2 kyr event (Kobashi et al., 2007), ∼ 7-16 ppb K −1 for onsets of D-O events (Huber et al., 2006) and Evolution of Greenland temperature and atmospheric methane concentration during the 8.2 kyr event. The measured δ 15 N signature allows for a quantitative temperature reconstruction and represents a temperature change of about 3 • C over Central Greenland which is accompanied by an 80 ppb decline in atmospheric CH 4 (Kobashi et al., 2007). recorded in Greenland are not necessarily representative for the temperature changes in boreal and even less in tropical wetland sources nor for precipitation changes.
While the origin of the abrupt climate variations can be explained by a shut-down or resumption of the North Atlantic thermohaline circulation (Vellinga and Wood, 2002;McManus et al., 2004;Ellison et al., 2006), the origin of atmospheric methane variations is much less clear and abrupt climate variations are still not well understood (Clement et al., 2008). Most ice core studies attribute a change in methane concentration to a change in methane emissions from both boreal and sub-tropical to tropical wet ecosystems (Chappellaz et al., 1997;Fischer et al., 2008). The interpolar concentration gradient (IPG) was reconstructed from ice core methane data from Greenland and Antarctica (Brook et al., 2000;Baumgartner et al., 2012). This gradient together with methane isotopic information indicates important sources during both warm and cold periods located in high northern latitudes with a tendency to a higher IPG during interstadials (Bock et al., 2010;Brook et al., 2000;Dallenbach et al., 2000). For the 8.2 kyr event, the interpolar gradient suggests that the CH 4 emissions must have been reduced in the early phase of the 8.2 kyr event in the Northern Hemisphere, followed by sub-tropical and tropical emission reductions (Spahni et al., 2003). However, a quantitative estimate of the response of boreal CH 4 sources on rapid climate changes based on the IPG is subject to large uncertainties due to the transient character of the changes and the simplified models used in this approach.
The dominant and most directly responding type of methane emitting ecosystems in the northern high latitudes, i.e., the latitudes with the strongest temperature response to rapid changes in the Atlantic Meridional Overturning Circulation (AMOC), are boreal peatlands (Dallenbach et al., 2000;Baumgartner et al., 2012). Northern high-latitude peat-land and permafrost soils are associated with large carbon stocks that serve as a substrate for microbial methane production (Gorham, 1991;Roulet, 2000). For the understanding of the covariations of atmospheric methane and climate during abrupt events it is, thus, crucial to investigate the carbon cycle of northern high-latitude peatlands (northern peatlands hereafter). The carbon cycle response of northern peatlands regulates the peatland methane emissions and thus directly affects the global methane concentration within its atmospheric mixing and lifetime of ∼ 2 and ∼ 10 yr, respectively (Lelieveld et al., 1998;Denman et al., 2007;Prather, 2012;Levine et al., 2011). Accordingly, the goal of this study is to investigate the CH 4 sensitivity of the northern peatlands to rapid climate cooling.
There are two different approaches to model methane emissions from wet ecosystems for present day: bottom-up models that represent the processes leading to methane emissions in a mechanistic way (e.g. Cao et al., 1996;Walter et al., 2001a;Wania et al., 2010;Ringeval et al., 2010) and top-down or inverse models which use observations of atmospheric methane concentrations and transport models to quantify methane emissions (e.g. Houweling et al., 1999;Pison et al., 2009;Bergamaschi et al., 2009). Combining both methods can help to narrow down the uncertainties of global wet ecosystem methane emissions (Spahni et al., 2011).
For past time periods, before the start of direct atmospheric methane monitoring in the 1980s, the spatially available CH 4 data is restricted to firn and ice core data from Greenland and Antarctica. Global inversions from atmospheric concentrations to emissions using atmospheric chemistry transport models are, thus, not applicable for past time periods. Ice-core based top-down approaches using simple box models (e.g. Fischer et al., 2008) can only explain the inter-hemispheric gradient, but not distributions within a hemisphere. Accordingly, biogeochemical process modelling of methane emissions and comparisons of resulting atmospheric concentrations to ice core data are the only way to estimate natural methane emission distributions.
Previous paleo-modelling studies have shown that simple methane models are capable of simulating emissions on longterm glacial-interglacial time scales Kaplan et al., 2006;Weber et al., 2010;Singarayer et al., 2011) and are also able to capture methane emissions during abrupt events (van Huissteden, 2004;Hopcroft et al., 2011), but find it difficult to replicate the magnitude of change in atmospheric CH 4 for abrupt events. Hopcroft et al. (2011) can not explain the full magnitude or rapidity of methane emission changes during D-O-events and believe their model to be too insensitive for abrupt events or to be missing important mechanism. van Huissteden (2004) finds doubled methane fluxes from northern European wetlands during interstadials, but attributes big uncertainties to the peatland area. The long-term simulations confirm the finding that atmospheric methane concentrations parallel temperature reconstructions from Antarctic ice cores also on glacial-interglacial time scales of the past 800 000 yr (Spahni et al., 2005;Loulergue et al., 2008). A drawback of previous paleo-modelling studies for methane is the limited representation of processes and interactions related to the interplay of vegetation and soil dynamics, soil temperature, soil hydrology, permafrost thawing, methane production, oxidation and transport from the soils to the atmosphere by a variety of mechanisms. For example, Kaplan et al. (2006) calculate the methane emissions as a fixed ratio of the heterotrophic respiration. Valdes et al. (2005); Kaplan et al. (2006);Singarayer et al. (2011) consider changes in the atmospheric CH 4 lifetime in order to explain glacial-interglacial changes in CH 4 . According to the latest chemistry model studies these changes in CH 4 lifetime appear to be very limited (Levine et al., 2011).
Here, we further develop and apply a process-oriented CH 4 module (Wania et al., 2009a(Wania et al., ,b, 2010 within the Bern version of the LPJ dynamic global vegetation model. Soil temperature at different depths, including freeze-thaw cycles and water table position is explicitly calculated and CH 4 transport by soil diffusion, ebullition and through plant aerenchyma is simulated mechanistically. The LPJ land model is forced with the climate output from idealised freshwater hosing experiments with the NCAR CSM1.4 climate model. Freshwater is added in the model to the preindustrial northern North Atlantic surface ocean. In response, the Atlantic Overturning Circulation collapses and surface air temperature in the Northern Hemisphere is reduced. We will discuss the results of our idealised simulations in relation to the 8.2 kyr event, but we emphasise that a direct comparison between model results and proxy data for the 8.2 kyr event is hampered as environmental conditions at 8.2 kyr ago and before the onset of the industrialisation were different. For the early Holocene, the initiation of boreal peatlands is reasonably well known from recent observations (Yu et al., 2010). Since peatland initiation peaks at 10 kyr BP we apply present day boreal peatland distribution (Tarnocai et al., 2009). A comparison of this distribution with the ice coverage 8200 yr ago (Peltier, 2004) shows that only very small areas of present day peatlands were covered by land ice at that time. Additionally, the warm Holocene climate conditions prevailed already for more than 1000 yr allowing for an initiation of peat expansion to interglacial conditions. The outline of this paper is as follows. In Sects. 2.1 and 2.2, we describe the LPJ-Bern global dynamic vegetation model Joos et al., 2004;Gerber et al., 2003;Strassmann et al., 2008;Stocker et al., 2011), the implementation and improvements to LPJ-WHyMe (Wania et al., 2009a,b) and a recently developed methane module (Wania et al., 2010). Section 2.3 presents the calibration of the model with modern site data and describes the experimental setup for control simulations and for simulations forced with climate input from an ensemble of freshwater hosing model experiments (Bozbiyik et al., 2011). Results are presented in Sect. 3 and include simulated changes in CH 4 emissions for the set of freshwater hosing experiments, the attribution of change in CH 4 emissions to temperature and precipitation variations, and sensitivity to model parameters.

LPJ-Bern
The LPJ-Bern (hereafter LPJ) model is a subsequent development of the Lund-Potsdam-Jena dynamic global vegetation model (LPJ-DGVM; that combines process-based, large-scale representations of terrestrial vegetation dynamics and land-atmosphere carbon and water exchanges in a modular framework. Vegetation is defined by plant functional types (PFTs) each with its own set of parameters describing growth, carbon uptake and mortality. LPJ has been applied previously in paleo studies (e.g. Joos et al., 2004;Gerber et al., 2004) and in simulations assessing the anthropogenic climate perturbation and the impact of human induced landuse (e.g. Joos et al., 2001;Strassmann et al., 2008;Stocker et al., 2011).
For this study we implemented the boreal peatland and methane module based on LPJ-WHyMe (Wetland Hydrology and Methane, Wania et al., 2009aWania et al., ,b, 2010) that derives its input data from the LPJ. It includes new features for grid cells with a (prescribed) partition of peatland: 8 soil layers, permafrost dynamics with freezing and thawing (including a soil temperature solver to simulate temperature as a function of depth), peatland hydrology (active layer depth and water table position), peatland vegetation as two additional PFTs (flood-tolerant C 3 graminoids and Sphagnum moss), a slow-down of decomposition under inundation and the addition of root exudates (Wania et al., 2009a,b).
The main differences between our methane routine and the original LPJ-WHyMe version 1.3.1 (Wania et al., 2010) is a more mechanistic ebullition mechanism that includes also the partial pressure of CO 2 for triggering an ebullition event (see Sect. 2.2). The carbon balance over all layers is now preserved after every gas diffusion time step, whereas in LPJ-WHyMe the carbon balance was closed at the end of the year with a correction factor.

Methane routine in peatlands of LPJ-Bern
The carbon module simulates carbon allocation in vegetation, above and below ground litter and fast and slow soil carbon pools, as well as soil organic matter dynamics . Heterotrophic respiration (HR) is calculated based on the size of the litter and soil pools, regulated by soil moisture and soil temperature at a depth of 25 cm. HR is then distributed to 8 peat layers according to the root distribution, where it is used to compute methane production. In each layer methane production is calculated as the product of HR, the ratio of CH 4 to CO 2 production under anaerobic 1966 S. Zürcher et al.: Impact of an abrupt cooling event on interglacial methane emissions conditions (f CH 4 /CO 2 ), and the anoxity of the soil for each soil layer (Wania et al., 2010).
Depending on the water table position, the soil temperature and the presence of O 2 , the produced methane can be oxidised to CO 2 , accumulated or transported to other layers or released to the atmosphere. The three implemented transport processes are plant-mediated transport, diffusion and ebullition. Plant-mediated transport is implemented as described in (Wania et al., 2010): the flood-tolerant C 3 graminoids (one of two PFTs growing in peatlands) are adapted to inundation by the presence of aerenchyma, gas-filled tissue, through which gas can be transported by diffusion. Plant-mediated transport depends on temperature, water content and the actual gradients of CH 4 , CO 2 or O 2 . The tiller radius r tiller defines the open area available for plant-mediated transport and is a tunable parameter.
Diffusion transports O 2 , CO 2 and CH 4 between the 8 soil layers. The temperature-dependent diffusion coefficients of all gases are set to their molecular diffusivities in air or water (Lerman, 1979;Broecker and Peng, 1974), depending on the water table position in the acrotelm (top 30 cm of soil). Peatland soil layers below that horizon are assumed to be completely water saturated (catotelm). Water can freeze if the soil layer temperature falls below 0 • C and all transport mechanisms are prohibited. CH 4 is not incorporated into the ice during freezing leading to an enrichment of CH 4 in the remaining water and an eventual release by ebullition.
The new ebullition mechanism calculates how much of the CH 4 , CO 2 and N 2 is in the gas phase and how much is dissolved. Following Henry's law (FechnerLevy and Hemond, 1996;Rosenberry et al., 2003;Tokida et al., 2005Tokida et al., , 2007 the amount of dissolved gases in each layer depends on partial pressure of CH 4 , CO 2 and N 2 , hydrostatic pressure and soil temperature. The available volume (V available ) within each layer i is given by the layer thickness and the porosity for fully water saturated conditions. The total amount of N 2 is assumed to be constant, since it is not simulated by our model. However, the partitioning of N 2 between water and soil air changes. Assuming the gaseous volume of N 2 to be 1 % (lower bound estimate, e.g. Shannon et al., 1996) of the available volume in a layer (liquid phase and gas phase), we calculate with Henry's law the amount of N 2 dissolved in each layer (Tokida et al., 2007). We then add the produced or transported CH 4 and CO 2 to the layer and recalculate the partial pressure for all gases (Eq. 1, Tokida et al., 2007;Kellner et al., 2006): where P x,i is the partial pressure of gas species x in layer i, n x,i is the total number of moles of gas x, V g,i is the gas volume, V w,i is the volume of liquid water, H x is the dimensionless Henry's law constant for gas x (Yamamoto et al., 1976;Wiesenburg and Guinasso, 1979;Weiss, 1970), R is the universal gas constant, and T i is the temperature in layer i (see Appendix A). We assume that during ebullition the gas volume, V g,i , is equal to 15 % and the liquid water volume, V w,i , 85 % of the available volume (FechnerLevy and Hemond, 1996;Rosenberry et al., 2003;Tokida et al., 2005). In the code, ebullition is triggered if the sum of all partial pressures (pN 2 + pCO 2 + pCH 4 ) (Eq. 1 with V g,i = 15 % · V available ) is larger than the sum of atmospheric and hydrostatic pressure. The number of moles of CO 2 and CH 4 released per ebullition event is assumed to be 1 % of the available volume times the CO 2 and CH 4 partial pressure divided by R T i , respectively. In other words, a fifteenth of the CO 2 and CH 4 amount in the gas phase is released per ebullition event. This assumption and that the gas volume is equal to 15 % of the available volume are not critical for the total emissions. Instead these parameters modulate the frequency and the magnitude of ebullition events on a short scale (see parameter tuning, Sect. 3.4).
Oxygen is transported into soils by diffusion through open soil pores and through the aerenchyma of the flood tolerant graminoids. Only a fraction (f oxid ) of the available oxygen is used to oxidise CH 4 to CO 2 . The remaining O 2 is assumed to be consumed by other electron acceptors. Again f oxid is a tunable model parameter.

Calibration and site validation
The implementation of LPJ-WHyMe into the LPJ-Bern, the code improvements and the adaption of the ebullition scheme calls for a re-evaluation of flux densities at the site level and a retuning of model parameters. We compare simulated CH 4 emissions to flux measurement at 7 sites. A detailed description of the sites and their environmental conditions can be found in Wania et al. (2010). Measured seasonal methane fluxes over a specific year are available for each site (see Fig. 2) (Bubier et al., 1998;Ding et al., 2004;Saarnio et al., 1997;Alm et al., 1999;Granberg et al., 2001;Johansson et al., 2006;Jackowicz-Korczynski et al., 2010;Shannon and White, 1994;Dise, 1993). We force LPJ with local temperature and precipitation for the sites in Abisko and Sanjiang and we use CRU-data (Climate Research Unit's (CRU) TS 2.1, climatology from 1901 to 2002, (New et al., 1999;Mitchell and Jones, 2005)) of the corresponding grid cell and year for the other sites for the calibration. LPJ is forced with monthly climate data, i.e., temperature (Fig. 4b), precipitation ( Fig. 4c) and cloud cover, which are interpolated to daily values. The climate records include seasonal and interannual variability and are repeated over the spin-up period to reach a stable steady state.
The three most important factors for total methane emissions and the attribution to each emission pathway in the model turned out to be (i) the CH 4 /CO 2 -factor f CH 4 /CO 2 , followed by (ii) the oxidation fraction f oxid and (iii) the tiller radius r tiller of peatland grasses. Other factors controlling methane emissions are the exudate turnover rate, the moisture response, which is also used to calculate decomposition Daily methane emissions measured at seven sites compared to LPJ model output analogue to a previous model calibration (Wania et al., 2010). All our model runs use the same optimised parameter set for all sites. The respective root-mean-square errors (RMSE) calculated from daily differences between spline fits through the measurements and model output over the period where the measurements exist are given in the top left corner of each box. Note that at the Degeroe site ebullition was not measured and was, therefore, also excluded in the model output. The Abisko record is the only one with daily measurements and only the spline fit is plotted for the purpose of better clarity. Model results are plotted as smoothing spline function with a cut-off period of two months (blue line), observations are plotted as daily values (black dots) and a smoothing spline function with a cut-off period of two months (red line). Note the different scales on the y-axes. rates, the leaf-to-root ratio, which influences the crosssectional area of aerenchyma available for plant-mediated gas transport, the tiller porosity which influences the area for plant-mediated transport and parameters regulating the frequency when ebullition occurs and the volume released in a single ebullition event. Except for the two new parameters of our ebullition routine, we use the same values as (Wania et al., 2010) for these six parameters that have proven to be less important.
The CH 4 /CO 2 -factor, linking heterotrophic respiration and CH 4 production, directly affects the total production of CH 4 and, thus, total annual emissions. The oxidation fraction is the percentage of available O 2 that is used to oxidise CH 4 in each layer. Its effect is similar to the CH 4 /CO 2 -factor, but weaker. The tiller radius of peatland grasses influences directly the plant-mediated transport of O 2 , CO 2 and CH 4 . These three parameters were varied over the range of values found in the literature and chosen to minimise the root-meansquare error between model results and data for the seven sites. Data and model results were smoothed using a spline approximation (Enting, 1987) with a cut-off period of two months. The spline covers the time period where the respective measurements exist. Nominal daily data were applied to compute that root-mean-square error for each site between splined measurements and splined model output. The root-mean-square error for all sites is computed from the time period weighted (spanned by measurements) site errors. Local methane emissions for the best parameter set in this sensitivity study are shown in Fig. 2, together with the site measurements and the RMSE as computed from the spline fits. These optimised parameter set is then used to simulate global CH 4 emissions in the model experiment. The same set of parameters is used across all sites although a better fit with observational data could be achieved when parameters would be adjusted for each site individually. We find a slightly improved agreement between model and observations compared to results presented by (Wania et al., 2010) (see Appendix B). This is likely related to the elimination of some coding errors and the improved representation of ebullition.
The optimised parameters are 0.17 g C (g C) −1 for the ratio of CH 4 /CO 2 production under anaerobic conditions, 0.5 for the oxidation fraction, i.e., the fraction of available oxygen used for methane oxidation, and 0.0035 m for the tiller radius. Further, we make the assumptions that ebullition occurs when 15 % of the available pore volume is gas and that 1 % of the available volume is released as gas in a single ebullition event. The influence of these parameters on modelled changes in methane emissions after a climate perturbation will be presented in the result section.
It is possible to get similar root-mean-square deviations between model results and observations with different combinations for the production ratio and the oxidation fraction. Very different fractions of oxygen consumption by methanogenesis have been reported, reaching from 0.2 to 1 (Wania, 2007;Frenzel and Rudolph, 1998;Strom et al., 2005). Therefore, we set our oxidation fraction according to literature values and in agreement with Wania et al. (2010) to the value 0.5. Changing parameters, governing plant-mediated transport or ebullition, changes the attribution to the different pathways, but does only weakly influence the total emission per unit area. To quantify the distribution better, either more site measurements would be needed or some further constraints like the 13 C signature of the total emissions.
Despite a rather well representation of the amplitude in CH 4 emissions in the course of the year a significant phase lag between observed and modelled CH 4 emissions is found for some sites as illustrated in Fig. 2. This may be partly attributed to the coarse resolution of the CRU data that drives the model for Boreas, Salmisuo, Degeroe, Michigan and Minnesota. In these cases CRU data are not necessarily representing the local conditions. For example, a delayed rise in spring temperature leads to a time lag in soil thawing and, thus, methane production and emissions. Effects like water runoff from other regions are not taken into account either. However, even if we use site climate data at Sanjiang, emissions show a phase lag. This indicates that the seasonal change in the different emission processes is not accurately reflected by our model for all sites. The mismatch may be related to a possible neglect of heterotrophic respiration at shallow soils that thaw early in the season. The methane routine distributes substrate for methanogenesis derived from the total heterotrophic respiration according to the root distribution. In case a layer is still frozen, the substrate is put into the next unfrozen layer above. Analysis of the model output for Sanjiang shows that the soil thawing in deep layers has not began when the measured emissions should start and the substrate for production is transferred to upper layers where it is almost completely oxidised. The HR is directly correlated to temperature in LPJ and the simulated temperature profile seems to have a lag respective to the site profile.
Our sensitivity simulations presented in Sect. 3.4, also show that the partitioning between the different emission pathways (diffusion, plant transport, ebullition) has an influence on the seasonality of the emissions and is sensitive on the choice of model parameters, however, the total annual emission is not. As illustrated in Fig. 3, the time-integrated CH 4 flux is well represented by the model across the different sites and a tight correlation (R 2 = 0.92) is found between simulated and measured cumulative site emissions. This indicates that the model has a good skill in representing annual emissions across a wide range of environmental conditions in the boreal zone. In conclusion, total annual emissions are simulated well across the different sites, while deviations be- tween simulated and observed seasonal cycles remain (see also Wania, 2007).

Input data, spin-up, and setup for transient experiments
Spatially-resolved input data include the global land mask, the soil type (Zobler, 1986) (always "organic" on peatland) and a map of present day peatlands, here on a resolution of 3.75 • longitude and 2.5 • latitude. We prescribe the distribution and fraction of peatland area per grid cell (Tarnocai et al., 2009). The Tarnocai NCSCD dataset (Tarnocai et al., 2007) includes histels and histosols that were rasterised and regridded to the model resolution (Fig. 4a), and scaled globally to match an area of 1.048 million km 2 in the permafrost affected region in North America (Tarnocai et al., 2009). This scaling results in a Northern Hemisphere peatland area of 2.81 million km 2 . This area is larger than the 2.06 million km 2 derived from the IGBP-DIS map as used in (Wania et al., 2010). The number of wet days per month is prescribed from the CRU climatology (New et al., 1999;Mitchell and Jones, 2005). Atmospheric CO 2 is set to a constant value of 279 ppm as found in ice cores for preindustrial conditions (Monnin et al., 2004). LPJ is forced with monthly data of temperature (Fig. 4b), precipitation (Fig. 4c) and cloud cover. Here, climate output from a control run and a 3-member ensemble of freshwater hosing simulations recently conducted with the Climate System Model, version 1.4, (CSM1.4) of the National Centre for Atmospheric Research (NCAR) (Bozbiyik et al., 2011) is applied to simulate methane emissions and emission changes in northern peatlands during a cooling event  (Tarnocai et al., 2009) for present conditions (a); Northern Hemisphere climate patterns averaged over 1961-1990CRU data (New et al., 1999Mitchell and Jones, 2005) for temperature (b) and precipitation (c), at the peatland grid cells. similar to the 8.2 kyr event with respect to temperature change over Greenland and Northern Hemisphere land area (Fig. 1). All CSM1.4 freshwater hosing simulations branch from the control simulations that features a stable climate and a constant freshwater input flux anomaly of 1 Sv (1 Sverdrup = 106 m 3 s −1 ) is applied over 100 yr in the Northern North Atlantic. Freshwater input starts 126, 263 and 305 yr after the start of the control in ensemble simulation s1, s2 and s3 and simulations are continued for several centuries. We run several ensembles with the same forcing to distinguish between forced changes and changes that are related to internal climate variability such as the Arctic Oscillation.
The monthly temperature and precipitation data of each CSM1.4 simulation are splined for each grid cell and each calendar month separately, using a 20-yr cut-off frequency (Enting, 1987) (Fig. 5). Anomalies for each month relative to the corresponding monthly mean of the CSM1.4 control run are computed. Monthly anomalies are then added to the CRU climatology (New et al., 1999;Mitchell and Jones, 2005). As the LPJ model is sensitive to the absolute input climatology, this procedure eliminates biases of the NCAR CSM 1.4 (Kiehl et al., 2006) for modern climate relative to the CRU dataset and, thus, corrects for climate-related biases in grid cell vegetation distribution and carbon stocks.
All LPJ simulations are driven with a combination of the repeated 31-yr cycle of CRU data (TS 2.1, detrended climatology from 1960to 1990, New et al., 1999Mitchell and Jones, 2005) as a baseline and monthly anomalies from the NCAR CSM 1.4 model under preindustrial conditions. The spin-up is done for all runs with the CRU climatology and the anomalies from the control run and is 1000 yr long. To accelerate the spin-up procedure, equilibrium soil carbon stocks are computed after 400 yr based on average input fluxes and decomposition rates. Then spin-up is continued for another 600 yr.
Simulated annual boreal methane emissions are scaled to 30 Tg CH 4 yr −1 for the end of the spin-up period, which is a reasonable value for the present day and the preindustrial (Spahni et al., 2011, and therein). The scaling factor is 0.6. Since the area of effective methane emissions is highly heterogeneous, the simulated emissions have to be scaled accordingly. The prescribed peatland map only indicates the location of peat and the total area, but not its small scale structure. Peat layers on hummocks are emitting much less than peat layers in lawns. Accordingly, a perfect match of the total modelled peatland emission with best estimates from observations cannot be expected.

Preindustrial emissions from boreal peatlands and interannual variability
There is significant year-to-year climate variability and boreal methane emissions vary accordingly in the control simulation (Fig. 6). The range in annual emissions lies between 27 and 37 Tg CH 4 yr −1 and the standard deviation of annual emissions around the long-term mean is ± 1.8 Tg CH 4 yr −1 in the control. The spread in annual emission of 27 to 37 Tg CH 4 yr −1 may be compared with the spread of 29 to 37 Tg CH 4 yr −1 for Northern Hemisphere extratropical wetlands found by Chen and Prinn (2006) or 25 to 34 Tg CH 4 yr −1 by Spahni et al. (2011). Clearly, boreal wetland emissions have the potential to influence interannual variations in atmospheric methane by a few ppb.
Multi-year average methane emissions from peatlands are simulated to typically range between 0 to 50 g C m −2 yr −1 Fig. 5. Temperature and precipitation anomalies in boreal peatlands simulated by NCAR CSM 1.4 in response to freshwater hosing in the North Atlantic. Top: anomalies represent the differences in 50 yr averages of monthly data before and after the freshwater release of ensemble run 1. Bottom: anomalies represent spatial-averages using peatland area of each grid cell as weights and with respect to the mean of the control simulation (2 months running average).
( Fig. 7a) for the control. In general, emissions per unit area are relatively high in Northern Europe and in grid cells along the southern boundary of the peat zone in North America. Relatively low emissions per unit area are simulated in the Canadian Arctic and parts of Siberia, reflecting a generally short growing season and harsh climatic conditions (Fig. 4).

Temporal and spatial changes in methane emissions after a climate perturbation
We next turn to the response to the simulation of an idealised abrupt cooling event. The additional freshwater flux applied in the NCAR CSM1.4 ensemble simulations results in a collapse of the Atlantic Meridional Overturning Circulation and of the associated northward heat transport into the North Atlantic region with no sign of recovery until the end of the simulations. In return, simulated temperature and precipitation in the North Atlantic and over the adjacent continents decrease relative to the stable CSM1.4 control run (Fig. 5). On annual and spatial average, using peatland area of each grid cell as weights, temperature over peatlands (40-90 • N) decrease by about 2.5 • C within decades after adding the freshwater and remain low until the end of the simulations (Fig. 5c). Largest temperature anomalies over peatlands are simulated for the British Isles and Scandinavia, whereas simulated temperature changes remain small in northern North America and the Siberian lowland, except in near-costal regions (Fig. 5a). A slight warming is even simulated for the major peat regions in North America. On average, precipitation drops over peatland areas by 3 to 5 mm per month (Fig. 5d), with largest changes again on the British Isles (Fig. 7b). While the moss PFT was slightly dominant before the start of the freshwater input (foliar projective cover of 55 % for the moss PFT and 45 % for the grass PFT), it suffers more under the dryer and cooler conditions (∼ 43 % both afterwards). The biggest changes are in the Atlantic regions and the retreat of the sum of both PFTs is about 10 %.
In response to the widespread decrease in temperature and precipitation, peatland CH 4 emissions are reduced by about 19 % compared to the control run. In accordance with the evolution of climate, emissions do not recover in the centuries after the perturbation (Fig. 6). The spatial pattern of response reflects the change in climate with largest decrease in emissions per unit area of up to 25 g C m −2 yr −1 in Northern Europe and small changes in Eastern North America and in the Siberian lowlands (Fig. 7b). A slight increase in emissions is even found south of Hudson Bay, where an increase in temperature and small precipitation changes are simulated after the freshwater perturbation.
The results are consistent across the three ensemble members. The signal from the freshwater hosing is much larger than internal variability and, therefore, we will rely on run 1 for further discussion. The decrease in emissions from boreal peatlands by 6 Tg CH 4 yr −1 corresponds to a decrease of 18 ppb in atmospheric methane (assuming a lifetime of 8.4 yr (Stevenson et al., 2006;Prinn, 1994;Lelieveld et al., 1998;Levine et al., 2011), and a unit conversion factor of 2.78 Tg CH 4 ppb −1 ).

Influence of changes in temperature and precipitation on emission changes
Next, we investigate the sensitivity of methane emission changes to changes in individual climatic drivers. This is done by factorial simulations, where either precipitation anomalies or temperature anomalies from the NCAR CSM1.4 simulations are applied. The results show that the changes in the amount of precipitation and in temperature are about equally responsible for the drop in the total methane emissions from the northern peatlands (Fig. 8a); total emission decreases by about 4 Tg CH 4 yr −1 in each simulation. The combined effect of temperature and precipitation is not the sum of the individual contributions to the CH 4 emission response. The main reason for this is that cells where the production and transport of methane is already inhibited by low temperatures cannot further be reduced in their activity by missing water and vice versa; grid cells with no activity due to too dry conditions are inactive independent of a further decrease in temperature. In conclusion, both temperature and precipitation changes influence CH 4 emissions significantly.

Governing processes and their sensitivity to model parameters
We now turn our attention to the importance of individual processes for methane emissions from boreal peatlands. Gross primary productivity (GPP) of the two plant functional types existing on boreal peatlands decreases by 14 % from 1460 to 1260 Tg C yr −1 as evaluated by average values before and after the freshwater pulse. This sustained decrease is of the same order as interannual variability in GPP of ± 80 Tg C yr −1 (± 1 sdv). Net primary productivity (NPP) drops by 12 % from 840 to 740 Tg C (interannual: ± 40 Tg C) and total (aerobic and anaerobic) heterotrophic respiration (HR) by 13 % from 840 to 730 Tg C (interannual: ± 30 Tg C). This decrease is driven by the general decrease in temperature and precipitation in roughly equal parts as inferred from the factorial simulations described above. In turn, the simulated total soil carbon inventory in boreal peatlands of 536 Pg C (interannual: ± 0.1), including the litter above and below ground, the exudates, and a fast and a slowly decaying soil carbon pool, remains virtually unchanged (532 ± 0.1 Pg C (−1 %)). The changes in the fast carbon pools only (exudates, litter and fast soil pool) are 2.6 %. The decrease in respiration is, thus, the result of a decrease in the mean decomposition rate of organic material and not a decrease in soil carbon inventory. Note that the temperature and precipitation changes in the area of large soil carbon stocks are moderate ( Figs. 9 and 5a,b). The soil carbon inventory is in a steady state under preindustrial conditions and realistic compared to available data (MacDonald et al., 2006;Yu et al., 2010;Tarnocai et al., 2009). The average yearly water table position over all peatland gridcells drops from 4 mm below the surface to 6 mm. The biggest changes are in the Atlantic regions and the adjacent areas with up to 110 mm change in the annual mean. In our model, production of methane is linearly scaled with the rate of heterotrophic respiration and total preindustrial emissions from peatlands are scaled to 30 Tg CH 4 yr −1 as   8. (a) Total annual methane emission from northern peatlands in the standard scenario (red) compared to factorial runs where only changes in temperature (cyan) or only in precipitation (blue) were taken in account. The total annual methane production rate (violet) and soil oxidation rate (yellow) is plotted as well for the standard scenario. (b) Contribution of the three different transport pathways to the total emissions. Results are obtained with the best fitting parameters (solid lines) and a run with a doubled tiller-radius (dashed lines). The total emissions and emissions by diffusion are identical for the standard run and the doubled tiller-radius run, whereas there is a redistribution of plant-mediated emissions and emissions by ebullition. Climate forcing data are from ensemble simulation s1 and model output is smoothed with a spline. 49 Fig. 8. (a) Total annual methane emission from northern peatlands in the standard scenario (red) compared to factorial runs where only changes in temperature (cyan) or only in precipitation (blue) were taken in account. The total annual methane production rate (violet) and soil oxidation rate (yellow) is plotted as well for the standard scenario. (b) Contribution of the three different transport pathways to the total emissions. Results are obtained with the best fitting parameters (solid lines) and a run with a doubled tiller-radius (dashed lines). The total emissions and emissions by diffusion are identical for the standard run and the doubled tiller-radius run, whereas there is a redistribution of plant-mediated emissions and emissions by ebullition. Climate forcing data are from ensemble simulation s1 and model output is smoothed with a spline. described in the method section. For the standard set of parameters, preindustrial total methane production amounts to 42 ± 2 Tg CH 4 yr −1 of which 29 % (12 ± 1 Tg CH 4 yr −1 ) are removed within the soil by oxidation. After the freshwater pulse, production drops on average to 35 ± 2 Tg CH 4 yr −1 of which again about 31 % (11 ± 1 Tg CH 4 yr −1 ) are oxidised and 69 % emitted (Fig. 8a). 5.0 % of total NPP and HR is converted into methane (42 Tg C/840 Tg C) under preindustrial conditions. These fractions stay almost the same (4.7 %) under the climate prevailing after the freshwater input. A simple approach using the decline in RH as a scaling for the dropdown in methane would lead to only a decline in methane emissions by about 13 % as mentioned above, while we find a decline of about 19 % with the coupled process modelling approach (Table 1). This suggests that a fixed scaling of CH 4 emissions to total heterotrophic respiration may not be adequate to simulate variations in CH 4 emissions under varying climate conditions such as reconstructed for Dansgaard-Oeschger events or glacial-interglacial cycles.
The total methane emission flux is the sum of plant transport, diffusion and ebullition. For the standard parameter choice, ebullition and plant transport are the dominant transport pathways (solid lines in Fig. 8b) and each contribute about 35 to 40 % to the total emission. The percentage contribution of these pathways to total emissions remains roughly constant after the climatic perturbation. The fraction of methane that is oxidised remains also fairly constant.
Next, we analyse the sensitivity of the simulated emissions and their changes to the choice of model parameters in the methane module (Table 1). In a suite of sensitivity experiments, a parameter is varied individually around its standard value. Each run is rescaled with the same factor that scales the standard parameter run to 30 Tg CH 4 yr −1 at the end of the spin-up. Parameters varied are the CH 4 /CO 2 production factor, the oxidation factor, i.e., the percentage of available O 2 that is used to oxidise CH 4 in each layer, and the tiller Table 1. Average methane emissions before perturbation (EBP) and after perturbation (EAP) from northern peatland; standard parameters: f CH 4 /CO 2 = 0.17 g C g C −1 , f oxid = 0.5, tiller radius r tiller = 0.0035 m, V gas (fraction of gas from actual water volume) = 0.15, V (fraction for ebullition) = 0.01 radius of the peatland grass PFT and the two parameters controlling the ebullition mechanism. The CH 4 /CO 2 production factor and the oxidation fraction (range 0.2 to 0.6) scale the total methane emissions from boreal wetlands almost linearly under preindustrial conditions. Changing the parameters governing methane transport to reduce a specific methane emission pathway usually leads to an increased flux by the other two pathways as a compensation. However, large reductions in emissions through one pathway would lead to different emissions. This has two reasons: On one hand CH 4 oxidation is different for each pathway, on the other hand high CH 4 amounts can build up in the soil in case of limited transport capacities. For illustration Fig. 8b shows a comparison between the standard run and a run with a doubled tiller radius. As shown by Wania et al. (2010), an increase in the tiller radius leads to more plant-mediated transport and in turn to a lower ebullition flux (Fig. 8b), whereas the diffusive flux to the atmosphere remains practically unchanged. A larger tiller radius increases O 2 diffusion into the soil and, thus, potentially favours oxidation. At the same time, a larger tiller radius also favours a larger methane efflux to the atmosphere (Wania et al., 2010). We find that the fraction of produced methane that is oxidised in the soil increases from 30 to 32 % when the tiller radius is doubled.
It is interesting to analyse how the drop in emission for the simulated climate anomaly depends on the choice of parameters. The reduction in methane emissions after the freshwater perturbation barely depends on the oxidation fraction ( Table 1). The reduction in emissions is only 16.9 % when the oxidation fraction is set to 0.2 compared to a reduction by −18.6 % for the standard value of 0.5. Correspondingly, the reduction in emissions becomes slightly larger if a larger fraction of available oxygen is used for methane oxidation. Changes in the tiller radius and in the parameters governing ebullition also have a small influence on the simulated change in emission (Table 1).

Discussion and conclusions
Various climate input scenarios and parameter sets were used to model methane emissions during a freshwater experiment under interglacial climate conditions, reflecting conditions to some degree similar to the 8.2 kyr event. How well does the simulated climate pattern fit with reconstructions of the climate anomaly about 8.2 kyr BP?
The Greenland temperature in our input data drops by almost 4 • C over 20 yr, whereas Kobashi et al. (2007) estimate a cooling by 3.3 ± 1 • C during the 8.2 kyr event from icecore measurements, revising previous less accurate estimates of 6 • C (Alley et al., 1997) and 7.4 • C (5.4-11.7 • C; Leuenberger et al., 1999). Kobashi et al. (2007) suggest a drying, especially surrounding the North Atlantic area, and an average northern (30-90 • N) temperature cooling by 1-2 • C during the 8.2 kyr event. This agrees with findings from many paleo-data from Europe (Wiersma et al., 2006) and also fits with our input data (see Fig. 5). Morrill et al. (2005) present data from paleoclimate proxy records for the 8.2 kyr event from 52 site measurements in Greenland, Europe, America and Africa. The magnitude of their reconstructed temperature drops in Northern Germany, Estonia, Greenland and Norway compare very well to our input data. Also Alley et al. (2005) simulate a similar temperature change pattern to our input (Fig. 6), while LeGrande et al. (2008) find a bit lower surface temperature drops in North America (−0.6 • C), Middle Europe (−1.4 • C), Scandinavia and Siberia compared to our input data. Tindall et al. (2011) simulate maximum temperature anomalies of −1.33 • C in Siberia, −4 • C in Scandinavia and no changes in Middle Europe and wide parts of North America. Further, they find generally a slight precipitation decline of 17 % at maximum, but much less in the area where we simulate peatlands. The precipitation change in our input data is in most grid cells between 1 % and 5 %. In summary we can conclude that even though we use climate input data from simulations at pre-industrial conditions, the patterns and amplitudes of changes in temperature and precipitation are comparable to what is known about the 8.2 kyr event. However, the seasonal temperature and precipitation variability might be different in the input data compared to the variability at 8.2 kyr BP.
For the period 8.2 kyr ago, the location and initiation of boreal peatlands is reasonably well known from peat data (Yu et al., 2010). We used the present day boreal peatland distribution in our simulations based on a map by Tarnocai et al. (2009). This approach appears to be justified, because the comparison of the recent peatland distribution with the ice coverage 8.2 kyr ago (Peltier, 2004) shows that only very small areas of present day peatlands were covered by land ice at that time. Warm Holocene climate conditions prevailed already for more than 1000 yr at that time allowing for an initiation of peat expansion to interglacial conditions (Renssen et al., 2012;Seppä et al., 2009). A sensitivity test in a simulation with ice masking out peatlands 8.2 kyr ago only shows a reduction of 6 % in the simulated emissions and no change to the sensitivity at all. The presence of remnants of the Laurentide ice sheet, thus, appears to have only a minor effect on our boreal CH 4 emissions.
Beside climate input data and ice sheet configuration other uncertainties arise from the preindustrial conditions used in our simulations. The soil carbon inventory in our simulations is in steady state under preindustrial conditions and realistic compared to available data (Yu et al., 2010;Tarnocai et al., 2009). Therefore, it is higher than at the beginning of the 8.2 kyr event. In addition the different solar insolation and CO 2 -level may have an influence on the seasonality, the climate response to the freshwater forcing, as well as on the vegetation productivity. The effects on the latter are of opposite direction. While higher solar insolation at high latitudes favours plant growth during the 8.2 kyr event, higher CO 2 levels during the pre-industrial period have a larger fertilisation effect on plant productivity.
Knowing the above limitations, we can compare our results from the freshwater experiment to the 8.2 kyr event. We find that the decrease of 18 ppb in our simulation could explain about 23 % of the 80 ppb reduction during the 8.2 kyr event observed in ice core records (Fig. 1). This would suggest that variations in emissions from boreal peatlands contributed significantly to the CH 4 variations recorded in ice cores, but that other boreal and especially tropical CH 4 sources must have varied much stronger. Evidence for an attribution of reduced CH 4 emissions to tropical wetlands is suggested in Spahni et al. (2003). Part of the CH 4 concentration changes could also be explained by atmospheric sink changes. However, Levine et al. (2011) find that the change in CH 4 from the Last Glacial Maximum to the pre-industrial era was almost entirely source driven. It is, thus, not very likely that the smaller climate changes at the 8.2 kyr event lead to drastic sink changes.
We find that peatland CH 4 emissions drop by 19 %, caused by a reduction in precipitation of only 3 mm per month and a cooling of about 2.5 • C averaged over 50 yr. In addition the emission reduction is not very sensitive to the selection of most parameters in our model. The driving factor for the total magnitude of the modelled emissions is the production. Changes in the model parameters that describe the different transport mechanisms only influence the seasonal timing of the emissions and the attribution to the different transport pathways. The type and rate of CH 4 transport by diffusion in pore water or plants is important on whether CH 4 is directly moved to the atmosphere or gets stored in the soil. In the latter case it stays until it reaches the threshold for ebullition, which is the fastest and most efficient transport pathway. Therefore, the corresponding model parameters have an effect on the intraannual variability of CH 4 emissions.
On top of the variability of CH 4 transport and storage, the CH 4 production also changes significantly from year to year in our model. The modeled interannual variability in CH 4 emissions is on the order of 10 Tg CH 4 yr −1 (27-37 Tg CH 4 yr −1 ), which is comparable to the variability found by Chen and Prinn (2006) and Spahni et al. (2011). This interannual variability amounts to a third of the average peatland emissions of 30 Tg CH 4 yr −1 and is also comparable in size to the average and sustained decline in CH 4 emissions found during our freshwater hosing experiment.
It is, thus, plausible that a decline in peatland CH 4 emissions occurred during the 8.2 kyr event similar to the decline simulated during the freshwater hosing event in our model and that this decline was at least partly caused by a direct, fast response of the existing peatland ecosystems to climate change. The observed emission decline during the 8.2 kyr event occurred on time scales of decades, similarly as the simulated decrease in our experiments, and again showing that existing peatlands are able to respond quickly and sensitively to changing climate boundary conditions. Thus, we conclude that a sustained change in peatland emissions on the order of 10 Tg CH 4 yr −1 does not require a much slower climate induced geographical relocation of peatlands or a change in peatland ecosystem structure, such as a transgression from fens to bogs, but can be entirely explained by CH 4 emission changes that are comparable to current year-to-year emission variations.
A dynamic relocation or aging of peatlands is not yet included in our model, but would become important for the simulation of glacial CH 4 emissions or during transient Holocene climate conditions on time scales of millennia, where peatland location, structure and ecosystem composition can change. For example, during glacial times northern boreal wetlands were most likely displaced southwards compared to the Holocene in response to the widespread cooler conditions and the ice cover by the Laurentide and Fennoscandinavian ice sheets (Weber et al., 2010). For the 8.2 kyr event itself, which lasted only about 150 yr, a peatland relocation in response to the cooling appears to be unlikely as the time scale for peat development is much longer (Yu et al., 2010;Spahni et al., 2012).
The total average CH 4 emissions of 30 Tg CH 4 yr −1 are based on results from Spahni et al. (2011), compatible with Chen and Prinn (2006) for current conditions, but are at the lower end compared to the survey of Zhuang et al. (2004). Zhuang et al. (2004) summarised the recent literature and found that emission estimates for the pan-Arctic region from eleven studies ranged from 31 to 106 Tg CH 4 yr −1 , although the upper end of their range is incompatible with ice core observations on the interhemispheric CH 4 gradient (Brook et al., 2000;Baumgartner et al., 2012). MacDonald et al. (2006) also suggests that boreal wetland CH 4 emissions might have been higher at 8 kyr BP than today although peat extent was lower. The reason for this is that the emissions from minerotrophic fens that dominated early during peat formation are much higher than from the currently predominant ombrotrophic sphagnum bogs that only developed over time. They suggest that the northern peatland complex was likely at 20 % of its current areal extent at the end of the YD and expanded to about 50 % by 8 ka. On the basis of current estimates of overall CH 4 production from northern peatlands, they may have contributed 12 to 27 Tg CH 4 yr −1 by 8 ka. However, taking the higher CH 4 emission rates of early fens into account, MacDonald et al. (2006) speculate that the emissions before the 8.2 kyr event may have been also considerably larger. Accordingly this would then also imply a larger drop in CH 4 concentrations during the 8.2 kyr event. Brook et al. (2000) and Baumgartner et al. (2012) estimate the total CH 4 source strength in the Northern Hemisphere to be 60-70 Tg CH 4 yr −1 for the early Holocene. This estimate includes all extratropical CH 4 sources (e.g., ruminants, wet soils, thawing of submerged permafrost, thermokarst lakes etc.). If we set the annual average peatland emissions to 50 instead of 30 Tg CH 4 yr −1 , we could explain a 30 ppb change during our freshwater hosing experiment. This could still only explain about 40 % of the observed CH 4 decrease during the 8.2 kyr event. Accordingly, a substantial contribution of tropical and other boreal sources to this event remains most likely.
The simulated CH 4 emission reduction (19 %) is generally within the range of other estimates and ice core observations. Kobashi et al. (2007) find that total CH 4 emissions decreased by 15 ± 5 % during the 8.2 kyr event, while Walter et al. (2001b) and Christensen et al. (2003) find that a global temperature change of ± 1 • C leads to about 20 % changes in methane emissions.
Changes in CH 4 emissions are still subject to considerable uncertainties with respect to the total peatland CH 4 emissions as well as the processes on the soil and ecosystem level. To constrain the parameters further, additional information from field data and potentially from isotopic information are required. Moreover, CH 4 emissions of peatland ecosystems also vary over time. As outlined above the temporal transgression of peatlands from minerotrophic fens to ombrotrophic bogs is not yet included in our model. Sowers (2010) measured a decreasing trend in δ 13 C throughout the Holocene beginning at −46.4 ‰ 11 kyr BP and decreasing to −48.4 ‰ at 1 ka. He suggests that the 2 ‰ δ 13 C drop is likely to be a combination of increased CH 4 emissions from Arctic lakes and wetlands, an increase in the ratio of C 3 /C 4 plants in wetlands and an increase in methanogenic communities utilising the CO 2 reduction pathway as compared to the acetate fermentation pathway. In particular, the latter process would be directly related to the transgression of fens to bogs and could explain the δ 13 CH 4 shift observed in ice cores (Sowers, 2010) in the first half of the Holocene.
Generally, model improvements could be made by a carbon allocation to the different layers instead of an overall total number for HR given to the methane routine, a more sophisticated parameterisation of the soil anoxity that regulates the production, a changing ratio of CH 4 production and HR in transient simulations and by including a dynamic peatland distribution. The simulations done for this study could be improved by choosing more specific boundary conditions for the 8.2 kyr event (such as an improved peatland distribution map, the climate influence of the ice sheet remnants and the choice of orbital parameters and CO 2 concentrations at the time of the 8.2 kyr event) instead of preindustrial conditions as used in these simulations. There are also uncertainties to which extent the NCAR climate input is appropriate to do a 8.2 kyr event-like simulation. To improve the model calibration more site measurements would be needed and some further constraints like the 13 C signature would give additional constraints for the attribution to the different processes. However, despite these limitations, our model experiments and sensitivity tests indicate that it is most unlikely for the CH 4 drop observed in ice cores during the 8.2 kyr event to be solely due to changes in peatland CH 4 emissions.

Ebullition routine
We assume the existence of three gas species CH 4 , CO 2 and N 2 in the bubbles (Tokida et al., 2007). All calculations are done for each soil layer separately. The available volume per m 2 , V available , is the layer height multiplied by the porosity and corrected for the frozen volume.
The soil in wetlands is divided in the acrotelm (defined as the upper 0.3 m which experiences a fluctuating water table) and the catotelm, the underlying permanently inundated layers with 1.7 m thickness. This assumption is justified as the model simulates only for very few sites ever an acrotelm to catotelm boundary below 30 cm when this boundary is not fixed. The respective porosities are ρ acro = 0.9 for the three first layers, ρ cato = 0.8 for the 5 deeper layers. P E is the environmental pressure and, therefore, the sum of the atmospheric pressure (101 325 Pa) and the hydrostatic pressure. Initial conditions: n CO 2 = n CH 4 = 0 (A2) V gas = 1 %V available (A3) V diss,N 2 = 99 %V available (A4) We assume the gaseous volume of N 2 to be 1 % (see Sect. 2.3) of the available volume. This assures a minimum gas volume independent of the variables CH 4 and CO 2 . With Eq. (1) (ideal gas law and mass balance) we can calculate the total amount of N 2 Now we add CH 4 and CO 2 as they are simulated in the model.
n CO 2 (t) = n CO 2 (t − 1) + n CO 2 (A6) n CH 4 (t) = n CH 4 (t − 1) + n CH 4 (A7) n CH 4 (t) = n CH 4 (t − 1) + n CH 4 (A8) For the ebullition criterion we assume that it is critical if the gas volume exceeds 15 % of the available volume If P gas ≥ P E , ebullition is triggered: We assume that a fifteenth of the total gas volume bubbles out. The assumption is not critical (see Sect. 3.4) V = 1 %V available (A12) Therefore, the change in the total amount of CO 2 and CH 4 is ) while n N 2 is kept constant. If we weight each site with the time length of the measurement we get an improvement of 30 % for the site evaluation.