Connect with us

Climate News Live

Methane, Monsoons, and Modulation of Millennial‐Scale Climate

Global Warming

Methane, Monsoons, and Modulation of Millennial‐Scale Climate

1 Introduction Well‐known variations in Earth’s orbital parameters have pronounced effects on its insolation (i.e., incoming solar radiation) and climate system. Also known as Milankovitch cycles, imprints of oscillatory changes in orbital precession, obliquity, and eccentricity are found in virtually all sufficiently resolved paleoclimatic records during and beyond the Cenozoic Era. Despite disagreement regarding causal…

Methane, Monsoons, and Modulation of Millennial‐Scale Climate

1 Introduction

Well‐known variations in Earth’s orbital parameters have pronounced effects on its insolation (i.e., incoming solar radiation) and climate system. Also known as Milankovitch cycles, imprints of oscillatory changes in orbital precession, obliquity, and eccentricity are found in virtually all sufficiently resolved paleoclimatic records during and beyond the Cenozoic Era. Despite disagreement regarding causal mechanisms (Abe‐Ouchi et al., 2013), it is well‐established that orbital forcing paced the glacial‐interglacial cycles of the late Pleistocene (Imbrie et al., 1992), with an average periodicity approximating 100,000 years (hereafter 100 kyr). Regional records of temperature and hydroclimate change also indicate shorter‐term variability (Caley et al., 2011) tied to precession (19‐ and 23‐kyr periodicities) and obliquity (41‐kyr periodicity). On even shorter timescales, arguably more pertinent to human and societal relevance (Castañeda et al., 2009), high‐resolution records reveal substantial millennial‐scale climate variability (1 to 10‐kyr periodicity) characterized by abrupt changes (e.g., Heinrich events Heinrich, 1988) and rapid swings (e.g., Dansgaard‐Oeschger cycles Dansgaard et al., 1984) superimposed on the more gradual, orbitally‐induced changes (Brook & Buizert, 2018). However, the relationship between millennial‐scale climate variability and orbital variations of insolation remains elusive.

Previously, orbitally‐paced changes in insolation have been hypothesized to modulate millennial‐scale activity in both atmospheric methane (Brook et al., 1996) and the Asian monsoons (Cheng et al., 2016). Studies have also postulated a strong coupling between monsoon intensity and methane at both precessional and millennial timescales, wherein increased monsoon rainfall is associated with increased methane production via wetland expansion in low‐latitude regions (Guo et al., 2011; Ruddiman & Raymo, 2003). Clarifying their relationship under changing background conditions is important for constraining extreme future scenarios as uncertainties in monsoon hydrology and resultant wetland changes feed back onto the future trajectory of methane (Brook & Buizert, 2018; Kirschke et al., 2013; O’Connor et al., 2010). However, the proposed mechanisms of how orbital forcing modulates millennial‐scale variability are inconsistent: Whereas Brook et al. (1996) find increased millennial‐scale methane activity with increased Northern Hemisphere summer insolation, Cheng et al. (2016) posit an antiphased modulation relationship where high Northern Hemisphere summer insolation corresponds to weakened millennial‐scale monsoon intervals.

Does orbital forcing modulate the magnitude of millennial‐scale variability? We define the term “modulate” following signal‐processing literature (Roder, 1931; Rial, 2000) wherein a higher‐frequency carrier signal (millennial fluctuations of the climate system) is modulated by a lower‐frequency input signal (orbital variations) to yield a resultant modulated signal (observed paleoclimate record). Thus, modulation implies a significant control of a longer timescale signal on shorter‐event amplitude or frequency. Here, we isolate the magnitude of millennial‐scale variability (MMV) of atmospheric methane, measured directly in ice cores, and that of composite Chinese speleothem δ18O, a purported proxy for East Asian monsoon intensity over the past 640 kyr, and quantitatively assess the extent to which primary orbital frequencies modulate the millennial band of climate variance. We stress that covariability over millennial timescales between atmospheric methane and speleothem δ18O (e.g., both records contain responses to Heinrich events; Marcott et al., 2014; Rhodes et al., 2015) does not preclude that the amplitude of these two parameters are differently modulated by external factors (Extier et al., 2018). Moreover, modulating frequencies need not comprise the maximal contribution of overall variability, although they need to be a significant contributor in order to be detected. In this work, we uncover a clear signature of precession and obliquity cycles in the MMV of methane, whereas we cannot find any imprint of insolation directly modulating the MMV of speleothem δ18O. We discuss the implications of this decoupling and outline plausible explanations for our finding, considering the impact of different interpretations of speleothem δ18O. Finally, we also find a long‐term trend in the MMV of speleothem δ18O coinciding approximately with the Mid‐Brunhes Event (∼430 ka), whereas no secular trend exists in the MMV of methane.

