Interactive comment on “ Methane variations on orbital timescales : a transient modeling experiment ”

Abstract. Methane (CH4) variations on orbital timescales are often associated with variations in wetland coverage, most notably in the summer monsoon areas of the Northern Hemisphere. Here we test this assumption by simulating orbitally forced variations in global wetland emissions, using a simple wetland distribution and CH4 emissions model that has been run on the output of a climate model (CLIMBER-2) containing atmosphere, ocean and vegetation components. The transient climate modeling simulation extends over the last 650 000 yr and includes variations in land-ice distribution and greenhouse gases. Tropical temperature and global vegetation are found to be the dominant controls for global CH4 emissions and therefore atmospheric concentrations. The relative importance of wetland coverage, vegetation coverage, and emission temperatures depends on the specific climatic zone (boreal, tropics and Indian/Asian monsoon area) and timescale (precession, obliquity and glacial-interglacial timescales). Despite the low spatial resolution of the climate model and crude parameterizations for methane production and release, simulated variations in CH4 emissions agree well with those in measured concentrations, both in their time series and spectra. The simulated lags between emissions and orbital forcing also show close agreement with those found in measured data, both on the precession and obliquity timescale. We find causal links between atmospheric CH4 concentrations and tropical temperatures and global vegetation, but only covariance between monsoon precipitation and CH4 concentrations. The primary importance of the first two factors explains the lags found in the CH4 record from ice cores. Simulation of the dynamical vegetation response to climate variation on orbital timescales would be needed to reduce the uncertainty in these preliminary attributions.


Introduction
Methane (CH 4 ) is a forcing factor for climate as a greenhouse gas (Wang et al., 1996).Vice versa, CH 4 concentrations are also influenced by to climate changes, which cause changes in the natural CH 4 production rate.The most important natural source of CH 4 is anaerobic decomposition of organic material in wetlands (Fung et al., 1991).Global wetland extent, together with more local factors influencing CH 4 production and release to the atmosphere such as temperature and vegetation (Gedney et al., 2004), is therefore key to understanding variations in the atmospheric CH 4 concentration in the period before anthropogenic influences.
Variations in CH 4 concentration have been substantial during the late Quaternary.Spahni et al. (2005) found an increase in the Epica Dome C (EDC) ice core record of ∼360 parts per billion by volume (ppbv) to ∼690 ppbv during the transition from Last Glacial Maximum (LGM; 21 000 yr before present) to the Early Holocene (11 000 yr before present).Long-term oscillations in the CH 4 record have been linked to variations in insolation caused by orbital forcing, particularly to the 23 000 yr (23-Kyr) precession component (e.g.Ruddiman and Raymo, 2003).A strong 41 000 yr (41-Kyr) obliquity component and a 100-Kyr oscillation are also apparent in the record (Loulergue et al., 2008).These oscillations lag their respective forcing considerably.Ruddiman and Raymo (2003) devised an age model for the Vostok ice Published by Copernicus Publications on behalf of the European Geosciences Union.
core by tuning changes in CH 4 to be synchronous with low latitude mid-July insolation.This insolation curve, which is dominated by precession, effectively lags 21 June 65 • N insolation by ∼1.6 Kyr.They also point out that the CH 4 signal has the same delayed phase as northern ice sheets at the obliquity timescale.
The exact nature of the link between insolation and variations in past atmospheric CH 4 concentration is a muchdebated topic (Schmidt et al., 2004).The most common view on CH 4 variations at orbital timescales is that Northern Hemisphere (NH) extra-tropical (boreal; defined as the area north of 30 • N) emissions are mainly influenced by changes in temperature, while tropical (defined as the area between 30 • S and 30 • N) emissions are governed by changes in precipitation, i.e. monsoon intensity (Blunier et al., 1995).The land-sea temperature contrast that drives monsoon systems is sensitive to changes in insolation, in particular the 23-Kyr cycle which therefore greatly influences the strength of monsoon systems (Kutzbach et al., 2008).
Given that the bulk of modern day wetland emissions originates from the tropical zones with a modest contribution from the extra-tropical NH, and negligible Southern Hemisphere (SH) emissions (Walter et al., 2001;Prigent et al., 2007), it seems likely that changes in atmospheric CH 4 concentrations are primarily forced by changes in monsoon intensity (Chappellaz et al., 1990).However, Crowley (1991) noted that past atmospheric CH 4 concentrations compared just as favorably to reconstructed mid-latitude temperatures as they do to simulated Asian monsoon intensity.Therefore, changes in boreal emissions may well be the cause of the observed 23-Kyr oscillation in atmospheric CH 4 .
Modeling studies have recently investigated global wetland area and associated CH 4 emissions during the LGM (Kaplan, 2002;Valdes et al., 2005;Kaplan et al., 2006;Weber et al., 2011), thus implicitly focusing on the effects particular to a glacial termination.Here we perform for the first time a much longer simulation, which covers the last 650 000 yr, with the aim to examine the climatic origin of orbital cycles in atmospheric CH 4 .Focus is on the precession and obliquity cycles.The questions that we aim to address are as follows.First, whether the oscillations in CH 4 observed in ice core data originate from the tropics, the boreal zone and/or specifically from the Indian/Asian monsoon area; second, which climatic parameters play an important role in each of these regions; and third, how we can interpret the lags that have been found in the measured data with respect to the orbital forcing.
Wetland CH 4 emissions over the last 650 000 yr are derived using output from a transient climate simulation with an atmosphere-vegetation-ocean model of intermediate complexity (CLIMBER-2).The climate model has coarse spatial resolution.However, it adequately simulates large-scale climatic features like the seasonal migration of the Inter Tropical Convergence Zone (ITCZ), monsoon systems and the freezing and thawing of the northern land masses, which are of importance for wetland formation and global CH 4 emissions.The climate model output is used to force an offline model that computes the global distributions of wetland area and CH 4 emissions (TRENCH).Section 2 describes the models used.In Sect. 3 we analyze simulated emissions and examine the relative contribution to variations in global emissions of the tropics, the boreal zone, and the Indian/Asian monsoon area.The influence of changes in temperature, vegetation and hydrological processes is analyzed for each of the three regions.In Sect. 4 we discuss the results, compare model output to observational data, and draw conclusions.

