In and out of glacial extremes by way of dust-climate feedbacks.

Mineral dust aerosols cool Earth directly by scattering incoming solar radiation and indirectly by affecting clouds and biogeochemical cycles. Recent Earth history has featured quasi-100,000-y, glacial-interglacial climate cycles with lower/higher temperatures and greenhouse gas concentrations during glacials/interglacials. Global average, glacial maxima dust levels were more than 3 times higher than during interglacials, thereby contributing to glacial cooling. However, the timing, strength, and overall role of dust-climate feedbacks over these cycles remain unclear. Here we use dust deposition data and temperature reconstructions from ice sheet, ocean sediment, and land archives to construct dust-climate relationships. Although absolute dust deposition rates vary greatly among these archives, they all exhibit striking, nonlinear increases toward coldest glacial conditions. From these relationships and reconstructed temperature time series, we diagnose glacial-interglacial time series of dust radiative forcing and iron fertilization of ocean biota, and use these time series to force Earth system model simulations. The results of these simulations show that dust-climate feedbacks, perhaps set off by orbital forcing, push the system in and out of extreme cold conditions such as glacial maxima. Without these dust effects, glacial temperature and atmospheric CO2 concentrations would have been much more stable at higher, intermediate glacial levels. The structure of residual anomalies over the glacial-interglacial climate cycles after subtraction of dust effects provides constraints for the strength and timing of other processes governing these cycles.

Mineral dust aerosols (hereafter called dust) are an integral part of Earth’s climate system. They have a direct effect by absorbing and scattering incoming solar and outgoing terrestrial radiation (1, 2) and an indirect effect by affecting cloud formation (3, 4) and biogeochemical cycles (5). A main indirect effect is the fertilization of ocean phytoplankton production by dust-mitigated input of iron to the ocean surface layer (6, 7). Dust impact on climate has varied over time, depending on the amount of dust particles in the atmosphere and therefore on aspects like vegetation, wind strength, and precipitation patterns (8–10). Dust impact is also spatially heterogeneous. Since the major dust sources are located in northern subtropical latitudes, the direct effect is mainly concentrated in these regions (11), although polar regions are sensitive to amplified high-latitude variability in dust concentrations (12). Potential carbon dioxide (CO2) drawdown through dust-induced iron fertilization of ocean biota is restricted to high-nutrient, low-chlorophyll (HNLC) regions, mainly the Southern Ocean and northern Pacific where iron is the limiting micronutrient (6, 13). This indirect effect of dust on climate via CO2 thus increases the climatic importance of secondary dust sources like South America and Oceania.

Over the past 800,000 y, Earth’s climate has varied over ∼100,000-y cycles of glacial−interglacial stages, with higher and lower temperatures and greenhouse gas concentrations during interglacials and glacials, respectively (14–16). Glacial global mean temperatures (GMT) dropped by about 4 ± 0.8 K while atmospheric pCO2 (CO2 partial pressure) decreased by 80 ppm to 100 ppm (16, 17). In contrast, glacial maxima dust levels were, on average, more than three times higher than interglacial levels (10, 18, 19), but with much higher ratios in high latitudes (ref. 12 and Table S1). Most of the GMT drop has been attributed to radiative forcing decreases from increased albedo due to equatorward ice extension and from decreased greenhouse gas concentrations; vegetation and atmospheric dust are thought to play secondary roles (20, 21). About half of the pCO2 decrease may be due to increased glacial ocean stratification, trapping carbon-rich waters in the deep layers away from the atmosphere (22, 23). Another ∼20-ppm reduction may result from increased atmospheric dustiness that would increase iron input to HNLC areas of the ocean, enhancing biological pump efficiency (24, 25). Other effects like temperature-dependent CO2 solubility in ocean water, carbon stored in the land biosphere, weathering rates, and ocean nutrient inventories may help explain the rest of glacial−interglacial changes in atmospheric pCO2 (26, 27). The increased atmospheric dust load during glacial times is thought to be due to a generally colder and drier atmosphere, which increased the number and strength of dust sources through reduced vegetation cover (28), and reduced washout of suspended particles during transport (7).

Results

Dust Deposition Data.