2 Methods

We focus on quantifying the MMV in the EPICA Dome Concordia (EDC) CH4 record (Loulergue et al., 2008) (Figure 1a) and the composite Chinese speleothem δ18O record (Cheng et al., 2016) (Figure 2a). The latter consists of several stalagmite δ18O records spliced together to form a composite from three caves in eastern China—Sanbao Cave (31.67°N, 110.43°E), Hulu Cave (32.5°N, 119.17°E), and Dongge Cave (25.28°N, 108.08°E). At first, we assume that the composite δ18O record purely reflects variations in East Asian monsoon rainfall intensity—the published interpretation (Cheng et al., 2016)—although there is ongoing debate regarding this interpretation (Beck et al., 2018; Zhang et al., 2018; Chiang et al., 2015; Clemens et al., 2018; Pausata et al., 2011), which we address in our discussion below. Both records are extremely well‐dated (U‐series dating for the speleothems, and ice‐flow modeling and gas chemistry combined with key stratigraphic datums for the ice core) and represent state‐of‐the‐art climate reconstructions in terms of their resolution over the past 640 kyr, spanning multiple realizations of Milankovitch cycles (see Cheng et al., 2016; Bazin et al., 2013, for details). We evaluate the MMV in these records by first filtering them to isolate their millennial components (Figures 1b and 2b). Next, we separate their positive and negative amplitude envelopes (Figures 1c and 2c) to independently assess the modulation of high‐value (e.g., Dansgaard‐Oeschger interstadials) versus low‐value (e.g., Dansgaard‐Oeschger stadials, Heinrich Events) millennial‐scale variability. We then calculate a continuous record of the MMV (Figures 1d and 2d), which reflects a combined metric of both positive‐ and negative‐amplitude variability, derived from the variance in the filtered records (Figures 1b and 2b). Finally, we utilize spectral analysis to determine periodicities where variance is concentrated (Figures 1e–1h and 2e–2h) and to identify the power of primary orbital frequencies in the MMV of atmospheric methane (MMVmethane) and that of composite Chinese speleothem δ18O (MMVChina). The methane and speleothem δ18O records are variable in their time‐step resolution as well as in their age‐model uncertainties; the ice‐core resolution becomes higher toward more recent values, whereas no observable trend is found in that of the speleothem record (although there is significant variability). We find that these effects have little impact on our overall results (Figure S1 and Figure S2). We also find that ice thinning or speleothem sampling resolution is unlikely to affect long‐term trends in our results (see Figure S2). Furthermore, our analysis is insensitive to the various ice core age models (Figure S1); for this work, we use the updated AICC2012 timescale (Bazin et al., 2013).


Orbital and millennial‐scale atmospheric methane variability over the past 640 kyr. (a) EPICA Dome C record of CH4, (b) its millennial‐scale variability calculated as the 10‐kyr high‐pass filtered record of the original time series, (c) high‐value and low‐value CH4 in the high‐pass filtered record, and (d) the magnitude of millennial‐scale variability (MMV) calculated as the centered rolling standard deviation of the high‐pass filtered record using 2‐kyr sliding windows (100‐year step). (e–h) Periodograms of corresponding time series using the Lomb‐Scargle methodology. Primary orbital frequencies are marked with dashed lines (19 and 23—precession; 41—obliquity; 101—eccentricity). Note different scaling for power spectral density. Strong peaks in the precessional band in the low‐CH4, high‐CH4, and MMV record indicate that insolation modulates the amplitude of millennial variations in atmospheric methane.


