Terrestrial methane emissions from Last Glacial Maximum to preindustrial

We investigate the changes in terrestrial natural methane emissions between the Last Glacial Maximum (LGM) and preindustrial (PI) by performing time-slice experiments with a methane-enabled version of MPI-ESM, the Max Planck Institute for Meteorology Earth System Model. We consider all natural sources of methane except for emissions from wild animals and geological sources, i.e. emissions from wetlands, fires, and termites. Changes are dominated by changes in tropical wetland emissions, with mid-to-high latitude wetlands playing a secondary role, and all other natural sources being of minor importance. 5 The emissions are determined by the interplay of vegetation productivity, a function of CO2 and temperature, source area size, affected by sea level and ice sheet extent, and the state of the West African Monsoon, with increased emissions from north Africa during strong monsoon phases. We show that it is possible to explain the difference in atmospheric methane between LGM and PI purely by changes in emissions. As emissions more than double between LGM and PI, changes in the atmospheric lifetime of CH4, as proposed in 10 other studies, are not required.


Introduction
The atmospheric concentration of methane undergoes major changes in the time between the last glacial maximum (LGM) and preindustrial (PI). Between LGM and 10 ka BP (before-present, with present = 1950 CE) atmospheric CH 4 , as reconstructed from ice cores, nearly doubles from~380 ppb at LGM to 695 ppb at 10 ka BP (Köhler et al., 2017), with very rapid con- 15 centration changes of about 150 ppb occurring during the transitions from the Bølling Allerød (BA) into the Younger Dryas (YD) and from the YD into the Preboreal (PB) / early Holocene ( Figure 1). Furthermore, while Holocene atmospheric CH 4 is very similar for 10 ka BP and PI (694 ppb, mean concentration for 300 a BP to 200 a BP), CH 4 decreases linearly by 15% at from 10 to 5 ka BP and increases again linearly towards PI. If we assume that the atmospheric lifetime of CH 4 did not change dramatically between the LGM and the present, these changes in atmospheric CH 4 would require large changes in CH 4 20 emissions.
The change in methane between LGM and PI has been investigated in a number of studies. Some have used box models to explain the methane changes observed in ice cores. Recently Bock et al. (2017), for example, pointed to tropical wetlands as the main driver of glacial-interglacial CH 4 change from a study of methane isotopes from ice cores. In addition there are studies with comprehensive models. Kaplan (2002) investigated wetland CH 4 emissions during the LGM and the present geological sources, wild animals, wildfires, permafrost soils and vegetation as further onshore (land) sources and geological and "other" as offshore (oceanic) sources.
We aim to assess the changes in the natural sources of methane from the LGM to the present in order to determine the factors driving the changes in atmospheric CH 4 . We use a methane-enabled version of MPI-ESM, the Max Planck Institute Earth System Model, to investigate changes in natural methane emissions for six time slices from the LGM to the present. In this model we include submodels for methane fluxes from wetlands, termites and wildfires, but of the other natural methane 65 fluxes listed above, many cannot easily be derived from the climate model state and therefore are neglected for now. As time slice experiments very likely are not helpful for looking into the BA-YD and YD-PB transitions, we neglect these for now, focusing instead on the longer timescale changes in methane.

Model and experiments
2.1 MPI-ESM 1.2 70 We use the Max Planck Institute Earth System Model (MPI-ESM) in version 1. 2 (Mauritsen et al., 2019), the version to be used in CMIP6. All experiments are performed in resolution T31GR30 . In comparison to the CMIP5 version , a number of errors were corrected in the atmosphere and ocean models, and the land surface scheme JSBACH Brovkin et al., 2013;Schneck et al., 2013) has been updated with a multilayer hydrology scheme (Hagemann and Stacke, 2015), the SPITFIRE fire model (Thonicke et al., 2010;Lasslop et al., 2014), and the improved 75 soil carbon model YASSO (Tuomi et al., 2009;Goll et al., 2015).

