Climate evolution across the Mid-Brunhes Transition

The Mid-Brunhes Transition (MBT) began ∼ 430 ka with an increase in the amplitude of the 100 kyr climate cycles of the past 800 000 years. The MBT has been identified in ice-core records, which indicate interglaciations became warmer with higher atmospheric CO₂ levels after the MBT, and benthic oxygen isotope (δ¹⁸O) records, which suggest that post-MBT interglaciations had higher sea levels and warmer temperatures than pre-MBT interglaciations. It remains unclear, however, whether the MBT was a globally synchronous phenomenon that included other components of the climate system. Here, we further characterize changes in the climate system across the MBT through statistical analyses of ice-core and δ¹⁸O records as well as sea-surface temperature, benthic carbon isotope, and dust accumulation records. Our results demonstrate that the MBT was a global event with a significant increase in climate variance in most components of the climate system assessed here. However, our results indicate that the onset of high-amplitude variability in temperature, atmospheric CO₂, and sea level at ∼430 ka was preceded by changes in the carbon cycle, ice sheets, and monsoon strength during Marine Isotope Stage (MIS) 14 and MIS 13.

1 Introduction

The last 800 kyr of the Pleistocene epoch is characterized by the emergence of dominant ∼100 kyr glacial–interglacial climate cycles (Pisias and Moore, 1981; Imbrie et al., 1993; Raymo et al., 1997; Clark et al., 2006). These climate cycles typically have long glacial periods punctuated by short interglaciations. Since ∼430 ka (i.e., starting with Marine Isotope Stage (MIS) 11), interglaciations have experienced warmer temperatures (Jouzel et al., 2007) and higher concentrations of atmospheric CO₂ (Lüthi et al., 2008) relative to earlier interglaciations of the last 800 kyr (Fig. 1). The transition to higher-amplitude interglaciations has also been recognized in deep-sea records of δ¹⁸O measured in benthic foraminifera (Lisiecki and Raymo, 2005) that identify lesser ice volume and/or warmer deep-ocean temperatures (Fig. 1).

Jansen et al. (1986) originally described this change in amplitude of interglaciations as a singular Mid-Brunhes Event, but Yin (2013) argued that it is more appropriately considered as a transition between two distinct climate states, thus referring to it as the Mid-Brunhes Transition (MBT). The change from low-amplitude to high-amplitude 100 kyr variability at ∼430 ka occurs during an interval of reduced eccentricity and corresponding precession (Fig. 1), but similar orbital forcing occurred at times before and after the onset of the MBT with no comparable response, suggesting that the MBT was an unforced change internal to the climate system. Mechanisms proposed for the MBT include a latitudinal shift in the position of the Southern Hemisphere westerlies that increased upwelling of respired carbon in the post-MBT Southern Ocean (Kemp et al., 2010) and a change in Antarctic Bottom Water (AABW) formation through insolation-induced feedbacks on sea ice and surface water density (Yin, 2013). However, several questions remain. (1) How and when was the MBT expressed in other components of the climate system? (2) Was the MBT a global or regional transition? (3) Did components expressing a transition change synchronously? Here, we address these questions by providing a statistical characterization of changes occurring over the last 800 kyr as recorded by a variety of paleoclimatic proxies with broad spatial coverage.

Figure 1 Orbital forcing and climate records for the last 800 kyr. (a, b, c) Precession, obliquity, and eccentricity (Laskar et al., 2004). (d) Deuterium-derived temperature from the European Project for Ice Coring in Antarctica (EPICA) Dome C ice core in Antarctica (Jouzel et al., 2007). (e) Atmospheric CO₂ from EPICA Dome C (EPICA-community-members, 2004; Lüthi et al., 2008). (f) Atmospheric CH₄ from EPICA Dome C (EPICA-community-members, 2004). (g) Global benthic oxygen isotope stack (Lisiecki and Raymo, 2005).

Figure 2 Site locations. Map indicating the locations of the cores used in this research and modern sea-surface temperature values. Each symbol represents a different proxy record. Diamonds indicate sea-surface temperatures. Circles indicate benthic ∂¹³C. Squares indicate dust.

2 Methods

2.1 Data collection