Orbital and millennial‐scale variability in the composite Chinese speleothem record over the past 640 kyr. (a) Composite Chinese speleothem δ18O, (b) its millennial‐scale variability calculated as the 10‐kyr high‐pass filtered record of the original time series, (c) negative‐value and positive‐value δ18O in the high‐pass filtered record, and (d) the magnitude of millennial‐scale variability calculated as the centered rolling standard deviation of the high‐pass filtered record using 2‐kyr sliding windows (100‐year step). (e–h) Periodograms of corresponding time series using the Lomb‐Scargle methodology. Primary orbital frequencies are marked with dashed lines (19 and 23—precession; 41—obliquity; 101—eccentricity). Note different scaling for power spectral density. Despite that the “raw” δ18O record contains strong concentrations of variance in the precessional band, the lack of peaks in the negative and positive events as well as in the MMV record suggest that precession does not modulate the envelope of millenial‐scale δ18O variations in the composite Chinese speleothem record.

The MMV in each record was calculated using the following procedure. The original time series was filtered using a Butterworth filter at a cutoff threshold of 10 kyr. Prior to filtering, we interpolated the original time series to a common step of 100 years (changing this value to 50 or 1,000 years did not impact our interpretations; see Figure S1). The MMV was calculated using a centered rolling standard deviation applied to the high‐pass filtered time series with a window of 2 kyr and a 100‐year step (i.e., essentially, a continuous window). The outcomes of our study remain unchanged if we reduced our cutoff threshold to 6 kyr or increased it to 12 kyr; we justify our choice of 10 so as to include the millennial band frequencies and exclude leakage from hemiprecession frequencies (Hagelberg et al., 1994).

We also compared the MMV in these records to that of sea‐level (Figure S3) from a recently compiled probabilistic stack (Spratt & Lisiecki, 2016) and also calculated the MMV in the following records: the EDC δ2H record from Antarctica (Jouzel et al., 2007) (Figures 4c and S4), the molybdenum X‐ray fluorescence (XRF) composite record from Cariaco Basin, Venezuela (Gibson & Peterson, 2014) (Figure S5), and the composite Borneo speleothem δ18O record (Carolin et al., 2016) from Malaysia (Figure S6). We normalized orbital eccentricity, tilt, and precession (ETP) parameters over the past 640 kyr for use in the wavelet coherence calculations (Figures 3b and 3d) to extract maximal phase and amplitude correlations with orbital forcing. Finally, all analyses evaluating the spectral characteristics of the time series in this study were performed using the Lomb‐Scargle periodogram, a method chosen so as to not generate spurious spectral characteristics (VanderPlas, 2017). Wavelet coherence between time series was performed in a Monte Carlo framework (n=1,000) using published codes (Grinsted et al., 2004) wherein each time series was linearly detrended. Python codes for our analyses are available at


Comparison of orbital forcing and MMV in atmospheric methane and composite Chinese speleothem δ18O. (a) MMVmethane (red) compared with precession (yellow), the predominant peak in the MMVmethane spectra (Figure 2h), and (b) wavelet coherence between MMVmethane and normalized eccentricity‐tilt‐precession (ETP) over the past 640 kyr. (c) MMVChina (violet) compared with obliquity (orange), the predominant peak in the MMVChina spectra (Figure 2h) and (d) wavelet coherence between MMVChina and ETP over the past 640 kyr. Primary orbital frequencies are marked with dashed lines (19 and 23—precession; 41—obliquity; 101—eccentricity), and lighter colors correspond to stronger coherence. The cone of influence, where edge effects could prevail, has been shaded. Black lines depict significance at the 10% level (number of Monte Carlo simulations = 1,000), and black arrows indicate phase, where those pointing to the right depict zero phase (upward arrows indicate that the phase of MMV leads ETP). Note that the axis for precession is inverted to show higher Northern Hemisphere summer insolation upward, and also note that the time series comparisons are with (a) precession and (c) obliquity, whereas the wavelet analyses use ETP to test sensitivity with all three aspects of primary orbital forcing (where precession is multiplied by −1 as low precession corresponds to high ETP and high Northern Hemisphere insolation). Note that coherence between precession and MMVmethane is persistent and in phase, whereas none of the orbital frequencies are in phase or persistently coherent with MMVChina.