The climate model
As a basis for simulating CH 4 emissions, the CLIMAte and BiosphERe model (CLIMBER-2, Petoukhov et al., 2000) was used.This is a coupled atmosphere-ocean-vegetation model of intermediate complexity that is suitable for long simulations due to its fast turnaround time.The model consists of a statistical-dynamical atmosphere, which resolves the large-scale flow that arises due to spatial temperature gradients.It parameterizes atmospheric transports caused by synoptic-scale variations (i.e.weather) in a sophisticated manner (Petoukhov et al., 2000), but it does not contain the weather events themselves.The other components are a 3basin zonally averaged ocean, including sea ice, and a terrestrial vegetation model VECODE (VEgetation COntinuous DEscription; Brovkin et al., 1997).The latter computes the fraction of the potential vegetation in a grid cell (i.e., grass, trees, and bare soil) as a function of the annual sum of positive day-temperatures and the annual precipitation.The atmospheric model of CLIMBER-2 has a resolution of 10 • in latitude and 51.43 • in longitude.The land surface model computes soil characteristics on the same grid (with the same resolution) as the atmospheric component.Soil moisture in CLIMBER-2 is based on a basic runoff scheme, taking the balance between precipitation and evaporation into account (Petoukhov et al., 2000).
With CLIMBER-2 we performed a climate simulation for the interval 650 Kyr BP to present, using orbital forcing (Laskar, 2004) as well as prescribing time-varying NH ice sheets and greenhouse gas concentrations (GHG).GHG concentrations (CO 2 and CH 4 ) were derived from EPICA Dome C, with the EDC-3 timescale (Parennin et al., 2007;Lüthi et al., 2008;Loulergue et al., 2008).For the last 500 yr we used values from several sources (Robertson et al., 2001).The total greenhouse gas forcing is dominated by CO 2 .Therefore, the simulated climate is only forced by the prescribed CH 4 concentration variations to a very limited extent.The volumes of the Laurentide and the Eurasian ice sheet were obtained from a 3-dimensional ice sheet model coupled to a model of deep-ocean temperature (Bintanja et al., 2005).The input for this model is formed by the marine benthic oxygen isotopes stack of Lisiecki and Raymo (2005).Their global isotope record was converted into a timeseries of ice sheet volumes using the LR04 time scale (Lisiecki and Raymo, 2005).The atmospheric model is forced by the ice area (surface albedo) and height (orographic effects) of ice sheets.Therefore, ice volumes were translated into a surface area and height of the Laurentide and Eurasian ice sheets based on the ICE-5G distribution for the LGM until present (Peltier, 2004), set on the spatial grid of CLIMBER-2.The correspondence between ice volume and area/height was extrapolated back in time in order to cover the whole period of 650 000 yr.The CLIMBER-2 simulation is described in detail by Weber and Tuenter (2011).This is the first simulation ever performed for the late Quaternary that takes variations in ice sheets and GHG concentrations into account.
Sea level effects due to changes in ice volume are not considered in this study.The increased capacity of newly exposed continental flats for wetlands could partially compensate for the loss of boreal wetland area due to ice cover.However, this effect is difficult to simulate within the present modeling set-up of CLIMBER-2, due to its coarse resolution.A previous model analysis (Weber et al., 2011) based on eight different General Circulation Model (GCM) simulations for the LGM, has shown that it is likely a minor factor in comparison to the main climate factors that are considered in the present study.We did test sea level effects in the present model set-up and performed a simulation where we identified areas with large continental flats (e.g.northern Siberia, the Sunda Shelf), and simulated CH 4 production from wetlands in these areas at times of low lea level.The main effects were a shift in source area of emissions from the boreal zone towards the tropics and a 7 % higher CH 4 production during the LGM and other glacials, in agreement with previous studies.The setup of this "sea level effect" run was speculative, however, mainly due to the low resolution of CLIMBER-2.Because the focus of the present paper is on analyzing climatic effects and the behavior of the different regions over time, we decided that the benefits of including this effect did not weigh up against the additional uncertainties.
Compared to GCMs, the main drawback of CLIMBER-2 is obviously its coarse resolution and associated parameterizations of small-scale dynamical processes.For the present purpose of estimating wetland CH 4 emissions over orbital timescales it is necessary to assess the model's sensitivity to the different forcings involved.These include both internal forcing mechanisms, such as ice sheets and greenhouse gases, as well as external orbital forcing.CLIMBER-2 is successful in simulating cold climates like the LGM (Ganopolski et al., 1998a), showing a similar sensitivity and response pattern as comprehensive GCMs to the separate impacts of large continental ice sheets and reduced greenhouse gases (Schneider von Deimling et al., 2006).Also, the strength of the NH monsoon during an orbitally forced period like the mid-Holocene lies within the range simu-lated by GCMs (Ganopolski et al., 1998b;Braconnot et al., 2007, Montoya et al., 2005).Limited model-data comparisons have been carried out for orbitally-forced variations in local proxy data representing the Indian monsoon (Ziegler et al., 2010a) and tropical temperature (Groot et al., 2011), which showed reasonable agreement both in timeseries and spectra.Leads and lags with respect to the orbital forcing (Tuenter et al., 2005) were found to be similar to those found in later simulations with more a complex model (Kutzbach et al., 2008).Many results found with intermediate-complexity models similar to CLIMBER-2 have been confirmed in a similar way by later simulations with GCMs (Weber, 2010).CLIMBER-2 thus seems capable of adequately representing climatic change in response to different forcings, albeit only at sub continental scale.