Observational Holocene and Last Glacial Maximum (LGM) dust flux data were recently compiled to create global dust deposition maps for these time slices (25). Fig. 1 shows zonal means and medians calculated from the results used to construct these maps. Holocene and LGM mean, medians, and ratios from these data are also given in Table S1. Northern Hemisphere (NH) dust fluxes typically exceed Southern Hemisphere (SH) ones by factors of 5 or more, reflecting north−south partition of dust source areas and limited cross-equatorial transport (Fig. 1). However, ratios of LGM to Holocene fluxes show considerable hemispheric symmetry, with ratios close to 3 from about 60°S to 40°N (Fig. 1). Likewise, both hemispheres exhibit strong LGM dust flux amplification at high latitudes.

Fig. 1.

Zonal mean and median dust deposition for the Holocene (HOL) and the LGM expressed as (A) rates and (B) ratios of these rates, calculated from the results presented in ref. 25. The blue points in B are LGM/HOL ratios for the individual dust deposition time series considered in Fig. 2. Also indicated are the latitudes for the Antarctic (EDC) and Greenland (NGRIP) temperature anomaly time series in Fig. 2.

Fig. 2.

SH and NH time series of temperature and scaled dust deposition rates for the last three glacial cycles. Shown are (A) Antarctic temperature anomaly relative to a Holocene mean (14); (B) Antarctic dust flux (ref. 56; brown line), scaled by the mean of the record and Southern Ocean, 232Th-normalized dust fluxes from the South Atlantic (ref. 30; orange line) and South Pacific (ref. 31; maroon line), scaled by interglacial values; and (C) Greenland temperature anomaly relative to a Holocene mean. Values from −122 ky to 0 ky are from ref. 32; values from −340 ky to −122.4 ky have mainly been estimated using Antarctic methane observations and an Antarctic methane−Greenland temperature relationship (ref. 33; see ), and (D) Chinese Loess dust fluxes scaled by interglacial values (29). This record is shown on the original reported time scale (dashed line) and on a time scale stretched by 4% (solid line) so as to better match up with the Greenland temperature data (vertical arrows). Thin lines in A–C are low-passed filtered time series used in the subsequent analyses (see ). All ice core records are on the Antarctic ice core chronology (AICC2012) timescale (57).

Here we deal with the time-varying role of radiative dust forcing and of indirect dust forcing from ocean iron fertilization over glacial/interglacial climate cycles. For this, we sought to identify specific, well-resolved, published dust flux time series that may be considered representative for each of these forcings. As shown in Fig. 1, the NH dominates global dust variability, and maximum radiative forcing is found between about 20°N to 60°N. We therefore chose a dust record (Fig. 2) from the Loess Plateau of northern China (ref. 29; Xifeng: 35.7°N, 107.7°E) as our proxy for the relative time variation of glacial/interglacial dust radiative forcing. The LGM/Holocene dust flux ratio from this record mirrors the zonal mean at this latitude (upper blue dot in Fig. 1), further supporting the representativeness of our choice.

To estimate the role of time-varying, ocean iron fertilization, we turned to the Southern Ocean, the main HNLC, iron-limited region of the global ocean. We chose two 232Th-normalized dust flux records from sediment cores (Fig. 2), one each from the South Atlantic (ref. 30; core Ocean Drilling Program 1090; 42.9°S, 8.9°E) and South Pacific (ref. 31; core PS75/76; 55.5°S, 156.1°W) to represent relative iron fertilization forcing variations in these areas. While South Atlantic dust fluxes are somewhat larger than the South Pacific ones (31), both records exhibit about the same LGM/Holocene dust flux ratios, and both ratios mirror zonal means at their respective latitudes (middle two blue dots in Fig. 1). This allows us to scale the two records by their respective interglacial levels, combine them to better resolve the Southern Ocean, and use the combined record as our proxy for glacial/interglacial iron fertilization forcing. Although mass accumulation rates from these two sediment cores extend much farther back in time, the 232Th normalization of dust fluxes that are useful for our purposes here can only be calculated for the last three glacial/interglacial cycles (24, 31). This sets the time scale for the analyses we report here. Also plotted in Fig. 2 is the scaled, high-resolution, East Antarctic dust record (14) from the European Project for Ice Coring in Antarctica, Dome C (EDC) ice core (75.1°S, 123.3° E).