Comparison of MMV in composite Chinese speleothem δ18O with atmospheric methane and Antarctic ice core δ2H. (a) MMVmethane (red) compared with MMVChina (violet) compared with obliquity (orange) and (b) wavelet coherence between the two records over the past 640 kyr. (c) MMVChina (violet) compared with MMVAntarctica (green) and (d) wavelet coherence between the two records over the past 640 kyr. Primary orbital frequencies are marked with dashed lines (19 and 23—precession; 41—obliquity; 101—eccentricity), and lighter colors correspond to stronger coherence. The cone of influence, where edge effects could prevail, has been shaded. Black lines depict significance at the 10% level (number of Monte Carlo simulations = 1,000), and black arrows indicate phase, where those pointing to the right depict zero phase (upward arrows indicate that the phase of MMV leads ETP). Note that MMVChina and MMVAntarctica are persistently coherent and in phase over the length of their records, whereas MMVChina and MMVmethane show little direct correspondence.

3 Results

Spectral analysis reveals the dominant influence of orbital forcing in atmospheric methane variability over the past 640 kyr. In the original (“raw”) methane record, the variance is mainly concentrated in the obliquity and eccentricity bands (Chappellaz et al., 1990; Loulergue et al., 2008), to nearly the same magnitude (Figure 1e). Peaks in the precessional bands are also present, although their magnitudes are comparatively small. Power in the millennial band becomes prominent in the filtered record, which, by construction, does not contain variance in periodicities beyond 10 kyr (Figure 1f; note change of scale in axis). Spectra of both high‐value and low‐value excursions contain separate and distinct peaks in the millennial bands (Figure 1g), potentially pointing to separate and distinct controls on uptake and release of CH4 during warm and cold periods (Siddall et al., 2010). However, over modulating timescales (>10 kyr), we find that the spectral peaks for both high‐value and low‐value CH4 are identical. These peaks are also very similar to those observed in the MMVmethane record (Figure 1h). Most notably, we find strong spectral power in both precessional bands as well as power in the obliquity band for the high‐value and low‐value records as well as in the MMV of methane (Figures 1g and 1h).

On closer inspection, we find that precession minima (i.e., Northern Hemisphere summer insolation maxima) correspond to large increases in MMVmethane (Figure 3a). Many of these increases are associated with the release and subsequent drawdown of atmospheric methane coinciding with precession minima, nearing the culmination of deglaciation (Figures 3a and S3; e.g., over the Bølling‐Allerød and Younger Dryas sequence). Conversely, precessional maxima (i.e., Northern Hemisphere summer insolation minima) correspond to reduced MMVmethane (Figure 3a). This relationship is insensitive to glacial or interglacial background state (Figure S3). The amplitude of these swings in MMV often follow the amplitude of precession changes; for example, a lull in MMVmethane is observed at ∼400 ka, when variability in precession is subdued. These changes appear to be threshold‐limited (e.g., low values of insolation limit MMVmethane to near‐zero values) and preclude a straightforward investigation of the sensitivity between MMVmethane and precession. Nevertheless, wavelet analysis reveals strong coherence and near‐zero phase between MMVmethane and orbital forcing at the 19 and 23 kyr bands over the last two glacial cycles, with varying but persistent power over the past 640 kyr (Figure 3b). Comparatively weak coherence and out‐of‐phase relationships are found at the 41 and 101 kyr bands, respectively (Figure 3b), despite the peaks in the periodograms (Figures 1g and 1h), suggesting that MMVmethane is not directly modulated by eccentricity and might be faintly influenced by obliquity. Thus, we conclude that precession directly and robustly modulates millennial‐scale variability in atmospheric methane, supporting the assertion of higher insolation facilitating increased variability (Brook et al., 1996).