CH 4 emissions model
The wetland and CH 4 emission model TRENCH (TRansient Emissions of Natural CH 4 ) is driven off-line by 100-yr averaged monthly output from CLIMBER-2.The monthly output permits the analysis of seasonal variations in wetlands and CH 4 emissions.TRENCH consists of two consecutive components: a wetland finding algorithm and a CH 4 emission algorithm.The wetland location algorithm follows the approach originally proposed by Kaplan (2002), and uses climate variables from CLIMBER-2 that are known to influence wetland formation.Wetland existence is determined based on soil moisture and soil temperature.Wetlands are assumed to exist only during months when the soil temperature is above freezing point.When, in addition, a soil moisture threshold is reached, a grid cell is assumed to support wetlands.This threshold is set at 5 % of maximum saturation.This is low compared to other model studies, but one has to consider here the size of the grid cell over which soil moisture is computed which is much larger than in CLIMBER-2 than in GCMs.The fraction of a grid cell covered by wetland is taken to be linearly related to the soil moisture value derived from CLIMBER-2 as follows: Where C wet is the fraction of wetland in a grid cell, a is a tuning factor, M is the soil moisture, and I is the fraction of a grid cell covered with land ice.Resulting grid cell fractions are multiplied with the land surface area in a grid cell to obtain wetland area in km 2 .We assume that orography does not play a role at the present spatial resolution, as the large area of grid cells ensures that they always contain sufficient flat terrain to support wetlands when climatic conditions allow this.The threshold value of soil moisture was chosen by trial and error so as to obtain a reasonable geographical distribution of wetlands over the globe.Known wetland complexes such as the Hudson Bay lowlands and the present day Beringia area, as well as the obvious tropical belt and monsoon areas, are represented.Likewise, the Tibetan Plateau for example does not sustain wetlands in the model.We certainly miss smaller-scale wetlands associated with deltas, groundwater-fed systems etc. due to the basic nature of our land-surface model (compare the discussion in Kaplan, 2002).The tuning constant a was chosen to obtain a plausible value for worldwide wetland extent in km 2 .A comparison to values found in literature for the Pre-Industrial Holocene (PIH; 1850 AD) is made in Sect.2.3.The CH 4 emission algorithm of TRENCH is based on Gedney et al. (2004), who propose a simple model for emissions dependent on wetland fraction, soil carbon content, and temperature.The temperature sensitivity is expressed in the form of a Q 10 factor.In the present model we compute timevarying emissions for each grid cell as follows: Here E is the CH 4 emission (in Tg CH 4 ) of a given grid cell at a certain time step and month, k is a tunable constant chosen so that global annual emissions are close to measured data, C wet is the fraction of wetland in a grid cell as calculated in the wetland distribution model (Eq.1), and V a vegetation factor which is defined below.The temperature sensitivity is expressed in the form of a Q 10 factor, with T the soil temperature in degrees Celsius.
The Q 10 is a temperature-dependant variable itself, because the efficiency of CH 4 production changes with temperature (Gedney et al., 2004).In literature, values for Q 10 vary between 1.7-16 (e.g.Walter and Heimann, 2001) but most studies converge at values around 2-4.The base value Q 10 (T 0 )-the Q 10 at 0 • C -was set at 2.0.This value was chosen because it results in a reasonable distribution of CH 4 emissions over boreal wetlands (north of 30 • N) and tropical wetlands (latitudes between 30 • N and 30 • S), in agreement with previous observational and modeling studies.
Wetland CH 4 production also depends on the vegetation productivity determined by available nitrogen and soil carbon content, which is basically the amount of decomposable organic material in the substrate (Christensen et al., 2003).
Because the dynamical vegetation response to climate variations is only crudely simulated by the present version of CLIMBER-2, we used the grass and trees vegetation cover fractions F grass and F trees , provided by VECODE, as the best estimate.We added the two vegetation fractions together to form the vegetation factor V , where F trees carries 60 % of the weight and F grass 40 %.Trees have more extensive root systems than grass, and vegetation related root systems are important for the transport of CH 4 to the surface (Rice et al., 2010).As this is a questionable assumption, we tested different weighting factors and found that the results presented in Sect. 3 were not sensitive to the exact weighting of trees and grass.