Since dust fluxes are expected to be climate-dependent with increased fluxes for cooler, drier glacial conditions (7, 28), we analyze the proxy time series discussed above together with relevant temperature time series to search for dust−climate relationships. Since dust radiative and iron fertilization forcing are mainly concentrated in the NH and the SH, respectively, we use respective NH and SH temperature time series for this. For the SH, we chose the high-resolution Antarctic record (ref. 14 and Fig. 2) derived from EDC ice core data. There is strong evidence for a tight relation of Southern Ocean iron fertilization to this record (30, 31). For the NH, we chose a high-resolution Greenland record (32) derived from the North Greenland Ice Core Project (NGRIP) ice core data (75.0°N, 42.3°W) back to 122.4 ky B.P. and extended it farther back to 340 ky B.P., mainly using Antarctic methane observations and an Antarctic methane−Greenland temperature relationship for the period 122.4 ky B.P. to 0 ky B.P. (ref. 33 and ). This temperature time series is shown in Fig. 2. Multimillenium-scale climate variability in Eastern Asia, and indeed throughout NH subtropics and midlatitudes, has been shown to be well correlated with such variability in Greenland (34–36).

Derived Dust Deposition−Temperature Relationships.

Visual inspection of Fig. 2 shows maximum dust fluxes for maximum glacial cooling. However, for intermediate glacial conditions, for example, for the period −110 ky to −75 ky, dust fluxes are much lower, approaching low interglacial levels. This nonlinear tendency can be captured quantitatively for the high-resolution Antarctic results (Fig. 2 ) for which dust deposition and temperature estimates are derived from the same ice core. We found a good best fit to dust deposition (D) vs. temperature anomaly (ΔT) for these results (Fig. 3) using the function(dashed line in Fig. 3, and Table S2).

Fig. 3.

Scatter plots of dust deposition flux versus temperature for SH and NH records plotted in Fig. 2. (A) Antarctic dust flux scaled by an interglacial mean and plotted against Antarctic temperature anomaly (purple dots). The dashed line is the best fit to this data using Eq. . The black dots with error bars in A–C are four-bin means with 1-sigma SD (see ), and the solid black lines are best fits to these means using Eq. . (B) Scaled Southern Ocean (SO) dust flux (combined from both records shown in Fig. 2) plotted against Antarctic temperature anomaly (blue dots). (C) Scaled Chinese Loess (CL) dust flux plotted against Greenland temperature anomaly on the stretched time scale of Fig. 2 (red dots). (D) Plots like A−C but with original (unscaled) dust flux data plotted on a logarithmic scale. Also shown in D are global mean and median dust fluxes for the LGM and the Holocene (black symbols to the left and to the right in the frame, respectively), based on the results in ref. 25 and tabulated in the first column of Table S1.

This good fit encourages us to seek dust deposition−temperature relationships by applying Eq. to the other data from Fig. 2. Fig. 2 shows individual dust deposition time series from the Southern Ocean and Chinese Loess that may be compared with Antarctic and Greenland temperature anomalies, respectively (Fig. 2 ). Time scale issues may arise when plotting data together from different data archives and geographical positions, and we here seek to minimize such issues. The Southern Ocean dust data have already been put on an EDC time scale (31). We find that, by stretching the original Chinese Loess time scale by 4%, we obtain a better matchup with the NGRIP/Greenland time scale (Fig. 2 ). The dust−temperature relationships for the Southern Ocean and Chinese Loess are shown in Fig. 3 . To further minimize possible issues connected with time scales, data uncertainty, and uneven data distribution across temperature ranges, we collect the dust deposition−temperature data points for each data set combination in Fig. 3 into four bins (see ). We then calculate the means (black dots in Fig. 3 ) and medians for each bin and fit Eq. to the bin means (solid black lines in the figure) and medians. Some properties of these fits are listed in Table S2. Also listed are properties of fits for sensitivity studies using the original dust deposition−temperature data points, data collected into 16 rather than 4 bins, four-binned data on the original Chinese Loess time scale (29), and four-binned data for an alternative Greenland temperature series (). For the high-resolution Antarctic data, the four-bin-based fit and the fit to all data are essentially identical (solid and dashed lines in Fig. 3, and Table S2), lending support to our binning procedure.