In stark contrast, we find no trace of precessional forcing in the modulation of millennial‐scale speleothem δ18O variability (Figure 2). Variance in the original (“raw”) δ18O record is dominated by precession (Figures 2a and 2e), yet power is absent in the precessional bands of both the high‐value and low‐value records (Figures 2c and 2g), as well as in the MMV record (Figures 2d and 2h). Instead, we observe power at a spectral peak near the obliquity band (∼43 kyr). A time series comparison of MMVChina and obliquity demonstrates that (1) the correspondence between them is not consistent over the past 640 kyr (Figure 3c) and (2) MMVChina peaks lead obliquity from 350–640 ka (Figure 3d). Hence, MMV in composite Chinese δ18O is not modulated by precession, and moreover, its relationship with obliquity is complex, as it leads obliquity‐driven insolation changes. Our results differ from Cheng et al. (2016), wherein the suborbital (i.e., millennial) component was calculated using the residual between normalized δ18O and normalized insolation, versus the discretized filters employed herein. We conclude that insolation does not modulate MMVChina and that Milankovitch cycles are not the “pacemaker of millennial‐scale events” in the suborbital component of composite Chinese speleothem δ18O.

Our analysis shows that the modulating agents of millennial‐scale variability in atmospheric methane and in Chinese stalagmite δ18O are disparate. Orbitally paced insolation, via precession, modulates MMVmethane but does not modulate MMVChina; the MMV time series are not correlated and also show limited coherence in the primary orbital bands (Figures 4a and 4b). On the other hand, we find a striking correspondence between MMVChina and the MMV in Antarctic δ2H (MMVAntarctica—see Figures 4 and S4 for details). Significant and in‐phase coherence (Figures 4c and 4d) between these records indicates a persistent and robust coupling over the past 640 kyr, with trends toward higher‐magnitude millennial‐scale variability in both records (particularly pronounced in MMVAntarctica; Figure S4d) starting around ∼430 ka (Figure 4c). Ultimately, the MMV in composite Chinese speleothem δ18O and Antarctic ice core δ2H share more commonalities in their modes of temporal variability than with that of atmospheric methane—where the latter is modulated by precession.

4 Discussion

What does the decoupling of orbital modulation of the MMV in methane and monsoon intensity imply? Studies over the past three decades have proposed that low‐latitude precipitation, via its influence on tropical wetlands, is a major driver of long‐term variability in atmospheric methane over glacial cycles (Chappellaz et al., 1990). Although this interpretation has been challenged (Crowley, 1991; Schmidt et al., 2004), it has also been used to tune ice‐core age‐models to precession cycles (Ruddiman & Raymo, 2003). Shifts in the intertropical convergence zone (ITCZ) and monsoon rainfall variability have been invoked as causes for orbital and millennial‐scale variability in methane (Brook & Buizert, 2018; Guo et al., 2011), with a secondary role for boreal sources (Chappellaz et al., 1997) and a minimal role for “geological” emissions including marine clatharates, submarine seeps, mud volcanoes, and permafrost thaw (Bock et al., 2017; Petrenko et al., 2017). Yet large uncertainties remain in pinpointing the relative contributions of various sources to past methane changes, and the subsequent significance for future ramifications remain hotly debated (Baumgartner et al., 2012; Bock et al., 2017; Etiope & Schwietzke, 2019; Kirschke et al., 2013; Lamarche‐Gagnon et al., 2019; Schuur et al., 2015; Petrenko et al., 2017).

Our results point to midlatitude and high‐latitude sources as potentially important drivers of methane variability during abrupt climate change events, given that precession modulates MMVmethane but not MMVChina. As we discuss below, this inference relies on the assumption that composite Chinese speleothem δ18O reflects some component of millennial‐scale variability in East Asian summer monsoon rainfall intensity or upstream changes in tropical vapor (Cheng et al., 2016). If taken at face value, the in‐phase behavior of MMVmethane with Northern Hemisphere summer insolation (Figure 3b) and decoupling with Northern Hemisphere subtropical monsoon rainfall (Figure 4a) suggests that boreal sources of methane are sensitive to precession and obliquity and contribute to the modulation of atmospheric methane during millennial‐scale events.