Model performance for the PIH and LGM
The calculated maximum global wetland extent for the PIH occurs in June and is 6.3 10 6 km 2 .During minimum inundation, which occurs in December, the total wetland area is 3.6 10 6 km 2 .Taking into account that modern wetlands are possibly ∼20 % smaller than wetlands during the PIH due to destruction by agriculture (Chappellaz et al., 1993), this agrees with contemporary estimates based on satellite data of 5.9 10 6 km 2 and 2.1 10 6 km 2 , respectively (Prigent et al., 2007).Boreal wetlands (annual average) make out 38 % of the global annual mean wetland area, which is consistent with estimates for the present (Lehner and Döll, 2004;Prigent et al., 2007).Simulated wetland extent peaks in high northern latitudes (45 • N-65 • N), while extensive wetlands are also found in the tropics (Fig. 1).
Simulated global annual emissions are set at 151 Tg CH 4 for the PIH by the value taken for the tunable constant k.We chose this target value based on estimates of contemporary wetland emissions by Houweling et al. (2000) and Chen and Prinn (2006).This amount of 151 Tg CH 4 originates for 23 % (35 Tg) from boreal wetlands and 77 % (116 Tg) from the tropics.This distribution is in agreement with many previous observational and modeling studies, which find that 15-37 % of the global emissions comes from boreal wetlands and 56-85 % originates in tropical wetlands (Cao et al., 1996;Walter et al., 2001;Valdes et al., 2005).Boreal wetland area peaks during early summer when snowmelt rates surge, but the average temperature is still low (Fig. 2).Together with the limited emission season this makes the higher latitudes a less efficient source of CH 4 than the tropics.Therefore, the peak in simulated wetland extent in the high latitudes only translates to a moderate level of emissions (Fig. 1).TRENCH was thus tuned to reproduce the basic features of wetland distribution and CH 4 emissions for the PIH.We did not tune on any aspect of the model's transient behavior.
In order to assess model performance, the results for the LGM are now compared to those of earlier studies using GCMs.Global wetland area during the LGM is at its maximum in June with 5.3 10 6 km 2 and at its minimum in December with 3.3 10 6 km 2 , a decrease of 17 % in the annual-mean compared to the PIH.Boreal wetlands make up ∼26 % of the annual global mean, shifting the main component of wetland extent towards the tropics.This change of distribution is caused by the presence of large continental ice sheets at high latitudes and a relatively dry Eurasian mainland during the LGM due to generally reduced precipitation.Global CH 4 emissions during the LGM are reduced to 118 Tg CH 4 yr −1 (21 %) compared to the PIH.Tropical emissions are reduced by 12 %, while boreal emissions drop by 58 %.Previous GCM-based modeling studies have found similar moderate reductions in global emissions of 16-29 % (Kaplan, 2002;Valdes et al., 2005;Kaplan et al., 2006) as well as stronger reductions of 29-42 % (Weber et al., 2011).The decrease is especially strong in the boreal areas, which are active in the NH summer.The tropical zone is active year-round and is relatively unaffected or even positively affected considering wetlands during glacials.
During the LGM a large cooling occurs in higher latitudes.For example, Wu et al. (2007) reconstructed a 10 • C surface cooling in the Northern Hemisphere (NH), while only minor to moderate cooling of 0-3 • C takes place in the tropics.This is accurately simulated by CLIMBER-2 (Fig. 3).The strong cooling of the boreal zone limits the ability of large areas to sustain wetlands during most of the year due to the temperature constraint (Sect.2.2).This temperature-effect is larger than the direct loss of wetland area by ice cover.Additionally, the decreasing temperature has a negative impact on vegetation growth.Figure 3 shows that vegetation cover was strongly diminished in the LGM with respect to the PIH in boreal latitudes, while there is little change in the tropics.