The binned means for the Southern Ocean data provide the basis for a well-defined fit using Eq. (Fig. 3). The binned means for the Chinese dust data show more scatter with respect to an Eq. fit (Fig. 3). Most of this scatter results from relatively elevated dust fluxes for interglacial temperatures (rightmost black dot in Fig. 3). This result can be traced to elevated Holocene dust fluxes compared with earlier interglacial periods (Fig. 2), perhaps related to Holocene human activity (37). In general, it is striking that climate relationships for dust fluxes spanning the globe and five orders of magnitude (Fig. 3) appear to follow similar nonlinear relationships with greatly enhanced fluxes for extreme glacial temperatures.

Dust Radiative and Iron Fertilization Forcings.

We use the above dust deposition−temperature relationships to derive climate dependencies of dust radiative and iron fertilization forcing under the assumption that both forcings vary in proportion to the dust deposition, and can therefore be expressed in the same form as Eq. . First, we will address dust radiative forcing, RF, relative to Holocene forcing. Fig. 1 shows that RF is concentrated in the NH, and we therefore relate RF to Greenland temperature and dust fluxes derived from Chinese Loess. Specifically, we here use low-passed filtered, Greenland temperature deviation relative to a Holocene mean, ΔT, as plotted in Fig. 2. As so defined, both RF and ΔT are zero for interglacial conditions, and application of Eq. yields a + c = 0. Then, for LGM conditions with RF = RF and ΔT = ΔT, application of Eq. yields a*exp(−b* ΔT) + c = RF. After solving for a and c from these two conditions, we obtain the desired RF expression,where ΔT is chosen to be −24 °C from Fig. 2 and values for b are taken from the fits of Eq. to the four-binned results from the Chinese Loess data (Table S2). Published estimates for RF vary considerably but center upon about −1 W⋅m−2 (38–40), a value we adopt here for our standard case.

For iron fertilization forcing, we choose temperature−dust relationships from the high southern latitudes, as the Southern Ocean is the main region where this process is relevant. We deal with iron fertilization in the context of the Danish Center for Earth System Science (DCESS) model (41) for which reduction of high-latitude new production (relative to that which would occur if phytoplankton there could make full use of all available nutrients) is expressed in terms of an efficiency factor (see equation 19 in ref. 41). For present, interglacial conditions, the value of this factor, here referred to as the iron fertilization factor (IF), was calibrated to be 0.36 by fitting to observed, high-latitude phosphate concentrations (41). As motivated above, we relate iron fertilization forcing, and therefore IF, to Antarctic temperature, specifically here to the low-passed filtered, Antarctic temperature deviation relative to a Holocene mean, ΔT, as plotted in Fig. 2. As so defined, IF = 0.36 and ΔT = 0 for interglacial conditions, and application of Eq. yields a + c = 0.36. For LGM conditions with IF = IF and ΔT = ΔT, application of Eq. yields a*exp(−b*ΔT) + c = IF. After solving for a and c from these two conditions, we obtain the desired IF expression,where ΔT is chosen to be −8 °C from Fig. 2 and values for b are taken from the fits of Eq. to the four-binned results from the Southern Ocean data (Table S2). In a recent DCESS model application (25), IF needed to be increased to 0.50 to emulate the initial LGM pCO2 drawdown in a complex ocean biogeochemical model forced by observed LGM dust deposition. Here we adopt this value as our standard case IF.

The use of Eqs. and for the standard case b values based on four-bin means and the low-passed filtered, temperature anomaly time series from Greenland and Antarctica (Fig. 2 ) leads to radiative and iron fertilization forcing time series over the last three glacial cycles (Fig. 4). As sensitivity studies, we also calculated these time series based on (i) different values of b obtained by fitting Eq. to four-bin medians (Table S2), (ii) different values of b obtained by fitting Eq. to the raw data and to 16-bin means and medians (Fig. S1 and Table S2), (iii) different values of b for radiative forcing obtained by fitting Eq. to bin means and medians using an alternative low-passed filtered, Greenland temperature anomaly time series () and application of that time series to construct alternative radiative forcing time series, (iv) radiative forcing calculated for 50% decrease/increase compared with our standard LGM value (RF = −0.5 and −1.5 W⋅m−2), and (v) iron fertilization forcing calculated for 50% decrease/increase of the difference between standard LGM and present-day values (IF = 0.43 and 0.57). Forcing time series for these alternatives are plotted in Fig. S2.