Several mechanisms involving boreal sources as major contributors to methane variations have been proposed, including variability in midlatitude wetlands and peatlands from changes in sea level and hydroclimate (Baumgartner et al., 2012; Bock et al., 2017; Macdonald et al., 2006; Ridgwell et al., 2012), high‐latitude permafrost thaw (O’Connor et al., 2010; Zimov et al., 2006), methane trapped under ice sheets (Lamarche‐Gagnon et al., 2019), and emissions from temperature‐sensitive thermokarst lakes (Lewkowicz & Way, 2019; Walter et al., 2007; Walter Anthony et al., 2018). We contend that precession‐triggered changes in local summer temperatures and associated hydrological anomalies in middle to high northern latitudes could generate shifts in methane sources through a combination of these mechanisms and thereby augment millennial‐scale variability in atmospheric methane, explaining the modulation via precession‐driven insolation. MMVmethane would be amplified during periods of boreal summer insolation maxima due to increased variability in available wetland, peatland, and permafrost sources as well as areas of ice melt resulting from more extreme summer temperatures and hydroclimate seasonality (Lewkowicz & Way, 2019). For example, Petrenko et al. (2017) used 14CH4 measurements to investigate the ∼250 ppb rise in atmospheric methane after the Younger Dryas and estimated that up to ∼20% might be derived from “geological sources,” including permafrost thaw, although it is worth noting that they do not discriminate between tropical or extratropical wetlands in the remnant contribution. Regardless, we hypothesize that this ∼20% rise (as a minimum) would not be as high during previous millennial‐scale events which occurred during boreal summer insolation minima—a scenario testable with additional data or by using biogeochemistry‐resolving model simulations of abrupt climate change forced under different background orbital conditions. However, approaches focusing on one or few events ought to be carefully scrutinized as not all millennial‐scale climate changes are equally created (Dansgaard et al., 1984; Heinrich, 1988). Over the late Pleistocene, we posit that millennial‐scale variability of atmospheric methane is more subdued during periods of precession maxima presumably because these midlatitude to high‐latitude sources would be less active owing to milder summer temperatures (Zheng et al., 2020).

Additionally, tropical Southern Hemisphere sources could play an important role in buffering the amplitude of millennial‐scale methane variability and contribute toward the modulation by insolation. During Northern Hemisphere stadial events, the ITCZ is thought to occupy a more southward position alongside an intensification of the austral monsoons (Kanner et al., 2012). These processes would potentially be stronger during periods of austral summer maxima (i.e., precession maxima) (Wang et al., 2006) and facilitate more wetland‐based sources of methane which would serve to offset coeval CH4 decreases from the Northern Hemisphere (Rhodes et al., 2015). The resultant damping of overall variability in atmospheric methane during these periods is consistent with MMVmethane being in phase with Northern Hemisphere summer insolation. Well‐dated hydroclimate records from the Southern Hemisphere tropics spanning multiple Milankovitch cycles could confirm this hypothesis if their MMV was coherent and in phase with local (austral) insolation at the precession and obliquity bands. Curiously, however, we note that MMVChina does not exhibit this behavior (Figures 2 and 3). If confirmed, spectral signatures of precession in the MMV of yet to be generated Southern Hemisphere records will necessitate the following outcome: Either composite Chinese speleothem δ18O is not representative of global paleomonsoon process; that is, the purported “give and take” hydrological balance of rainfall (invoked alongside ITCZ shifts) between the hemispheres over millennial‐scale events is not symmetric (Cheng et al., 2016; Caley et al., 2011; Wang et al., 2007), or speleothem δ18O is altogether decoupled from regional Asian monsoon rainfall over millennial timescales (Zhang et al., 2018).