Wetland methane emission model
The present-day area that wetland methane emissions originate from is highly uncertain. The generation of methane in the soil is dependent on plant composition, carbon content and carbon quality, essentially ecosystem properties, as well as the degree of anoxia in the soil, which depends on soil structure and water content, essentially hydrological properties. As there is no 80 better estimate of the methane-generating area available, we determine the surface inundation and assume that this is a useful approximation of the areas where methane is generated.

Dynamic inundation model
We use an approach based on the TOPMODEL hydrological framework (Beven and Kirkby, 1979) to determine inundation extent dynamically. TOPMODEL is a conceptual rainfall-runoff model, based on the compound topographic index (CTI) 85 χ i = ln (α i / tan (β i )) with α i a dimensionless index for the area draining through point i and β i the local slope at that point.
TOPMODEL determines the local water table z i in point i in relation to the grid cell mean water tablez: with χ i the local CTI index in point i,χ the grid cell mean CTI index, and f a parameter describing the exponential decline of transmissivity with depth. From Eq. 1 we determine the grid cell fraction with a local water table depth z i ≥ 0. Since inundated 90 areas become unreasonably large in some locations, we limit the valid range of CTI values by introducing the constraint χ i ≥ χ min following Stocker et al. (2014), with χ min constant in space and time. We assume this to be the inundated and therefore methane-emitting area A inun , unless soils are frozen. In these cases we determine the fraction of liquid water in the soil f liq from the soil temperature T soil (limiting f liq to 0.1 ≤ f liq ≤ 1 for numerical reasons), and determine the inundated area as A inun = A inun × f liq , reducing the inundated area under freezing conditions, as frozen soils emit less methane, similar to Gedney and Cox (2003).
To determine the grid cell mean water table position z, we determine the layer saturation Ψ k = Θ k /Θ f c for each soil layer k by dividing the volumetric moisture content Θ k by the field capacity Θ f c . Starting from the bottom of the soil column, z is located in the first soil layer l with layer saturation Ψ l less than the saturation threshold Ψ thres . The final water table position 100 then is with z b,l the bottom of soil layer l, and ∆z l the height of soil layer l.
Values for f , χ min , and Ψ thres were determined from sensitivity experiments. In the experiments described here, we use f = 2.6 , χ min = 8.5, and Ψ thres = 0.95. Furthermore, comparison with remote sensing data (Prigent et al., 2012) showed that 105 inundated area in grid cells with a mean CTI indexχ ≤ 5.5 is negligible. Inundation is therefore only determined for grid cells withχ > 5.5.
We use the CTI index product by Marthews et al. (2015) for the CTI index at a resolution of 15 in all present-day land areas, while we determine CTI index values for shelf areas that are below sea level at present, but above sea level under glacial conditions, from the ETOPO1 dataset (Amante and Eakins, 2009) using the topmodel library for R (Buytaert, 2011). In order 110 to reduce storage requirements, we approximate the distribution of CTI values within a model grid cell by fitting a gamma distribution, following Sivapalan et al. (1987).

Wetland methane production and transport
We use the methane transport model by Riley et al. (2011)  Adaptations are described in the following. In the grid cell fraction determined to be inundated by the inundation model, soil organic matter (SOM) is decomposed under anaerobic conditions in the YASSO soil carbon model , assuming a reduction of decomposition by a factor of 0.35 (Wania et al., 2010) in comparison to the aerobic case. As YASSO 120 is a zero-dimensional representation of soil C processes, we distribute the decomposition product to the soil layers according to the root distribution from Jackson et al. (1996). Partitioning of the anaerobic decomposition product into CO 2 and CH 4 is temperature-dependent, as in the original Riley model, with a baseline fraction of CH 4 production fCH 4 = 0.35 and a Q 10 factor for fCH 4 of Q 10 = 1.8 with a reference temperature of 295K.
For each grid cell, the methane model determines CH 4 production and transport for two grid cell fractions, the aerobic (non-125 inundated) and the anaerobic (inundated) fraction of the grid cell. If the inundated fraction changes, the amounts of CO 2 , CH 4 and O 2 are conserved, transferring gases from the shrinking fraction to the growing fraction, proportional to the area change.
While vegetation in JSBACH is determined for vegetation tiles, allowing a fractional coverage of plant functional types, the relevant properties (root distribution, SOM decomposition) are aggregated to grid cell level for the methane transport model for performance reasons. Previous sensitivity experiments showed that differences to a tile-resolving formulation are negligible.