Simulated emissions in the transient run
Figure 4 shows timeseries of the forcings used in the transient run, namely ice volume, greenhouse gas variations and insolation at 65 • N as an illustration of the orbital forcing.
We also show the global mean values of the main climatic parameters, as derived from CLIMBER-2, determining wetland extent, CH 4 emissions, and the emission temperature (defined as the average soil temperature of grid cells emitting CH 4 ), as well as the vegetation factor.All parameters show variations on multiple (orbital) timescales.Simulated temperature and vegetation changes, together with the varying wetland distribution, were used to compute total annual and global CH 4 emissions.These were found to vary between 116 and 167 Tg yr −1 .In Fig. 5 the global CH 4 emissions are compared to the atmospheric CH 4 concentrations derived from EDC.The EDC-3 curve shows pronounced sub-orbital variability, which is absent from the simulated curve.This is because CLIMBER-2 lacks internally generated variations on these timescales.Simulated climate and vegetation are purely driven by the applied forcings, which vary on orbital timescales.The longer-timescale variations (>3 Kyr) in both curves correspond very well, with a correlation of 0.79.As wetlands are assumed to be the largest natural source, wetland CH 4 emissions for a large part determine atmospheric concentrations during this period.Therefore we directly compare measured concentrations and simulated wetland emissions here, without taking lifetime variations into account.Spectra of the two time series were determined using Analyseries (Pailliard et al., 1996).Both spectra (Fig. 6) indicate a strong 100 Kyr signal and significant influence of 41 Kyr oscillations.Precession related features (i.e.23 Kyr) show a stronger signal in the model results than in the measured data.Differences in the relative importance of spectral peaks in measured concentrations and simulated global emissions could be due to lifetime variations that cannot be accounted for in the present model setup.As a check on our approach, the atmospheric CH 4 lifetime was tentatively diagnosed from the measured concentrations and simulated total wetland emissions.The modeled variations in wetland emissions are of the order of 30 % (50 Tg), which can be brought into agreement with the absolute variations in observed CH 4 concentrations of a few hundred ppb using reasonable assumptions on CH 4 lifetime variations between about 5.2 and 9.3 yr.For example, using a conversion factor between global burden (Tg CH 4 ) to surface mixing ratio (ppb) of 2.78 Tg ppb −1 (Denman et al., 2007) we find, for a PIH concentration of 676 ppb, a global burden of about 1880 Tg.Assuming a constant biomass burning source of 40 Tg CH 4 yr −1 (Fischer et al., 2008) and additional covarying sources, like ruminants, termites and geological sources, amounting to 30 Tg CH 4 yr −1 during the PIH (Denman et al., 2007), we find a lifetime of 1880/(150 + 70) = 8.6 yr for the PIH and that lifetimes vary between ∼5.2 yr during periods of low CH 4 concentrations, and ∼9.3 yr during periods of high CH 4 concentrations.This is consistent with data for the LGM and PIH (Fischer et al., 2008) and model studies that include atmospheric chemistry (Lelieveld et al., 1998).As these values are reasonable, our conclusion of close correspondence between the measured and simulated timeseries is supported.
In order to further analyze the origin of simulated wetland emissions we separately consider three regions (Fig. 7): the areas influenced by the Indian/Asian summer monsoon, the boreal area north of 30 • N, but excluding the East Asian monsoon grid box, and the tropical area, 30 • S-30 • N, excluding the remaining monsoon-dominated grid boxes.The Indian/Asian monsoon is analyzed separately as this area is generally assumed to be an important source of wetland emissions (Kutzbach et al., 2008).Moreover, this area is found to exhibit different behavior from the equatorial tropics where shifts in the Intertropical Convergence Zone (ITCZ) determine whether wetlands can exist or not.The bulk of global emissions (72 % during the PIH) are emitted from the tropics.The boreal zone and the monsoon areas contribute the remainder (18 % and 9 % respectively).The spectra for the individual regions (Fig. 6) show that they each have a typical timescale behavior, as mentioned above.Emissions from the Indian/Asian monsoon area are primarily controlled by precession, with little to no influence from 100 Kyr and 41 Kyr.Boreal emissions mainly reflect 41 Kyr and 100 Kyr signals, while tropical emissions contain a strong precession signal and significant power at the 41 Kyr and 100 Kyr frequencies.
The EDC-3 record for CH 4 and the time series for simulated global emissions, as well as those of the separate regions, were band pass filtered using Analyseries (Pailliard et al., 1996) at 21 (±2.5)Kyr and 41 (±3.5)Kyr.The filtered signals were then cross correlated to the astronomical precession index and obliquity parameter, respectively (Laskar, 2004).The measured EDC-3 derived CH 4 concentration is found to lag the orbital forcing by 2.0 Kyr on the precession timescale and by 5.5 Kyr on the obliquity timescale (Table 1).For the simulated global emissions lags are 1.3 Kyr and 5.3 Kyr, respectively, (Fig. 8).We consider this a very close match given that climate models are generally not successful in reproducing lags found in measured data (e.g.Ziegler et al., 2010b).Lags are longest for emissions from the boreal zone -reflecting the large influence of the ice sheets -and smallest for emissions from the monsoon areas, with intermediate lags for tropical emissions (Table 1).(Bintanja et al., 2005), insolation (Laskar, 2004), and GHG (Loulergue et al., 2008), and the resulting wetland extent, vegetation factor and mean emission temperature (lower graphs).