Numerous studies have highlighted complications in strictly interpreting Chinese speleothem δ18O as monsoon rainfall amount. However, even if speleothem δ18O were not a representative metric of the global paleomonsoons (Cheng et al., 2016; Caley et al., 2011) and instead reflected changes in the position of westerlies, or upstream vapor changes in the Indian monsoon or tropical domains over millennial or longer timescales (Chiang et al., 2015; Hu et al., 2019; Pausata et al., 2011), the remarkable similarity between MMVChina and MMVAntarctica indicates that disparate and distant regional ocean‐atmosphere features of the climate system are modulated in the same manner (Figures 4c and 4d). The lack of precession in these MMV records suggests that endogenous processes modulate the magnitude of millennial‐scale climate variability monitored by these records.

Does our analysis refute a tropical driver for the precession forcing of MMVmethane? If MMVChina does not capture East Asian monsoon variability or reflect concurrent changes in tropical rainfall, and if the mechanism behind the observed precession modulation of MMVmethane arose from the tropics, then the MMV calculated in records of tropical precipitation from the monsoon and ITCZ domain ought to contain precessional peaks. Such long, well‐resolved, and absolutely dated records of tropical rainfall do not yet exist to reject this hypothesis. As a preliminary investigation, we found that the MMV in a record of ITCZ changes from the Cariaco Basin (Gibson & Peterson, 2014), off Venezuela (shorter in length compared to the composite Chinese speleothem record), does not contain power in the precessional band as a modulator (Figure S5) potentially withdrawing northern South America as an area from where the mechanism for modulation arises. Speleothem δ18O from Borneo, although spanning only ∼150 kyr, also does not contain such a peak (Figure S6). We stress, however, that long, late Pleistocene reconstructions of hydroclimate for the Congo Basin and much of South America, including the Amazon basin—key regions for wetlands (Kaplan, 2002)—do not yet exist and could contain precessional signatures. Yet due to the above inferences as well as the correspondence and lack of precession in MMVChina and MMVAntarctica, we speculate that boreal midlatitude and high‐latitude sources drive the precessional modulation of MMVmethane. Instead, the similarity between the MMV in composite Chinese speleothem δ18O and Antarctic ice core δ2H suggests that MMVChina is modulated by high‐latitude yet internal sources of climate variability.

The bipolar seesaw paradigm links millennial‐scale climate change between Antarctica and Greenland via the Atlantic Meridional Overturning Circulation (AMOC) (Barker et al., 2011; Brook & Buizert, 2018; Siddall et al., 2010). Barker and colleagues generated a synthetic record of millennial‐scale Greenland temperature variability (GLT_syn), beyond the past glacial cycle, derived from Antarctic ice core δ2H, assuming stationarity in this paradigm (Barker et al., 2011). Accordingly, MMVAntarctica is virtually identical to the MMV in the synthetic series (Figure S7), which leads to the need for reconciliation: If MMVmethane is modulated by high‐latitude precession via changes in local temperature and hydroclimate, why is this not reflected in
urn:x-wiley:grl:media:grl60531:grl60538-math-0001, by extension from MMVAntarctica? One possibility is that the bipolar seesaw is not stationary (Siddall et al., 2010) over the past 640 kyr, implying asymmetry in hemispheric responses on millennial timescales, although more highly resolved records are needed to explore this hypothesis. However, another reconciling explanation is offered by the delayed response of ice sheets to external forcing. The waxing and waning of ice sheets offer more inertia to insolation change (Abe‐Ouchi et al., 2013) and can set their own millennial‐scale amplitude which varies from event to event (Lynch‐Stieglitz et al., 2014), either triggered by AMOC changes (Marcott et al., 2011) or otherwise (Li et al., 2010; Petersen et al., 2013). Thus, whereas boreal sources of methane can respond rapidly to direct insolation forcing (Lewkowicz & Way, 2019) and are modulated by precession, internal ice‐ocean‐atmosphere dynamics modulate abrupt climate change and subsequent downstream impacts, even if both processes can start rapidly (Barker et al., 2011; Li et al., 2010). Such teleconnections internal to the climate system can also explain why MMVChina and MMVAntarctica share many common traits. Considering that these records are coherent and in‐phase at the obliquity band (Figure 4d), despite not being in‐phase with obliquity‐forced changes in insolation (Figure 3d) nor having power at the precessional bands (Figures 2h and S4h), it is strongly suggested that internal interactions of the ice‐ocean‐atmosphere system set the tempo and magnitude of their millennial‐scale climate variability.