We compiled all available published records of sea-surface temperature (SST), benthic marine carbon isotopes ratios (δ¹³C), and dust accumulation (dust) that met our selection criteria and closely represented a global distribution as attainable (Fig. 2). Each data set has an average temporal resolution of <5 kyr, does not include any large age gaps, and spans much or all of the entire time period of consideration to limit biasing of the younger parts of the record. Lisiecki (2014) placed all of the δ¹³C records on the LR04 age model. Published SST records that were not on the LR04 age model were placed on it in one of two ways. If the original data had depth and benthic δ¹⁸O data, the SST record was placed on LR04 using the ager script in MATLAB as part of the ARAND software package (Howell et al., 2006). When only benthic δ¹⁸O records were available, the SST records were placed on LR04 by selecting corresponding tie points in the δ¹⁸O data series using the AnalySeries version 2.0 software (Paillard et al., 1996). Because some dust records could not be placed on the LR04 age model, certain statistical analyses of them (e.g., phase–lag relationships) are likely not robust, but the overall variance in them is preserved. Each record was then interpolated to a time step (Δt) of 2 kyr. With each record having an average resolution <5 kyr, this Δt allows for the preservation of higher-frequency variability while limiting the number of interpolated data points.

We used empirical orthogonal function (EOF) analysis to characterize the dominant modes of variability and robustly demonstrate global and regional signals of the SST, δ¹³C, and dust records. We then used spectral analyses of each resulting principal component (PC) to characterize their periodicity, phase, and amplitude.

2.2 Sea-surface temperatures

We used 11 SST records that span the entire 800 kyr time period and four additional records that span 8–758 ka. Inclusion of these four shorter records does not change our conclusions. The SST records cover the Pacific (n=9), Atlantic (n=5), and Indian (n=1) oceans (Fig. 2, Table 1). We note that Shakun et al. (2015) reconstructed a global SST stack for the last 800 kyr using 49 records, but only seven of these spanned the entire 800 kyr. Comparison of our SST PC1 based on 15 records to the Shakun SST stack shows excellent agreement (Fig. S1 in the Supplement).

2.3 Carbon isotopes (δ¹³C)

We analyzed the global set of δ¹³C records compiled by Lisiecki (2014) (n=26; Fig. 2), and separately analyzed the records in the Atlantic (n=14) and the Pacific (n=4) basins, thus distinguishing between the dominant water masses within each basin and removing the muting effect of the more negative Pacific values on the more positive Atlantic. Similar to SSTs, Lisiecki (2014) reconstructed a global δ¹³C stack for the last 3 Myr using 46 records, but only 18 of these spanned the last 800 kyr. Comparison of our δ¹³C PC1 to the Lisiecki δ¹³C stack shows excellent agreement (Fig. S2).

We then looked at regional and depth stacks of the δ¹³C records in the Atlantic basin to characterize changes in the dominant water masses on orbital timescales. Regional stacks were broken into North Atlantic (>20^∘ N; n=4), equatorial Atlantic (20^∘ S to 20^∘ N; n=14), and South Atlantic (>20^∘ S; n=8). We also created stacks for the deep North Atlantic (depth >2000 m; n=4) and intermediate North Atlantic (depth <2000 m; n=3). All included records were averaged to create the stack and each stacked record was interpolated to a 2 kyr time step. Stacking improves the signal-to-noise ratio of the δ¹³C records, making regional stacks useful in identifying circulation changes and comparing circulation responses with other climate records (Lisiecki, 2014).

Table 1 Data compilation. All data sets are used in these analyses with associated locations, proxy type, references, and digital object identifier when available.

2.4 Dust

We analyzed eight proxy records of dust that span the entire 800 kyr time period, and then separated them by hemisphere (Northern Hemisphere had three; Southern Hemisphere had five) to characterize hemispheric differences (Fig. 2). The various proxies for dust include Fe mass accumulation rates, weight percent of terrigenous material and Fe, flux of lithogenic grains, and grain size analysis. We standardized each record before analysis to account for these various proxy types and their differing range in values, thus allowing for comparison of their relative amplitudes of variation.

2.5 Empirical orthogonal function (EOF) analysis

We used EOF analysis to objectively characterize the climate variability recorded by the proxies across the MBT. Analyses of covariance between the data were conducted using the EOF script as part of the ARAND software package (Howell et al., 2006). The results provide both the dominant variability as a time series (principal component) and a spatial distribution of variance contribution (factor loadings). The records for SST and δ¹³C were kept in their original values of degrees and per mil, respectively, to preserve the original variance. Dust records were standardized to a mean value of zero and unit variance so that each record provided equal weight to the EOF. Statistical significance of all EOFs was determined through segmented linear regression analysis. All resulting break points occur on or after the second EOF and are thus considered significant.

2.6 Spectral analysis