Methane emissions from wildfires
To determine methane emissions from wildfires, we use the biomass burned, diagnosed from the SPITFIRE fire module (Thonicke et al., 2010;Lasslop et al., 2014), as well as information on vegetation composition from the dynamical vegetation model.
We use the methane emission factors from Kaiser et al. (2012), mapped to the JSBACH plant functional types, to determine the fraction of burned biomass emitted as methane. Therefore changes in fire-related methane emissions are completely deter-135 mined by changes in fire carbon emissions. Fire occurrence in the SPITFIRE model is determined as a function of flammability (higher under dryer/warmer conditions) and ignition probability, with ignition probability a function of lightning frequency and population density. We are currently limited to a fixed lightning distribution reflecting modern conditions, and we are assuming a population density of zero for all time slices earlier than preindustrial. Therefore the main factors affecting fire-related methane emissions are carbon content and moisture conditions.

140
For PI and PD we use an estimate of population density to determine the ignition probability, with ignition probability increasing with population density. However under very high population densities it is assumed that fire suppression increases, thus decreasing fire probability, and thus fire methane emissions, for very high population densities.

Methane emissions from termites
Methane emissions from termites are determined following the approach developed by Kirschke et al. (2013) and elaborated by tropical broadleaf evergreen tree, tropical broadleaf deciduous tree, and C4 grass, we determine M termite from the annual gross primary production GPP using M termite = 1.21 × exp (GPP × 0.0008) (Kirschke et al., 2013). From M termite we determine

Model experiments
We  The time slice experiments were initialised from a -so far unpublished -transient model experiment from 26 ka BP to PI with prescribed ice sheet extent from the GLAC-1D ice sheet reconstruction (Tarasov et al., 2012;Briggs et al., 2014;Ivanovic et al., 2016). This model experiment was initiated at 26 ka BP and run transiently from then to PI, i.e. the year 1850 CE. Ice sheet extent, as well as bathymetry and topography (Meccia and Mikolajewicz, 2018) and river routing (Riddick et al., 2018) were continuously updated throughout the deglaciation. For the assessment of wetland methane emissions, the wetland area can to some extent be measured directly from satellites.

195
Remote-sensing products of surface inundation are available, for example by Prigent et al. (2012) and Schroeder et al. (2015).
To assess the quality of the modelled surface inundation, we rely on the Prigent et al. (2012) data set. However, four points need to be considered when comparing these data to model results: 1. The remote-sensing process is unable to penetrate snow cover, so snow-covered areas are considered non-inundated.
2. The remote-sensing product shows all inundated areas, including areas flooded as a result of anthropogenic processes, 200 such as the creation of reservoirs and rice-paddies, which are not considered in the model.
3. Remote-sensing may be unable to penetrate dense forest canopy, implying that inundation estimates may be biased in forested areas.
4. Not all methane-generating areas have water tables above the surface. Water tables in northern peatlands, for example, tend to be below the surface for part of the year, especially in the summer.