Finally, we note the occurrence of a trend toward higher‐amplitude millennial‐scale variability in the latter half of the MMV in the Chinese and Antarctic records (Figure 4c). This trend of increasing MMV in both records is independent of the time‐step of the raw data sets (see Figure S2). The onset of this trend coincides with the Mid‐Brunhes Event (∼430 ka), when ice sheets increased in size and the 100‐kyr cycle became more prominent (Wang et al., 2003). Changes in the carbon reservoir as well as the effect of insolation on the Southern Hemisphere have been invoked to explain this event (Yin & Berger, 2010; Wang et al., 2003). According to our analysis, stronger glacial‐interglacial cycles coincide with stronger‐magnitude millennial‐scale climate variability. As insolation does not trend over the last 400 kyr, this observation provides an independent line of evidence that insolation does not modulate MMVChina nor MMVAntarctica, and instead, processes involving the ice‐ocean‐atmosphere system may modulate them. Curiously, we find that there is no trend in the millennial‐scale activity of methane before or after this event, which reaffirms our hypothesis that MMVmethane is modulated by changes in Northern Hemisphere summer insolation linked to precession cycles.

5 Conclusions

We provide a new framework to isolate and study millennial‐scale variability and address its modulation in well‐dated late Pleistocene records. Our findings denote that precession directly modulates the timing and magnitude of atmospheric methane variations over millennial timescales but not that of East Asian monsoon variability, as recorded in speleothem δ18O (Figures 1 and 2). If the composite Chinese speleothem δ18O record indeed reflects either East Asian monsoon intensity or upstream hydroclimate changes from the tropics, then the lack of precession in MMVChina suggests that future work should focus on midlatitude or high‐latitude methane sources as the origin for insolation‐forced modulation of MMVmethane. On the other hand, if the subtropical composite δ18O record is decoupled from variability in tropical sources of methane, then we cannot rule out the possible contribution of the tropics and must await similarly detailed reconstructions from these regions. Finally, the strong correspondence between the MMV in Antarctic EDC ice‐core δ2H and the composite Chinese speleothem δ18O record (Figure 4) indicates that internal feedbacks involving the cryosphere, and not orbital forcing, modulate the amplitude and frequency of millennial‐scale climate variability.


We are grateful to Warren Prell for commenting on a previous version of this manuscript and to Baylor Fox‐Kemper who provided input on the statistical analysis. We thank Yoshimi Kubota, Joseph Orchardo, and Nick Scroxton for fruitful discussion. We also thank two anonymous reviewers and Associate Editor Sarah Feakins of Geophysical Research Letters for their suggestions to improve this manuscript. K. T. acknowledges the Brown University Presidential Fellowship and the Department of Geosciences at the University of Arizona. This research was supported by NSF grants AGS‐1703123 to K. T. and J. W. P. and OCE‐1634774 to S. C. C. The MMV data sets calculated by the authors are attached as a supporting information spreadsheet; all data used herein have been previously peer‐reviewed and available: sea level from Spratt and Lisiecki (2016), the molybdenum XRF composite record from Gibson and Peterson (2014), the composite Borneo speleothem δ18O record from Carolin et al. (2016), the composite Chinese speleothem δ18O record from Cheng et al. (2016), atmospheric methane from Loulergue et al. (2008), and Antarctic δ2H from Jouzel et al. (2007).

    Subscribe to the newsletter news

    We hate SPAM and promise to keep your email address safe

    Click to comment

    Leave a Reply

    Your email address will not be published. Required fields are marked *

    Top Stories

    To Top