We used the Blackman–Tukey technique in the ARAND software package for spectral analysis of each PC (Howell et al., 2006). Analyses were conducted using all data points within the time interval of interest, boxcar windowing of the input data, and hamming spectral filter. Multiple tests were conducted for the 8–800, 450–800, and 8–350 ka time slices. These intervals characterize the dominant frequency of variability over the entire 800 kyr record, and for the pre- and post-MBT intervals, respectively. The removal of the 350–450 ka interval limited the influence of MIS 11, MIS 12, and Termination V (T5), as these were shown to potentially bias the spectral power. Furthermore, these selected intervals result in time series of equal length to limit biasing of longer records. Additional tests were conducted using wavelet analyses that characterize the change in spectral power as a time series. Complementary spectral analyses were conducted on CO₂ and CH₄ records from the European Project for Ice Coring in Antarctica (EPICA) Dome C ice core (EPICA-community-members, 2004; Jouzel et al., 2007), and benthic δ¹⁸O using the LR04 stack (Lisiecki and Raymo, 2005). Cross-spectral analyses were conducted for the PCs against mean insolation values to determine phase and coherency of each PC. Mean insolation values were calculated for each of the dominant periodicities (eccentricity, obliquity, and precession) with the data derived from AnalySeries (Laskar et al., 2004; Paillard et al., 1996).

2.7 Variance tests

We used f tests to test for variance changes across the MBT for each principal component from the EOF analysis as well as for CO₂, CH₄, and the LR04 δ¹⁸O records. Analyses were conducted in MATLAB using the vartest2 script. This approach assumes the null hypothesis that the pre- and post-MBT distributions of the time series of each climate component have the same normally distributed variance. If the resulting variance values reject this hypothesis of no statistical difference, then the pre- and post-MBT time series are determined to have undergone a significant change in variance across the MBT. We interpret the change in variance to reflect a change in the amplitude of each climate signal.

3 Results

3.1 CO₂, CH₄, and benthic δ¹⁸O

Time series of the greenhouse gases CO₂ and CH₄ and of the LR04 stack of benthic δ¹⁸O suggest an increase in their interglacial values across the MBT (Fig. 1). Spectral analyses of the LR04 stack and atmospheric CO₂ indicate a small post-MBT increase in the 100 kyr band, whereas results for CH₄ indicate a decrease (Fig. S3). All three records show an increase in the precessional band (19–23 kyr). Variance tests suggest that δ¹⁸O and CO₂ have a statistically significant increase in variance across the MBT, while CH₄ variance decreases (Table S1).

3.2 Sea-surface temperatures

EOF analysis of global SSTs over the last 758 kyr identifies two statistically significant principal components (Fig. 3a). The first and second principal components (PC1 and PC2, respectively) account for 69 % of the total variance, with PC1 explaining 49 % alone. While some degree of regional variability in each record exists, factor loadings indicate that each record positively contributed to PC1 with a larger contribution coming from high-latitude records. Thus, PC1 is representative of a global SST signal. SST PC1 demonstrates a stepwise increase in variance starting at 436 ka, with an increase of interglacial temperatures, while showing no significant change in the lower limit glacial values, which is one of the defining characteristics of the MBT. The highest spectral density is in the 100 kyr frequency band throughout the entire time period (Fig. S3d). Wavelet analysis (Fig. 4a) shows a significant increase in the 100 kyr frequency band at 580 ka that reaches its maximum spectral power during MIS 11 and persists throughout most of the remaining interval, albeit with decreasing intensity after ∼250 ka. Variance f tests reveal a significant increase in amplitude from the pre- to post-MBT SSTs (Table S1). These results thus confirm that there was a stepwise global transition of SSTs from lower- to higher-amplitude interglaciations as previously inferred from individual records.

Figure 3 Principal components. Plots of the first (PC1; blue) and second (PC2; red) principal components from our EOF analysis of each climate variable. Percent variance is explained by each PC represented by the numbers with the corresponding color. (a) Sea-surface temperatures. (b) Dust records. (c) Global ∂¹³C. (d) ∂¹³C of the Atlantic. (e) ∂¹³C of the Pacific.

Figure 4 Wavelet analysis. Wavelets of four of the first principal components. (a) Sea-surface temperatures. (b) Dust records. (c) Global ∂¹³C. (d) ∂¹³C of the Atlantic. Red colors represent higher spectral power. Blue colors represent lower spectral power. Statistical significance is highlighted by the thin black line. Milankovitch periods are highlighted by the dashed horizontal lines.

Variance calculations on proxies of bottom water temperature (Elderfield et al., 2012) and on the Antarctic EPICA ice-core deuterium record (EPICA-community-members, 2004), a measure of Antarctic atmospheric temperature, also indicate statistically significant increases in variance across the MBT (Table S1). In both proxies, the time series indicate an increase of interglacial temperature values while showing no significant change to the lower limit glacial values, similar to PC1 of SSTs (Fig. 5).