205
In order to make model output and remote-sensing data comparable, we therefore mask all snow-covered areas in the model output, and we use data on rice-growing areas by Monfreda et al. (2008) to mask all rice-growing areas from both the remotesensing data and the model output. After these modifications, modelled inundated areas for the present-day period (mean over  (Fig. 2). For the glacier-free NH extratropics (NXT, here defined as north of 30°N), the seasonality of inundation is similar in observations and model, but the summer peak in inundation is larger in the model (2.5 × 10 6 km 2 for the JJA mean) than in the observations (2.3 × 10 6 km 2 ). major differences are apparent: 1) The annual maximum inundation is more localised in the observations, while it is less clearly defined and reaching lower maximum values in the model, and 2) after removal of the rice-growing areas the model does not show a significant inundation maximum in India, very likely due to an underestimate of the Indian monsoon precipitation in the low model resolution. We thus judge the methane generating areas produced by the model as reasonable, keeping in mind the likely low bias of the remote-sensing inundation product.

220
As described above, we use the inundated area to determine the methane emitting area. To evaluate the inundated areas leading to the wetland emissions, it has to be kept in mind that NXT emissions mainly are from the summer season, implying that the JJA (June, July, August) mean inundation is relevant for these, while the seasonality of TRO emissions is much less pronounced, implying that the annual mean inundation is relevant. In the following we therefore assess the effective inundated area, defined as the annual mean inundated area in tropical latitudes (TRO, between 30°N and 30°S), and the JJA (June, July,

225
August) mean inundated area in the glacier-free NH extratropics (NXT, north of 30°N). For the present-day climate state, the effective inundated area is 1.5 × 10 6 km 2 in TRO and 2.6 × 10 6 km 2 in NXT (differences to the numbers shown above due to the removal of rice-growing areas in the comparison to observations).

Natural methane emissions
So far it has not been possible to directly measure the quantity -surface methane fluxes -that we aim to assess in this 230 publication on appropriate scales. Methane flux measurements exist for single sites of meter scale, mainly using measurement chambers, and for slightly larger scales, using eddy-covariance towers, but so far the scales relevant for global scale modelling, the model grid-cell to global scales, have not been covered by direct methane flux measurements (Melton et al., 2013;Saunois et al., 2016;Poulter et al., 2017). For assessment of our model experiments we therefore need to rely on global assessments (Saunois et al., 2016), and we can gain some additional insight from atmospheric inversions (Bousquet et al., 2011). Furthermore, wetland methane emission estimates from atmospheric inversions (Bousquet et al., 2011) show that the ma- dinal distribution of modelled PD wetland methane emissions therefore is well within the range obtained from atmospheric inversions.

250
Overall, the PD state is rather similar to the PI state assessed in the following section, with very small differences in the distribution of emissions ( Figures A1 and A2), but generally higher methane emissions. At 287.4 K, the global mean temperature in the PD climate state is 0.5 K warmer than preindustrial (Table 1). Precipitation is similar, leading to negligible differences in the effective inundation. With 1140 PgC, 8% larger than PI, the global soil C stock is also rather similar. However, vegetation productivity is enhanced in comparison to PI, due to warmer temperatures and higher CO 2 concentrations. The net methane 255 emissions in PD climate are 29% larger than PI (Table 2), with wetland methane emissions 33% larger, with a larger increase (+42%) in TRO than in NXT (+16%). Fire emissions are 18% larger than PI, termite emissions increase by 66%, and the soil uptake increases by 140%. The latter increase is largely due to the higher atmospheric concentration of CH 4 , which drives additional methane into the soils in comparison to the lower-CH 4 PI state. The larger fire emissions are mainly due to higher population densities in the 2000s than in 1850, although the very high population densities in eastern North America, Europe, 260 and southern Asia are assumed to drive an increase in fire suppression in the SPITFIRE model (Lasslop et al., 2014). Thus fire emissions are decreased here, despite the general increase in most other places. Termite emissions are higher in the modern climate due to an increase in GPP under higher CO 2 , while wetland emissions largely increase due to the higher temperatures of the modern climate, with CO 2 -fertilisation playing an additional role.