Fig. 4.

Reconstructed dust forcing time series for the last three glacial cycles. Shown are (A) dust radiative forcing (RF) anomaly relative to preindustrial and (B) iron fertilization parameter (IF) for ocean new production in the DCESS model (41) high-latitude zone. This is how the model deals with iron fertilization strength in the Southern Ocean (SO; see Dust Radiative and Iron Fertilization Forcings).

Forced Simulations.

We carried out a number of DCESS model simulations over the last three glacial cycles using the dust radiative and/or iron fertilization forcings shown in Fig. 4 and Fig. S2. The radiative forcing is applied to model low-latitude to midlatitude and high-latitude atmospheric boxes, since large dust depositions extend into high NH latitudes during glacial periods (ref. 12 and Fig. 1). Iron fertilization forcing is applied to the surface layer of the high-latitude ocean that represents the Southern Ocean in the model (). Also, in preparation for our simulations, we made model hindcasts for a range of climate sensitivities and forced by the estimated total radiative forcing anomaly for the period AD 1765–2012 (Fig. S3). Compared with observed atmospheric and ocean warming, the hindcasts tracked the observations best in both atmosphere and ocean for a pCO2-doubling, climate sensitivity of 2.5 K. This value is also consistent with recent work on the LGM (21), and we use the model version with this sensitivity in our simulations here for simplicity (but see ).

Fig. 5 shows model simulation results for the combined radiative and iron fertilization forcing derived using standard case, four-bin b values (Fig. 4). Here we present results for the last glacial cycle only, for which direct temperature estimates are available from both the Greenland and Antarctic ice cores such that the results are likely more reliable. We also made combined radiative and iron fertilization forcing simulations over the last three glacial cycles (Figs. S4 and S7) for all forcing time series shown in Fig. S2 and found results similar to those reported below. The pCO2 drawdown from the combined forcing exceeds 25 ppm at the LGM and only comes into play for the coldest glacial conditions (Fig. 5). Most model drawdown is associated with higher new production in the high-latitude zone, due to iron fertilization, combined with carbonate compensation (ref. 25; red line in Fig. S5). However, some drawdown can be explained by increased CO2 solubility in cooler ocean surface layers, a product of dust radiative forcing (blue line in Fig. S5). A comparison of this result with the observed pCO2 anomaly (42) and the residual from subtraction of this result from the anomaly shows that the lowest pCO2 values during the coldest glacial times stem from the dust−climate feedbacks considered here (Fig. 5). Other processes must be responsible for the 40- to 50-ppm pCO2 drawdown during the initial transition from interglacial to glacial conditions as well as for a comparable pCO2 increase during the latter part of glacial terminatins.

Fig. 5.

Comparison with data and data-based reconstructions over the last glacial cycle of a model simulation for combined dust radiative and iron fertilization forcing. (A) Simulation pCO2 anomaly from 278 ppm, preindustrial level (blue line), pCO2 anomaly from observed, Antarctic ice core pCO2 (ref. 42; black line), and the residual anomaly after subtracting model results from observations (red line). (B) Simulated anomaly from a Holocene mean of GMT (blue line) and residual anomaly after subtracting model results from a GMT anomaly reconstruction (red line). The GMT anomaly reconstruction is based on sea surface temperature reconstructions with 3-ky resolution (ref. 43; black line), adjusted for agreement with a comprehensive GMT reconstruction at the LGM (ref. 17; blue asterisk).The dotted line is daily mean, summer incoming radiation at 65°N.

Global cooling that results from the combined dust forcing exceeds 1 K at the LGM and is only significant for the coldest glacial conditions (Fig. 5). Most of the model cooling is associated with dust radiative forcing, including the ocean cooling effect explained above (blue line in Fig. S5). However, some of the cooling stems from decreased atmospheric greenhouse gas forcing from iron fertilization-induced pCO2 drawdown (red line in Fig. S5). We compare simulation cooling from the combined forcing to a GMT reconstruction, based on a global mean sea surface temperature reconstruction (43) that we scaled by the factor 1.84 so as to match estimated LGM cooling (ref. 17; 19 ky B.P. to 23 ky B.P.). We then construct residual cooling without dust by subtracting simulated cooling from the GMT reconstruction. Results show that the coldest glacial maxima temperatures are a product of the dust−climate feedbacks considered here (Fig. 5). Initial cooling during glacial inception appears to be coupled to a strong negative anomaly of summertime solar radiation at high northern latitudes from orbital forcing (Fig. 5). However, temperatures do not drop low enough then to induce a strong dust−climate feedback. However, as the next major negative radiation anomaly occurs around 70 ky B.P., the system starts out cold enough for the anomaly to induce a strong dust−climate feedback leading to extreme cold conditions around 65 ky B.P. (Fig. 5). Such a direct connection between orbital forcing and dust−climate feedbacks is not as clear for explaining the apparent dominant role of dust−climate feedbacks in LGM cooling.