3.3 Dust

The EOF analysis of the global dust records identifies two statistically significant principal components, with PC1 representing 56 % of the total variance and PC2 15 % (Fig. 3b). All records but the one from the Chinese Loess Plateau (CLP) reflect increased dust accumulation due to increased aridity and/or wind strength during glaciations, whereas higher dust accumulation in the CLP record reflects increased summer Asian monsoon strength, which is an interglacial signal (Sun and An, 2005). Accordingly, factor loadings for the dust records are all positive for PC1 except for the CLP.

Figure 5 Temperature records. (a) Deuterium-based temperature record from EPICA Dome C in Antarctica (light yellow; Jouzel et al., 2007). The darker yellow line is a 15-point moving average. (b) The first principal component of our sea-surface temperature analysis (red). (c) Bottom water temperature derived from Mg ∕ Ca measurements at ODP 1123 (light blue; Elderfield et al., 2012). The dark blue line is a 15-point moving average.

In contrast to the change in variance seen in temperature, CO₂, and CH₄ during MIS 11, variance tests of the dust PC1 suggest a stepwise increase in variance during MIS 12, with subsequent glaciations having higher amplitudes (Table S1). Separating the records by hemisphere shows that the increase in glacial amplitude starting at MIS 12 occurs in the southern PC1 but not in the northern PC1 (Fig. 6). Similarly, the signal during MIS 14 present in the global PC1 is absent in the northern PC1, suggesting that the northern control on dust accumulation was skipped during that glacial.

Spectral analysis of the global PC1 indicates dominant power in the 100 kyr frequency band that increases in spectral power across the MBT (Fig. S3b). Furthermore, wavelet analysis of PC1 demonstrates an increase in the spectral power of the 100 kyr band at ∼600 ka with its highest power during MIS 11 (Fig. 4b), similar to the SST PC1. The 100 kyr frequency remains statistically significant throughout the 100–600 ka interval.

3.4 δ¹³C

The first principal component of the global δ¹³C ($\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{G}}$; PC1) explains 58 % of the total variance (Fig. 3c). EOF analysis of δ¹³C records from the Atlantic basin ($\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{ATL}})$ yields two statistically significant PCs, with PC1 and PC2 explaining 58 % and 13 % of the total variance, respectively (Fig. 3d). EOF analysis of δ¹³C records from the Pacific ($\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{PAC}})$ yields one statistically significant principal component (PC1 is 81 % total variance) (Fig. 4e).

Figure 6 Dust principal components. The first principal components of our dust analysis for the global (yellow), north (red), and south (blue) records. Vertical gray boxes highlight specific glacial (dark gray) and interglacial (light gray) periods. The numbers indicate the associated MIS of each box.

Both the global and Atlantic PC1 exhibit a strong 100 kyr frequency that is persistent from 680 to 180 ka (Fig. 4c, d). Unlike SST and dust, however, $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{G}}$ and $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{ATL}}$ demonstrate a stronger 100 kyr power prior to MIS 11 with its highest power throughout MIS 13 and 12 (510–460 ka). Spectral analysis shows a decrease in power of the 100 kyr frequency band from pre- to post-MBT (Fig. S3f, g). Variance tests show that the pre- and post-MBT intervals for $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{G}}$ and $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{ATL}}$ are statistically different with higher variance during the pre-MBT (Table S1). Spectral analyses and variance tests of $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{PAC}}$ PC1 are similar to $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{G}}$ and $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{ATL}}$ PC1s. The only difference between the three PC1s is that there is less variance recorded in $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{PAC}}$ (Fig. 3e). We interpret this muted signal to be a result of three factors: the large size of the Pacific relative to the Atlantic, less mixing between water mass end members such as the positive NADW and more negative AABW, and ocean circulation aging the carbon isotopes over time leading to more homogenized water masses in the Pacific.

Factor loadings for $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{ATL}}$ PC1 are all positive, suggesting that the time series is representative of the entire Atlantic basin. In contrast, $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{ATL}}$ PC2 yields negative values for all but the intermediate North Atlantic records and does not show strong 100 kyr spectral power. As such, these results suggest that PC2 exhibits the dominant mode of variability recorded in the benthic δ¹³C of North Atlantic waters shallower than 2000 m depth. Curry and Oppo (2005) show that NADW formation to below ∼2000 m is reduced in the North Atlantic during glacial times. The sites with positive factor loadings in PC2 are located at depths <2000 m, and therefore each site should remain consistently bathed in NADW through glacial–interglacial cycles. We thus interpret PC2 as a record of changes in the isotopic values of the North Atlantic carbon reservoir rather than circulation changes.