Preindustrial methane emissions 265
The climate in our PI experiment is very similar to the one described by Mikolajewicz et al. (2018). The global mean nearsurface air temperature is 286.9 K ( Table 1). The annual mean temperature in the TRO area T TRO is 294.5 K, while T NXT , the annual mean temperature in the NXT area, is 275.2 K. The NH ice sheet area is limited to Greenland, with the ice sheet having an area of 1.8 × 10 6 km 2 in our model setup. Under these climatic boundary conditions, we obtain total net methane emissions of 181 T gCH 4 yr −1 (Table 2), with wetlands contributing 167 T gCH 4 yr −1 (Fig. 4a), fire and termites 15 T gCH 4 yr −1 270 and 7.0 T gCH 4 yr −1 , respectively ( Fig. 5a and b), while the soil uptake is 7.3 T gCH 4 yr −1 (Fig. 5c).
Wetland emissions, the dominant natural component of the terrestrial methane fluxes, mainly originate in TRO (110 T gCH 4 yr −1 ), while emissions from NXT are 56 T gCH 4 yr −1 ( Table 2). The main factors determining wetland methane fluxes, apart from temperature, are the emitting area and the soil carbon stock that the soil respiration (and thus methane production) is derived from. In the PI state the global effective inundated area is 4.0 × 10 6 km 2 (Fig. 4b), of which 2.7 × 10 6 km 2 are located in NXT, 275 while 1.3 × 10 6 km 2 are in TRO (Table 1). For soil carbon, on the other hand, the global stock is 1054 PgC (1Pg = 10 15 g), with most of the soil carbon (588 PgC) located in NXT, while it is 439 PgC in TRO.
Methane emissions from fires (Fig. 5a) closely follow the fire distribution, with most fire methane emissions coming from subtropical Africa and South America, although some emissions also originate in North America, Southern Europe and South Asia. Termite emissions, on the other hand mainly originate from tropical regions, especially southern Asia, with minor con-280 tributions from subtropical regions on all continents (Fig. 5b). Methane uptake by upland soils (Fig. 5c), finally, is distributed widely with no large regional variations.

Wetland methane emissions
Under LGM boundary conditions the global mean temperature is 4.4 K colder than under PI conditions (Table 1). Extensive glaciers cover the NH extratropics and sea level is lower, leading to a 15% increase in total land area, although total glacier-free 285 area is nearly identical (  Figure 6 a, Figure A3 a), while NXT effective inundated area I NXT decreases by 6%.
The global soil C stock is 617 PgC, substantially smaller (-41%) than at PI, with the decrease smaller in TRO (-33%) than in 290 NXT (-49%). As a result of these climate changes, wetland methane emissions decrease by 51% (Table 2, Figure 6 e, Figure   A3 e), with a TRO emission decrease of 47%, while NXT emissions decrease by 59%, with the majority of the latter emissions coming from areas in East Asia adjacent to the Yellow Sea and North America south of the Laurentide ice sheet. The wetland CH 4 emissions therefore decrease nearly everywhere (Figure 6 e), with one major exception: The continental shelf areas that are exposed due to the lower sea level become significant sources of methane, especially in Indonesia, but also in Africa and