Discussion

Our results show that dust−climate feedbacks can explain the final push into extreme glacial conditions for both GMT and pCO2, thereby explaining about one-fourth of the total interglacial−glacial change for both properties. These conclusions are found to hold up for a range of alternate analyses and forcing function options (Figs. S2, S4 and S7). Furthermore, our results indicate that these feedbacks play an important role in the initial stages of glacial terminations.

Models for simulating glacial−interglacial climate and/or CO2 evolution have often formulated dust forcing using Antarctic dust observations (39, 44–46). However, large glacial amplification of the Antarctic dust signal (ref. 12 and Fig. 1) reflects transport processes, making the direct, far-afield application of these observations problematic. Here we took a different approach by seeking and applying dust deposition observations from hotspots for dust radiative and iron fertilization forcing: NH subtropical latitudes and the Southern Ocean. Observations from these regions show much less glacial amplification than near the poles but do exhibit similar nonlinear relationships with greatly enhanced fluxes for extreme glacial temperatures (Figs. 1 and 3).

Here we used a constant equilibrium climate sensitivity of 2.5 K, determined by an excellent fit of DCESS model simulations for this sensitivity to modern-day observations of both atmosphere and ocean warming (Fig. S3). Some recent analyses have indicated lower climate sensitivity during glacial maxima than during interglacials but with interglacial values considerably higher than 2.5 K (47, 48). We have more confidence in our present-day (interglacial) result since it is based on fits to direct observations rather than on more uncertain paleoreconstructions, and we find it unlikely that glacial sensitivities were much lower than 2.5 K, given ice albedo feedbacks (47, 48). Recent best estimates for LGM cooling range from −2.7 K to −6.2 K (17, 21, 47–50). For our comparisons here (Fig. 5 and Figs. S4 and S7), we use the value of −4 K from ref. 17. We prefer this result since it was based on a comprehensive analysis using extensive proxy data from both land and ocean in combination with an ensemble of climate model simulations, in a paper focused entirely on LGM cooling. We note, however, that our main conclusion that dust−climate feedbacks can explain the final push into extreme glacial cooling is insensitive to the particular LGM cooling choice. Our choices for climate sensitivity (2.5 K) and glacial maxima cooling (−4 K) are both on the low side for such properties in other work cited above. However, these values are internally consistent in that commensurate increases in them would still yield results supporting our specific conclusions on the role of dust−climate feedbacks. For example, when we reran our standard case model using a climate sensitivity of 3 K, glacial maxima results showed a 7% increase in dust-driven, pCO2 drawdown but a 35% increase in dust-driven cooling. Such cooling would represent about one-fourth of overall glacial maxima cooling if that overall cooling were −5 K.

Residual GMT and pCO2 anomalies after subtraction of dust effects (red lines in Fig. 5) invite further interpretation of the strengths and timings of other mechanisms acting over glacial−interglacial cycles. GMT drops initially at glacial inception in response to decreased summer radiation at high northern latitudes that would have led to equatorward extension of sea ice and snow cover with associated cooling from increased albedo. Somewhat later, there is a large and rapid pCO2 drop that further enhances the cooling. The speed and size of this initial drop, at least one-third of the overall interglacial/glacial change, and a mirror image pCO2 rise in the middle of glacial terminations point toward relatively abrupt decreases/increases in deep ocean ventilation (22, 23) as important brackets for the glacial inception/termination cycle. The remaining slow drift to lower GMT and pCO2 over glacial time, punctuated by higher-frequency variability and the dust−climate feedbacks, may reflect the consequences of the growth of continental ice sheets via albedo increases (also from vegetation changes) and increased CO2 dissolution in the ocean from cooling.