During MIS 13, all three δ¹³C PC1s (global, Atlantic, and Pacific) demonstrate high positive values. This excursion, first recognized in individual records by Raymo et al. (1997), clearly stands out relative to other δ¹³C interglacial values recorded throughout the last 800 kyr. The MIS 13 excursion is even more apparent when compared against other proxy records such as atmospheric CO₂, SST, and CH₄ (Fig. 7). This high-amplitude change in δ¹³C values is similar to the changes recorded in other proxies during MIS 11, yet precedes the MBT by one glacial cycle. Removal of the MIS 13 interval from variance tests results in no statistical difference in variance before and after the MBT, suggesting a large effect of the carbon isotope excursion on these calculations.

Figure 7 Global ∂¹³C proxy comparison. Comparison of the global ∂¹³C first principal component (PC1; black) to (a) EPICA Dome C CO₂ (yellow; EPICA-community-members, 2004; Lüthi et al., 2008), (b) sea-surface temperature PC1 from this research (blue), and (c) EPICA Dome C CH₄ (red; EPICA-community-members, 2004).

3.5 δ¹³C gradients

Figure 8 shows regional stacks of δ¹³C from the deep (>2000 m) and intermediate (<2000 m) North Atlantic and the deep South Atlantic. As discussed, the intermediate North Atlantic (INA) signal is predominantly controlled by changes in the carbon reservoir over orbital timescales. In contrast, the deep North Atlantic (DNA) is controlled by changes in the relative influence of isotopically more positive NADW and isotopically more negative AABW, as well as any δ¹³C changes to the reservoir that feeds the deep basin from shallower and surficial waters (i.e., INA). Subtracting the INA from the DNA record (i.e., depth gradient) removes the influence of reservoir changes, with the residual time series reflecting only the relative influences of AABW and NADW on the isotopic values of carbon in the deep North Atlantic. This is supported by comparing the North Atlantic depth gradient time series against the South Atlantic stack (Fig. S4). Both time series demonstrate good correlation for the entire time interval (r²=0.58), but even more striking is the similarity in δ¹³C values, with both time series showing similar variability and range in δ¹³C space. The isotopic similarity between the two records suggests adequate removal of reservoir influences with the North Atlantic depth gradient, thus reflecting changes in dominant water mass influence (i.e., circulation). We also note that the correlation between the two records increases starting at MIS 15 (∼530 ka).

Figure 8 Regional ∂¹³C stacks. Stacked records of benthic ∂¹³C separated into three regions: intermediate North Atlantic (orange), deep North Atlantic (red), and deep South Atlantic (blue). All plots are shown in ∂¹³C space to highlight different isotopic values.

Figure 9 MIS 13 and 5e contour plots of ∂¹³C. Contour plots of the ∂¹³C values in the North Atlantic basin for the interglacials MIS 13 and MIS 5e. Red colors represent more positive, enriched values. Blue colors represent lower, depleted values. The plot was created using Ocean Data Viewer.

Figure 10 Latitudinal ∂¹³C gradient. (a) North Atlantic regional ∂¹³C stack plotted in ∂¹³C space (red) authigenic ε_Nd (yellow; Howe et al., 2017). (b) Latitudinal gradient of Atlantic ∂¹³C regional stacks (North Atlantic minus South Atlantic; blue). Lower values demonstrate increased similarity between the records. (c) South Atlantic regional ∂¹³C stack plotted in ∂¹³C space (black). Vertical gray bars indicate glacial periods. Numbers represent Marine Isotope Stage numbers for interglacials.

The depth gradient does not show the prominent MIS 13 excursion that was present in the original DNA stack (Fig. 8), suggesting that the excursion is likely due to a change in the carbon reservoir (represented by the INA) and not related to ocean circulation. Figure 9 shows contour δ¹³C plots of the Atlantic basin for MIS 13 and MIS 5e. Although there is some uncertainty in the these plots due to limited spatial coverage, they show a clear enrichment of the entire basin during MIS 13 relative to average post-MBT interglacial conditions, as represented here by MIS 5e. The global and Pacific δ¹³C PC1s also show the MIS 13 δ¹³C excursion, suggesting a change in the global carbon reservoir.