Factor analyses
The next step is to analyze the climate and vegetation factors that control variations in simulated CH 4 emissions and the simulated lags.First we define a total modification factor F tot as follows: Where E denotes annual emissions, t is time and E PIH denotes annual emissions during the PIH.According to Eq. ( 2 there are three individual factors determining the total emission changes, namely changes in wetland cover (F c ), in vegetation cover (F v ) and in temperature (F t ).To compute individual modification factors for each grid box, one term in Eq. ( 2) was set at its transient value and the other two terms were set at PIH values.The results were then divided by the emissions of that grid box for the PIH.Eq 5 shows this fac-torial calculation for wetland cover (F c ): Where k is the unaltered tunable constant, C(t) is the transient wetland cover, V PIH is the vegetation cover and Q T /10 10 PIH is the temperature dependence, both fixed at their PIH values.For each grid box F tot is equal to F c• F v• F t .The results from the factor analysis are not meaningful at the grid Table 1.Lags with respect to the orbital forcing in ice cover, greenhouse gases and EDC-3 atmospheric CH 4 , in the TRENCH simulated global and regional emissions and in the global modification factors.Lags are computed by searching for the optimum correlation between the forcing (precession or obliquity) and a given variable, band pass filtered at 21 Kyr or 41 Kyr, respectively.Zero phase is set at minimum precession and maximum obliquity, which each correspond to an insolation maximum for the NH.box size.Therefore, we only discuss the results as averaged globally and for the three separate zones defined above using a weighted average.By weighing the modification factors of each grid box according to their PIH emissions in the averaging process, skewing effects of extreme modifications in low-emission cells such as in high boreal latitudes are countered.The factor analysis is carried out over that part of the globe that could in principle support wetlands at all times, excluding areas that are covered by ice during some time steps.Focus is thus on climatic effects, not on the geographic controls (see also Weber et al., 2011).
Timeseries of the modification factors (Fig. 9) show that in our simulation temperature is the dominant control both for tropical and boreal emissions, with vegetation as a secondary control.During cold periods wetland extent is a reducing factor in the boreal zone, but an amplifying factor in the tropics, possibly related to reduced evapotranspiration.However, the effect is small in both zones and almost cancels out in the global mean.Wetlands coverage is only an important factor in the Indian/Asian monsoon region, together with vegetation changes.Spectral analyses (Fig. 10) show that precession-related variability is dominated by vegetation changes in the tropics and in the Indian/Asian monsoon region, while in the To summarize, at the precession timescale vegetation seems to be the dominant control of global CH 4 emissions, with temperature as a secondary control.The dominance of vegetation control originates primarily from the tropics.Our results also show a clear impact of vegetation on emissions in the Indian/Asian monsoon area, but these constitute only a small part of the global emissions.At the obliquity timescale temperature and vegetation play equally important roles.The simulated temperature changes are almost completely in phase with ice volume minima/greenhouse gas maxima.Vegetation changes lag summer insolation maxima by ca.1Kyr for precession to 2-4 Kyr for obliquity (Table 1).This reflects the strong influence of hydrological processes, which respond quasi-instantaneously to the orbital forcing, on vegetation.Compared to insolation, the lags of both temperature and vegetation effects tend to be longer in the boreal zone than in the tropics (Table 1).The combined effect of temperature and vegetation explains the lags that are found between simulated CH 4 emissions and the orbital forcings.

Summary, discussion and conclusions
We have simulated wetland CH 4 emissions over the last 650 000 yr using a simple wetland distribution and CH 4 emissions model coupled off-line to the atmosphere-oceanvegetation climate model CLIMBER-2.The resulting simulated global emissions show a close similarity to the measured EDC-3 timeseries of atmospheric CH 4 concentrations, both in spectra and in lags with respect to the orbital forcing.
In this study we find strong indications that global emissions are dominated by emissions from the tropics.However, in contrast to the general assertion that monsoon intensity is of primary importance (Chappellaz et al., 1990;Blunier et al., 1995) we find in our model that temperature and vegetation changes are the main drivers.As hypothesized by Crowley (1991), monsoon precipitation and global emissions do appear to co-vary, but no causal link was found.On the global scale, wetland extent was found to exert little control on emissions.Tropical temperature variations are found to be small, about 2-4 • C on glacial-interglacial timescales in agreement with proxy data (Farrera et al., 1999).However, such variations are significant as methanogenesis is much more sensitive to temperature changes at the high average temperatures in the tropics than at the lower average temperatures at high latitudes (compare Eq. 2).
Based on the present results we hypothesize that temperature and vegetation are the dominant controls of orbitalscale variations in global wetland CH 4 emissions, while wetland extent plays only a minor role.This is supported by a recent modeling study, based on an eight-member ensemble of GCM simulations for the LGM (Weber et al., 2011), where it was also found that temperature and vegetation are the main controlling factors.In their study wetlands exhibited southward shifts in the cold LGM climate, related to a southward shift of the mid-latitude storm tracks and the ITCZ.This caused the effect of changes in wetland distribution to be of regional importance only, and to cancel out in the global mean.Here we find a shift from the boreal to the tropical zone.This larger-scale pattern is likely related to the coarse resolution of CLIMBER-2 and the fact that the atmospheric dynamics are simplified, so that storm tracks are not resolved.This is a shortcoming of the present simulation.However, it does not seem to affect our main conclusion on the primary role of temperature and vegetation.Orbital cycles in vegetation have been found in many studies based on pollen data, for example for tropical South America (Hooghiemstra et al., 1993), Southeast Asia (Tsukada, 1996), Europe (Tzedakis et al., 2009) and North America (Whitlock and Bartlein, 1997).These studies indicate large shifts in vegetation presence and cover type.The vegetation model used in the present study only very crudely resolves such shifts.Nevertheless, we believe that our results on the role of vegetation are qualitatively valid although precise quantitative estimates are not possible using CLIMBER-2.
Finally, our results provide a plausible explanation for the lags found in the measured EDC-3 curve with respect to the orbital forcing.Atmospheric CH 4 concentrations have previously been related to monsoon precipitation (Ruddiman and Raymo, 2003).Here we find that lags in the CH 4 record are very likely due to temperature and vegetation effects rather than monsoon precipitation.The former two factors have a phasing in-between insolation maxima and ice sheet/greenhouse gas minima.This results in lags that are very close to those found in measured CH 4 data, both at the precession and obliquity timescale.In contrast, simulated monsoon precipitation for the late Quaternary has zero lag at the precession band (Weber and Tuenter, 2011).We conclude that process-based modeling of proxies is a necessary step to substantiate a reinterpretation of proxy records (Ziegler et al., 2010b;Clemens et al., 2010) and to resolve this mismatch between simulated and measured records.We stress here that our modeling approach is basic and needs to be followed up by more advanced modeling exercises (e.g.Singarayer et al., 2011), preferably with different models (Wolff, 2011).The attribution to the vegetation factor in our work suggests that improved simulations of the global dynamical vegetation response to climate variations on orbital timescales could provide a better understanding of the observed CH 4 concentration variations.Notwithstanding this caveat, we believe that our hypothesis of the dominance of vegetation and temperature effects may hold, as it seems a very basic outcome not crucially sensitive to the details of our modeling approach.

Fig. 1 .
Fig. 1.TRENCH simulated PIH wetlands during June, the month of maximum wetland extent (green line), and CH 4 emissions (black line), both integrated zonally and by 10 • latitude belts.

Fig. 2 .
Fig. 2. Seasonal variation of high (55-75 • N) latitude boreal wetland extent (black line) and emission temperature (defined as the average soil temperature of grid cells emitting CH 4 ; green line) during the PIH.

Fig. 7 .
Fig. 7. Definition of the tropical (red), boreal (blue) and Indian/Asian monsoon (pink) zones.Emissions from these three zones are separately analyzed.

Fig. 8 .
Fig. 8. (top) TRENCH simulated global CH 4 emissions, band pass filtered for 21 Kyr (±2.5 kyr) and the precession index, and (bottom) TRENCH simulated global CH 4 emissions, band pass filtered at 41 k (±3 Kyr), and the obliquity parameter.Lags between forcing and response are indicated for each timescale.

Fig. 10 .
Fig. 10.Spectra of the modification factors (see Sect. 3.2) for the entire globe, and for the tropical, boreal and Indian/Asian monsoon areas separately.
CLIMBER-2 simulated PIH temperature (solid black line) and LGM temperature (solid green line) for the NH summer months, and annual-mean grid-box vegetation factor (Sect. 2.2) during the PIH (dashed black line) and LGM (dashed green line).Note that the vegetation factor is a weighted average of tree cover (60 %) and grass cover (40 %).All data are integrated zonally and by 10 • latitude belts.