Materials and Methods

Direct temperature estimates for the NGRIP Greenland ice core used here are only available back to 122.4 ky B.P. (32). For earlier times, we adopt Greenland temperature estimated as follows (33): For the period 128,700 B.P. to 340,000 B.P., this temperature was derived from a proxy based on Antarctic ice core methane data using the relation T = −51.5 + 0.0802 [CH4(ppb)] from a linear regression of Greenland temperature estimates on Antarctic methane for the period 150 B.P. to 122,400 B.P. For the remaining period of 122,400 B.P. to 128,700 B.P., data from a variety of climate archives indicate that Greenland warming lags that of Antarctica, with rapid warming commencing around 128.5 ky B.P. in the northern North Atlantic and reaching full interglacial levels by about 127 ky B.P. (51). Guided by these results, the interglacial Greenland ice core temperature at 122.4 ky B.P. was extended back to 127 ky B.P., and a linear interpolation was applied between that value then and the (methane-based) temperature at 128.7 ky B.P.

We applied a low-passed filter to the original, ice core-based temperature and dust time series of Fig. 2 and then used these low-passed filtered data in the subsequent analyses. For this, we applied a fifth-order, low-pass Butterworth filter with cutoff period of 3,000 y. This was chosen to filter out millennium-scale variability in the records to concentrate on processes we deal with here. To help minimize potential issues with time scale mismatches, data uncertainty, and uneven data distribution across temperature ranges, we collected the dust deposition data points for each data set combination in Fig. 3 into four bins. For the NH dust data, we defined the following bins of Greenland temperature anomaly, ΔT, for data grouping: ΔT ≥ −6 °C, −12 °C ≥ ΔT < −6 °C, −18 °C ≥ ΔT < −12 °C, and ΔT < −18 °C. Bin means (and medians) were then assigned to the temperatures ΔT = 0, −9, −15 and −21 °C, respectively, for plotting and fits using Eq. . For the SH dust data, we defined the following bins of Antarctic temperature anomaly, ΔT, for data grouping: ΔT ≥ −2 °C, −4 °C ≥ ΔT < −2 °C, −6 °C ≥ ΔT < −4 °C, and ΔT < −6 °C. Bin means (and medians) were then assigned to the temperatures ΔT = 0, −3, −5 and −7 °C, respectively, for plotting and fits using Eq. . Some results of these fits are plotted in Fig. 3 and listed in Table S2. Furthermore, as a sensitivity study, we redid the binning procedure and fits to Eq. but using 16 bins defined in a way analogous to the above ().

We use the DCESS model, a relatively simple but well-tested and flexible Earth system model of intermediate complexity (41, 52, 53). It features components for the atmosphere, ocean, ocean sediment, land biosphere, and lithosphere and has been designed for global climate change simulations on time scales from years to millions of years. See for a short model description. For the present application, we hold volcanic input and weathering constant to preindustrial values (41). Weathering may depend on properties like temperature, pCO2, land ice abrasion, and sea level change but has been suggested to be stable over glacial−interglacial cycles (54). Furthermore, we hold land biosphere fluxes and reservoir sizes constant to preindustrial values (41). The original DCESS land biosphere module is too simple to deal correctly with glacial−interglacial cycles. Furthermore, there is a tendency toward approximate compensation over these cycles of loss/gain of land biosphere carbon with gain/loss of permafrost carbon (55).

Our simulations over three glacial cycles (340 ky B.P. to 0 B.P.) were started using preindustrial initial conditions (41). There was a slight pCO2 drift over the simulations due to slight net imbalances between (constant-held) weathering and volcanic net inputs and net carbon outputs from (calculated) carbonate and organic carbon burial down out of the ocean sediment. As a consequence, pCO2 levels at the preindustrial end of all of the simulations were slightly higher (by 9.5 ppm to 10.6 ppm) than the observed value of 278 ppm such that GMT values there were also slightly higher (by 0.23 °C to 0.25 °C) than the standard DCESS model preindustrial level of 15 °C (41). To address these slight drifts, we linear-detrended all simulation results so as to bring final pCO2 and GMT simulation levels into agreement with initial (prescribed preindustrial) levels. Since our results are presented as anomalies from preindustrial/Holocene conditions, they are essentially unaffected by this procedure.