Next, we evaluate the latitudinal gradient between the South Atlantic signal and the DNA signal in order to further assess the relative influence of the more negative AABW δ¹³C values on North Atlantic δ¹³C values (Fig. 10). Lisiecki (2014) interpreted weaker gradients during glaciations to reflect shoaling of NADW and greater penetration of AABW, which could result from reduced NADW formation or stronger AABW formation. Figure 10b shows a stepwise drop in mean values beginning in MIS 12 (∼436 ka), suggesting a weakening of the gradient due to greater similarity between North Atlantic and South Atlantic glacial and interglacial δ¹³C values.

4 Discussion

Our new analyses demonstrate that there was a statistically significant increase in variance in atmospheric CO₂, Antarctic temperature, global SSTs, and bottom water temperature at 436 ka. These changes are consistent with a transition between two distinct climate states associated with higher-amplitude interglaciations starting with MIS 11, supporting the notion of a MBT as defined by Yin (2013). The same climate variables mentioned above also show an increase in spectral power in the 100 kyr frequency band after the MBT. On the other hand, the dust analyses suggest that the transition to greater variability was experienced in the Southern Hemisphere in the glacial periods starting with MIS 12.

4.1 MIS 13 carbon isotope excursion

The PC1 of $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{G}}$ shows a strong correlation with the CO₂ record for most of the last 800 kyr (Fig. 7a). The exception is during MIS 13, when CO₂ levels were still at pre-MBT levels, while $\mathit{\delta }{{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{G}}$ shows an anomalously high enrichment relative to other interglacial values. This is further illustrated by δ¹³C contour plots showing that the Atlantic basin was enriched in δ¹³C during MIS 13 relative to the MIS 5e (Fig. 9).

We evaluated records of biologic activity in various locations of the Atlantic and Pacific oceans to assess potential sources and sinks in the carbon system during MIS 13. Ba ∕ Fe from the Antarctic zone (AZ) records the sedimentary concentration of biogenic Ba and is thus a proxy of organic matter flux to the deep ocean south of the Polar Front (Jaccard et al., 2013), whereas alkenone concentrations from the Subantarctic zone (SAZ) indicate export productivity to the deep ocean in the region north of the Polar Front (Martínez-Garcia et al., 2009). Based on these proxies, Jaccard et al. (2013) argued that there were two modes of export productivity in the Southern Ocean (SO), where high/low export occurs in the AZ during interglaciations/glaciations, and low/high export occurs in the SAZ during interglaciations/glaciations. They attributed the increase in SAZ export productivity to iron fertilization from increased dust accumulation in the SAZ associated with intensified SO westerlies during glacial periods. Our Southern Hemisphere dust PC1 record supports this hypothesis in showing that high values of dust accumulation correlate with increased values of SAZ export productivity over the last 800 kyr (Fig. S5). We note, however, that the increase in dust starting at MIS 12 does not have an associated decrease in glacial CO₂ values, suggesting that if iron fertilization contributed to lower CO₂ levels, it had an upper limit beyond which additional dust fluxes had little effect.

The antiphase relationship between export productivity between the SAZ and AZ requires a mechanism to increase organic matter productivity in the AZ during interglaciations as suggested by the Ba ∕ Fe signal (Fig. S5c). In the modern SO, vertical mixing and upwelling drive the delivery of nutrient-rich waters necessary for biologic activity to the surface ocean. Wind-driven upwelling is associated with SO westerlies which shift poleward during interglaciations (Toggweiler et al., 2006). Thus, any reduction of upwelling would result from a more northerly position or decrease in strength of the westerlies; a further decrease in nutrient-rich surface waters in the AZ during glaciations likely resulted from increased SO stratification (Sigman et al., 2010; Jaccard et al., 2013). We note, however, that Jaccard et al. (2013) find no AZ export productivity during MIS 13, whereas all other interglaciations over the last 800 kyr show some evidence for it (Fig. S5c). This skipped interglaciation in export productivity suggests some combination of a change in the position/strength of the SO westerlies or stratification of the AZ that limited the delivery of nutrient-rich deep waters to the surface as compared to other interglaciations of the last 800 kyr.

Figure 11 Marine Isotope Stages 15 to 13 and the carbon isotope excursion. (a) Summer insolation at 65^∘ N (red; Laskar et al., 2004). (b) First principal component of Atlantic ∂¹³C (black). (c) EPICA Dome C CO₂ (yellow; EPICA-community-members, 2004; Lüthi et al., 2008). (d) Detrended sea-level equivalent from Shakun et al. (2015) (blue). Derived from $\partial {{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{sw}}$ calculations. Negative numbers indicate lower sea level and increased ice volume. (e) Chinese Loess Plateau grain size indicating relative Asian summer monsoon strength (brown; Sun and An, 2005). (f) Quartz/calcite ratios from site U1313 in the North Atlantic as a measure of ice-rafted debris (light blue; Naafs et al., 2012). Dark gray bars highlight the interglacials (MIS 15 and MIS 13) between ∼630 and ∼470 ka. The light gray bar highlights MIS 14.

The PC1s of δ¹³C (global, Atlantic, and Pacific) demonstrate that the global ocean was enriched in heavy carbon during MIS 13 relative to any other interglaciation of the last 800 kyr (Fig. 3). In contrast, atmospheric CO₂ concentrations were ∼ 240 ppm during MIS 13, similar to other pre-MBT interglacial levels (Fig. 1). Ba ∕ Fe records of organic export productivity from the AZ that acts as a sink for light carbon indicate no increase during this interglaciation, while Ca/Al records from the SAZ indicate increased preservation and thus a deeper lysocline and lower dissolved inorganic carbon (Jaccard et al., 2010). The question thus becomes the following: if the ocean is heavily enriched in δ¹³C during MIS 13 while CO₂ and export productivity remained at low levels, what reservoir contained the isotopically light carbon?

Paleoclimate records from the CLP indicate greater precipitation during MIS 13 relative to the other interglaciations (Liu, 1985; Yin and Guo, 2008). This greater precipitation has been attributed to increased monsoon activity recognized throughout monsoonal areas of the Northern Hemisphere and persisting through MIS 15, 14, and 13 (Yin and Guo, 2008; Guo et al., 2009). Biogenic silica measurements from Lake Baikal exhibit continuously high terrestrial productivity in central Asia throughout MIS 15 to MIS 13 (Prokopenko et al., 2002), whereas sea-level reconstructions indicate that ice volume during MIS 14 was considerably less relative to other glacial maxima of the last 800 kyr (Fig. 11d) (Elderfield et al., 2012; Shakun et al., 2015). Thus, the smaller ice sheets of MIS 14 would likely have had a lesser effect on displacing forested areas of the Northern Hemisphere, allowing greater terrestrial carbon storage to potentially persist through a glacial cycle (Harden et al., 1992). We thus suggest that the increased monsoonal precipitation and smaller ice volume during MIS 14 would have combined to increase land biomass that continued into MIS 13. The Northern Hemisphere thus had the potential to store light carbon in the terrestrial reservoir resulting in the enriched δ¹³C MIS 13 signal seen in the ocean basins (Yin and Guo, 2008).

4.2 Ocean circulation changes in the Atlantic basin

One explanation for the glacial–interglacial variations in atmospheric CO₂ invokes a dominant role by the Southern Ocean in storing and releasing dissolved inorganic carbon (DIC) in the deep Southern Ocean, with deep-ocean sequestration of atmospheric CO₂ occurring through decreased upwelling and vertical mixing of AABW (Sigman et al., 2010). Expansion of Southern Ocean sea ice can also lower atmospheric CO₂ by insulating upwelled water from the atmosphere, thus reducing outgassing, and by increasing the volume of AABW and its capacity to hold DIC (Stephens and Keeling, 2000; Ferrari et al., 2014). According to this framework, pre-MBT interglaciations with lower CO₂ would be associated with greater sea-ice extent and a larger volume of AABW, whereas post-MBT interglaciations with higher CO₂ suggest reduced sea-ice extent and AABW volume. Glacial values of CO₂ remain relatively constant throughout the last 800 kyr (Fig. 1), suggesting that the change in relative AABW volume before and after the MBT only occurred during interglaciations.

Figure 12 Average interglacial ∂¹³C contours. Contour plots of the average interglacial ∂¹³C values in the Atlantic for (a) pre-MBT included MIS 13, (b) pre-MBT excluding MIS 13 (enriched carbon isotope excursion), and (c) post-MBT. Red colors indicate higher ∂¹³C values. Blue colors indicate lower ∂¹³C values. The boundary between the two water masses (NADW and AABW) is indicated at the 0.25 ‰ contour (Curry and Oppo, 2005).

This mechanism is consistent with ice-core evidence for greater sea-ice extent during pre-MBT interglaciations (Wolff et al., 2006) and with modeling results that show that interglacial AABW formation decreased after the MBT through insolation-induced feedbacks on sea ice and surface water density (Yin, 2013). Moreover, based on the Ba ∕ Fe proxy of organic matter flux to the deep ocean south of the Polar Front, Jaccard et al. (2013) argued that the deep Southern Ocean reservoir was larger prior to the MBT.

Our analyses of changes in Atlantic δ¹³C over the last 800 kyr further support an important role of AABW in causing the post-MBT increase in interglacial CO₂. In particular, the steeper latitudinal gradient between North and South Atlantic δ¹³C records before the MBT reflects greater northward penetration of AABW, whereas the post-MBT decrease in gradient suggests greater southward penetration of NADW (Fig. 10b). These gradient changes are further illustrated by contour plots of average interglacial δ¹³C values in the Atlantic which show that prior to the MBT, AABW penetrated north of the Equator, increasing the δ¹³C gradient (Fig. 12a), in contrast to remaining south of the Equator after the MBT, decreasing the gradient (Fig. 12c). Removal of MIS 13 and its associated enriched carbon isotope excursion further highlights the greater volume of AABW in the pre-MBT interglacial Atlantic (Fig. 12b). We note that a record of the water mass tracer ε_Nd from 6^∘ N (Howe and Piotrowski, 2017) is in good agreement with our North Atlantic regional δ¹³C stack (Fig. 10a), with both records suggesting that changes in volume of the interglacial AABW occurred south of the Equator. This reorganization of the dominant interglacial water masses in the Atlantic basin across the MBT, perhaps resulting from insolation-induced feedbacks (Yin, 2013), would lead to a greater release of deep-ocean CO₂ during the post-MBT interglaciations, with corresponding warmer interglaciations (Fig. 5). An alternative explanation for the observed decrease in latitudinal gradient could be changes in the isotopic composition of AABW across this time period. However, modeling results of long-term carbon fluctuations across this interval suggest that changes in the burial rate of organic and inorganic carbon caused the δ¹³C depletion – the opposite signal necessary to create the increased similarity between northern- and southern-sourced waters (Hoogakker et al., 2006). Thus, it is more likely explained by changes in AABW influence north of the Equator.

Cross-spectral analysis of pre-MBT North and South Atlantic δ¹³C stacks indicates in-phase coherency between the records at the eccentricity and obliquity frequencies. Similar tests for the post-MBT δ¹³C stacks exhibit coherency at eccentricity, obliquity, and precession frequencies, with the South Atlantic stack leading the North Atlantic by ∼23^∘ (7 kyr) in eccentricity, ∼18^∘ (2 kyr) in obliquity, and ∼36^∘ (2 kyr) in precession (Fig. S6). All phase relationships overlap within uncertainty, suggesting that South Atlantic δ¹³C leads North Atlantic δ¹³C by 2–7 kyr following the MBT. This lead by the South Atlantic is most apparent during terminations (Figs. 9, 12) and is most likely related to deglacial mechanisms for ventilation of respired CO₂ from the deep Southern Ocean such as enhanced wind-driven upwelling or the melting of sea ice in response to the bipolar seesaw (Cheng et al., 2009).

5 Conclusions

Using statistical analyses of multiple climate proxies, we have further characterized the Mid-Brunhes Transition as an increase in interglacial sea-surface and Antarctic temperatures, atmospheric CO₂, and CH₄ beginning with MIS 11. At the same time, our new analyses also document a number of changes in other components of the climate system that began as early as MIS 14 that suggest a more complex sequence of events prior to the MBT, although their relationship to the MBT remains unclear. Figure 13 highlights key features in the sequence of events beginning with an increase in Asian summer monsoon strength during MIS 15 that persisted through MIS 14 and into MIS 13. The strong monsoon strength during MIS 14 is associated with a weak glaciation, which in combination would have been conducive to a build-up of Northern Hemisphere land biomass. A continued strong Asian summer monsoon during MIS 13 associated with greater precipitation would have further sequestered land biomass and provided a reservoir for light carbon, resulting in the oceans becoming unusually enriched in δ¹³C as recorded in the global benthic δ¹³C carbon isotope excursion. MIS 12 was associated with the return of large ice sheets, collapse of the Asian summer monsoon, and the first increase in amplitude of Southern Hemisphere dust. A decrease in the latitudinal gradient of interglacial Atlantic δ¹³C at the MBT suggests a reorganization of the water masses in the basin and reduction in the size of interglacial AABW, thus possibly explaining the increase in interglacial values of atmospheric CO₂ with corresponding increases in interglacial SSTs and CH₄. This evidence for a change in AABW is consistent with modeling results that suggest that the MBT was forced by insolation (Yin, 2013).

Figure 13 Schematic representation of the sequence of events leading to the Mid-Brunhes Transition. Corresponding Marine Isotope Stages are located on the left side of each row. Boxes in a row indicate synchronous events.