295
Asia. At 29 T gCH 4 yr −1 , they contribute about 35% of the total wetland CH 4 emissions at the LGM.
For 15 ka BP, the global mean temperature change is −2.8 K, relative to PI. NH ice sheet extent is 24% lower than at LGM, but still extensive. The total land area is 12% larger than PI due to the lower sea level, but A nonglac in only larger by 1%.
A NXT is thus reduced by 10% (Table 1), while A TRO is increased by 11%. The change in T TRO is −2.0 K (Table 1), while it is −3.8 K for T NXT . Precipitation decreases by 6% in the global mean, with a a 4% decrease in TRO and a 10% decrease in 300 NXT. Precipitation in NH Africa is slightly increased due to a stronger West African monsoon. I TRO thus is larger by 29%, while I NXT is 8% smaller than PI (Table 1, Figure 6 b, Figure A3 b). Global soil C is at 815 PgC (-23%), with TRO C stocks at 398 PgC (-9%), while NXT stocks are at 392 PgC (-33%). Total wetland methane emissions decrease by -22% as a result, with TRO emissions decreasing by 17% and NXT emissions of by 33% ( Table 2, Figure 6 f, Figure A3 f). In contrast to the LGM situation, there is an increase in CH 4 emissions from NH (sub-) tropical Africa Figure 6 f) to 19 T gCH 4 yr −1 (+58%), due to 305 wetter conditions in the Sahel area. The exposed shelf areas emit about 38 T gCH 4 yr −1 overall, 29% of the total emissions.
For 10 ka BP, our model indicates a global mean temperature change of −0.7 K (Table 1). Glacial area is much reduced in comparison to the LGM, but remains of the Laurentide ice sheet still cover parts of north-eastern Canada, leading to a lower sea level than at PI. The total land area thus is 3% larger than at PI (A nonglac -1%), with A TRO 3% larger due to lower sea level and A NXT 4% smaller due to the remaining ice sheet coverage. The temperature decrease is larger in TRO (−1.0 K) than in 310 NXT (−0.6 K). Precipitation is near PI levels in the global mean (-1%), with a 6% increase in TRO, and a 1% decrease in NXT.
Precipitation in NH Africa is strongly increased due to a strong West African monsoon. I TRO is increased by 35%, mainly in due to the wetter conditions in north Africa, while I NXT is decreased by 12% in NXT (Table 1, Figure 6 c, Figure A3 c).
Global soil C is at 983 PgC, quite near the PI total stock (-7%), with a TRO soil C stock of 449 PgC (+2%) and a NXT stock of 510 PgC (-13%). As a result, wetland CH 4 emissions are very similar to PI (Table 2), with +7% in TRO wetland emissions area, where the West African monsoon is strongly increased, leading to more precipitation, increased inundated area, and more biomass and soil C. Emissions from NH Africa are 37 T gCH 4 yr −1 (+208%), an increase larger than the total increase in TRO emissions. Emissions from the (small) exposed shelf areas are at 8 T gCH 4 yr −1 (5% of the total wetland CH 4 emissions). NXT emissions are smaller than PI in North America and Europe, but they are larger than PI in northern Asia, due to the summer 320 warming from the changed insolation at 10 ka BP.
At 5 ka BP, global mean temperature change is at −0.2 K (Table 1). Ice sheet areas are as in the PI state, thus TRO and NXT areas are unchanged. T TRO is slightly lower (−0.6 K), but T NXT is very similar to PI. Precipitation changes are very small in the global mean, with a a 4% increase in TRO and a 2% increase in NXT, with the West African monsoon slightly stronger than PI. I TRO increases by 15%, while I NXT decreases by 1% (Table 1, Figure 6 d, Figure A3 d). Global soil C stocks are 1043 PgC, 325 slightly smaller than preindustrial (-1%), with an increase by 4% in TRO, especially in the southern Sahel region, while NXT is 4% lower than PI. Total wetland methane emissions increase by 2%, with TRO and NXT wetland emissions both increasing by 2% (

Methane emissions from wildfires
For all time slices before PI we assume that no humans were present, leading to a generally decreased probability of fire ignition in comparison to PI and PD. For the LGM (Fig. 7 a) fire CH 4 emissions ( Table 2) are 73% smaller. As biomass is reduced strongly under the cold and low-CO 2 conditions of the LGM, fire-related C emissions are also reduced. At 15 ka BP (Fig. 7 b) fire emissions are 53% lower than PI, while they are 40% smaller at 10 ka BP (Fig. 7 c). Generally, the spatial 335 pattern of emission changes at 15 and 10 ka BP mainly reflects precipitation changes: Enhanced emissions occur in areas where precipitation is reduced, enhancing vegetation flammability. In the Sahel area, this relationship is different, though. Here, the enhanced rainfall leads to an increase in vegetation cover, especially grass cover. As a result, more biomass is available for combustion, leading to enhanced emissions. At 5 ka BP, finally, fire emissions are reduced by 45%. As climate is already relatively similar to the PI situation, the main reason for the fire emission reduction here is the smaller ignition probability due 340 to the absence of humans.

Methane emissions from termites
Termite emissions are mainly determined by gross primary productivity (GPP) in tropical and subtropical areas. The lower atmospheric CO 2 and temperature under LGM conditions decrease GPP everywhere. Therefore termite CH 4 emissions are reduced by 58% relative to the PI level (Table 2). For 15 ka BP, there also is a general reduction in termite emissions, with 345 4.6 T gCH 4 yr −1 in total (-34%). However, the enhanced rainfall in the Sahel area leads to an increase in termite methane in this area. The latter is similar at 10 ka BP, where total emissions are 13% smaller in comparison to PI. The enhanced productivity in the Sahel therefore more than compensates the decrease in termite methane from the Amazon and African rain forests. At 5 ka BP, finally, termite emissions are slightly smaller than PI (-7%), with minor decreases in the rain forest areas and a slight increase in the Sahel. can transport gases between plant roots and leaves.
In our experiments, we find that the soil uptake of methane is to a large extent determined by the gradient of methane between soil and atmosphere. Thus higher atmospheric concentrations of methane directly lead to a larger soil uptake of methane. Under LGM conditions, the atmospheric CH 4 concentration is 370 ppb, slightly less than half the PI concentration. Consequently soil methane uptake decreases by 68% compared to PI (Table 2), with decreased temperatures being an additional factor (Fig. 9 a).

360
At 15 ka BP (atmospheric CH 4 of 464 ppb), the soil uptake is 52% smaller, while it is changed by -14% at 10 ka BP (688 ppb) and -28% at 5 ka BP (579 ppb). Spatially, the change in methane uptake is rather uniform, showing a similar reduction in uptake in most locations (Fig. 9). The exception to this is, once again, the Sahel area, which shows an increase in methane uptake most pronounced for 10 ka BP ( Fig. 9 c), but also for 5 ka BP (Fig. 9 d). For these time slices the increase in vegetation cover in the Sahel region leads to a localised increase in methane uptake. of net emissions, and soil uptake reduces the emissions by between 2.5% at 15 ka BP and 7.0% in the PD state (Table 2).
In the modelled emissions, we are missing two components of the natural methane cycle: Wild animals and geolog-  (Saunois et al., 2016), and estimates for other time slices are even less confident, but might be of the order of 15 − 20 T gCH 4 yr −1 for times before significant human influence (Chappellaz et al., 1993). In principle, these emissions should somehow be related to the net primary productivity, as this would determine the carrying capacity of the ecosystem, implying smaller fluxes in the glacial than in the Holocene. Adopting the ice-core based estimate by Petrenko et al. (2017) for the geological fluxes, we can thus hypothesise these to be 9 ± 9 T gCH 4 yr −1 , while wild animals might add 15 ± 10 T gCH 4 yr −1 . The total unaccounted fluxes might therefore be of the order of 24 ± 19 T gCH 4 yr −1 .

385
To compare the net fluxes to the reconstructed atmospheric CH 4 concentrations from ice cores, we determined the implied methane emissions (Fig. 10). We converted the methane concentrations into a methane burden, using a conversion factor of 2.767 T gCH 4 ppb −1 (Dlugokencky et al., 1998). With a tropospheric lifetime of 9.3 yrs (range given 7.1 − 10.6 yrs), an approximation for the present-day situation (Saunois et al., 2016), we then determined the methane flux required to obtain the CH 4 concentration reconstructed for all time slices except the present-day period, which is dominated by anthropogenic 390 emissions. From this flux we subtracted the unaccounted sources, as described above, to determine the implied emissions.
Uncertainties from the unaccounted fluxes are represented as error bars, however the uncertainty in tropospheric lifetime is not considered here, but would roughly add another 15%.
Comparing the modelled net emissions to the implied emissions (Fig. 10), the modelled fluxes are within the range of uncertainty for all time slices except for 15 ka BP and 5 ka BP, with modelled net emissions larger than the implied fluxes for 395 these time slices. The net emissions increase by more than 100% going from 20 ka BP to 10 ka BP and PI, we can thus explain the methane increase from LGM to Holocene with CH 4 emissions only, not requiring changes in methane lifetime. However, we so far cannot explain the Holocene changes in atmospheric CH 4 , decreasing between 10 ka BP and 5 ka BP, and increasing subsequently. We assume that this is due to an overestimate of the West African monsoon and its impact on African methane emissions, as a general reduction in the West African monsoon would lead to decreases in TRO emissions for 15 ka BP, 10 ka 400 BP, and 5 ka BP, bringing model results more in line with the implied emissions determined from ice core CH 4 . However, this is speculative at this point and would require further experiments.

Conclusions
In this assessment we considered all natural emissions of methane, with the exception of emissions from wild animals and geological sources. In our experiments we found that it is possible to explain the difference between LGM (20 ka BP) and PI 1. Global mean temperature and CO 2 : Higher atmospheric CO 2 concentrations increase NPP and thus also soil carbon available for anaerobic decomposition to CH 4 . Similarly, higher global mean temperature also increases NPP and soil 410 C decomposition, and it furthermore increases the ratio of CH 4 to CO 2 production in anaerobic decomposition. Thus, higher atmospheric CO 2 and higher global mean temperature lead to larger wetland emissions of CH 4 . This affects emissions from wildfires and termites in a similar way, as fire C release is dependent on biomass and termite biomass is dependent on GPP and thus CO 2 and temperature.
are covered by the ice sheets, which is especially important for wetland methane emissions from North America. At the same time large ice sheets lower sea level, enlarging tropical wetland area as the continental shelf is exposed and becomes a significant source of methane. This is mainly relevant in the tropics, as high latitude shelf areas exposed under glacial conditions, for example the Laptev sea shelf, experience extremely cold conditions in glacial climate, leading to negligible methane emissions. Exposed shelf areas in Indonesia, Africa, and South America, on the other 420 hand, emit significant amounts of methane. Thus lower sea level leads to larger emitting areas and thus higher emissions of methane.
3. The West African monsoon: During the time slices when the West African monsoon is stronger than at present, i.e. at 15 ka BP, 10 ka BP, and 5 ka BP, precipitation in the Sahel region is significantly enhanced in comparison to the PI state, leading to an increase in vegetation cover, productivity and biomass burning. As a result, methane emissions at 425 these times are stronger than at present, leading to a significant increase in (sub-)tropical CH 4 emissions, with all natural methane sources increased.
For methane emissions from wildfires, a further factor influencing the emissions is the human population density, as this strongly affects the fire probability in the SPITFIRE model employed in JSBACH (Fig. 10). The soil uptake of methane, on the other hand, is strongly dependent on the atmospheric concentration of methane (Fig. 10).

430
The changes in methane from LGM to the present are dominated by changes in tropical wetland emissions, with mid and high latitude wetland emissions being a significant but secondary factor, gaining in importance as the high latitudes become ice-free. In total the wetland emissions account for 93 − 96% of the net CH 4 flux, and all other methane sources are of minor importance.
Code and data availability. The primary data, that is the model code for MPI-ESM, is freely available to the scientific community and can be accessed with a license on the MPI-M model distribution website. In addition, secondary data and scripts that may be useful in  Table 2. Methane emissions for all time slices in T gCH 4 yr −1 . Shown are soil uptake, total wetland emissions, fire and termite emissions, net emissions, and wetland emissions from TRO and NXT.