Slow dynamics of the Northern Hemisphere glaciation

2.1. Timescale

[5] A new, absolutely dated magnetobiostratigraphic timescale (Table 1) was developed to redate older time series (published before 1994) and erect an accurate time frame for the 2–4 Myr interval. The magnetostratigraphic dates exhibit an excellent agreement with astronomically tuned dates, the average error is ∼25,000 years (kyr) [Mudelsee, 2005]. Biostratigraphic dates are likely less accurate because of stratigraphic uncertainties and site-specific dependences of the habitats.

2.2. Ice Volume and Temperature Proxy Records

[6] The NHG database (Table 2) includes recently published marine δ¹⁸O time series [Lisiecki and Raymo, 2005]. Twenty-five moderate to high-resolution isotope records (spacing less than 11 kyr or one half precession cycle) and twenty lower-resolution records (≥11 kyr, or with hiatuses) provide excellent global coverage. Records from benthic foraminifera (29, thereof 18 high-resolution) document changes in global ice volume and bottom water temperature (BWT); planktonic δ¹⁸O data (16, thereof 7 high-resolution) record ice volume and sea surface temperature (SST). Only four temperature proxy time series (from faunal assemblages and Mg/Ca techniques) are available for the 2–4 Myr interval (Table 2).

3.1. Regression

[7] We employed a simple parametric model of a climate change, a “ramp” [Mudelsee, 2000], to infer the start, end, and amplitude of the increase in mean δ¹⁸O and temperature from 2–4 Myr. Superimposed on this long-term (>1 Myr) trend are fluctuations, which come mainly from short-term (≤41 kyr) variations in ice volume and temperature. Such variations are partly induced by changes in the Earth's orbital parameters obliquity and precession [Berger et al., 1984]. Salinity effects [Whitman and Berger, 1992] and measurement uncertainties add some additional “noise.” To estimate the resultant uncertainties of the NHG parameters, start, end, and amplitude, we used extensive computer simulations that take data noise into account.

[8] Fits of the ramp regression model to the δ¹⁸O and temperature time series used a weighted least squares criterion for determining the amplitude combined with a brute force search for determining start and end [Mudelsee, 2000]. The increase in variability (time-dependent standard deviation, σ(t)) across the NHG was quantified by means of a running window. This allows not only to determine the weighting for the regression, but also to measure how much the size of the δ¹⁸O fluctuations increased as the ice sheets grew.

[9] We estimated the persistence time, τ, of δ¹⁸O fluctuations by the decay period of the autocorrelation function [Mudelsee, 2002] fitted to the standardized regression residuals (a residual is given by the difference between the δ¹⁸O value and the ramp model value, divided by σ(t)). To detect outliers (positive and negative δ¹⁸O extremes), we adopted a ±3 σ(t) threshold.

[10] We performed computer simulations [Efron and Tibshirani, 1993; Mudelsee, 2000] as follows to derive parameter uncertainties: (1) Take fitted ramps, σ(t) and τ to produce a simulated time series with a random number generator. (2) Fit ramp to simulated data to obtain simulated parameters (start, end, amplitude). (3) Repeat steps 1 and 2 to produce 400 copies of simulated parameters. (4) Take standard deviation of each set of copies as uncertainty estimate.

[11] Comparing the mean of each set of copies with the original result showed that all estimations carried out had negligible bias. Because statistical uncertainties of start and end (Table 3) are clearly larger than timescale uncertainties (Table 2) the regressions do not require timescale simulations.

3.2. Calibration of δ¹⁸O Changes

[12] Calibration formulas are required to calculate temperature or salinity changes into δ¹⁸O changes, and vice versa. We used the calibration of temperature-related δ¹⁸O changes from modern and experimental data [Chen, 1994], Δδ¹⁸O_T/ΔT = −0.234‰/°C. An estimation of the statistical error (1 − σ) of this formula yields 0.003‰/°C, which we inflated in a conservative approach to a value of 0.01‰/°C (S. Mulitza, Research Center Ocean Margins, personal communication, 2004) to take “species noise” and violations of the actualism assumption into account. The calibration of salinity-related δ¹⁸O changes, Δδ¹⁸O_S/ΔT = 0.05‰/°C, assumes a linear relation between salinity and temperature changes and was developed for interpreting western equatorial Pacific δ¹⁸O amplitudes [Whitman and Berger, 1992]. Systematic errors caused by applying the calibration to other sites are likely small owing to the low signal proportion of Δδ¹⁸O_S.

4.1. NHG Start and End

[13] The results exhibit considerable variations in NHG start and end among the δ¹⁸O records (Figures 1 and 2and Table 3). The reason for this scatter is that the global ice volume signal of the NHG (δ¹⁸O increase), common to all records owing to an ocean mixing time of less than ∼2 kyr [Broecker et al., 1988, 2004], is obscured by noise, mainly of climatic origin (regional temperature and salinity fluctuations). We assume that the variations in estimated start and end cancel out when averaging the results from many records. This is a reasonable assumption, given climatological evidence for regional differences in temperature sensitivity [Ravelo et al., 2004] and the fact that our NHG δ¹⁸O database is so large. The average start and end values of the NHG ice volume increase are 3606 ± 60 kyr and 2384 ± 64 kyr, respectively, for the high-resolution records (Figure 1 and Table 3). The low-resolution records (Figure 2) yield an end date of 2510 ± 125 kyr that agrees within error bars, and a start date of 3350 ± 120 kyr that reflects minor systematic uncertainties. Sources of this error can be midterm deviations from the ramp form at around the start (see below) or underestimated influences of hiatuses in the low-resolution records.

4.2. Ice Volume Signal

[14] Extraction of the NHG ice volume signal within the interval 2384 to 3606 kyr cannot be achieved by averaging results from individual records because the noise components temperature and salinity would not cancel out but reflect an overall cooling. Evaluating the size of the ice volume amplitude alone on the basis of the δ¹⁸O results (Table 3) suggests an upper limit of around 0.4‰. This is indicated by values from records from several, low-latitude locations (925 b, 572 p, 625 p, 758 p, 806 p; see Figure 1 for notation), because a warming at all these sites across the NHG seems unlikely. More detailed information has to come from independent temperature proxy records from the same locations and water masses as the δ¹⁸O records. The problem is that these temperature records are scarce, which is the reason why we used all available records (including those with a low time resolution). The SST record (cold and warm season averaged) from DSDP 572 [Hays et al., 1989], based on ostracoda assemblages, shows a minimal cooling (Figure 3 and Table 3) during 2384–3606 kyr, equivalent to a δ¹⁸O increase of Δδ¹⁸O_T = 0.03 ± 0.12‰. Adopting the linear relation of Whitman and Berger [1992] between temperature and salinity changes, the cooling relates to a salinity-related δ¹⁸O decrease, Δδ¹⁸O_S = −0.01‰. Subtracting Δδ¹⁸O_T and Δδ¹⁸O_S from the 572 p amplitude of 0.36 ± 0.06‰ yields an NHG ice volume signal of Δδ¹⁸O_I = 0.34 ± 0.13‰. The high-resolution SST record from ODP 806 [Andersson, 1997], based on foraminifera assemblages, reveals a remarkable warming of 0.85 ± 0.17°C, equivalent to Δδ¹⁸O_T = −0.20 ± 0.04‰ and Δδ¹⁸O_S = 0.04‰, converting the low 806 p amplitude (Table 3) into an ice volume signal of Δδ¹⁸O_I = 0.43 ± 0.06‰. The BWT record from ODP 806 [Lear et al., 2003], based on Mg/Ca thermometry, reveals despite its low resolution that, while the surface waters warmed at the Ontong Java Plateau (western equatorial Pacific) across the NHG, the deep waters cooled (Δδ¹⁸O_T = 0.24 ± 0.12‰, Δδ¹⁸O_S = −0.05‰, and Δδ¹⁸O_I = 0.25 ± 0.13‰). The Mg/Ca–derived BWT record from DSDP 607 [Dwyer et al., 1995], has Δδ¹⁸O_T = 0.15 ± 0.07‰, Δδ¹⁸O_S = −0.03‰, and Δδ¹⁸O_I = 0.41 ± 0.09‰. The Mg/Ca–derived BWT decrease of 1°C between 3.7 and 2.5 Ma (not continuously measured) at site ODP 1085 (D. H. Andreasen, Rutgers University, personal communication, 2004), finally, translates into Δδ¹⁸O_I = 0.35‰. It is remarkable that the few Pliocene temperature records available so far, derived from distributed geographical sites and water masses with different proxy techniques, reveal, within error bars, the same ice volume signal of the NHG between 3606 and 2384 kyr, namely, an increase of Δδ¹⁸O_I = 0.39 ± 0.04‰ (weighted mean). The equivalent sea level lowering [Mix and Ruddiman, 1984] is 43 m, with an uncertainty of somewhat more than 5 m.

4.3. Climate Variability

[15] The time-dependent δ¹⁸O standard deviation, σ(t), measuring the size of shorter-term fluctuations in ice volume, temperature, and salinity, shows an increase across the NHG, from 3200 kyr and ∼0.2‰ to 2560 kyr and ∼0.3‰ (Figure 4). No significant differences were found between benthic and planktonic records. A qualitatively similar long-term increase in variability had previously been recognized in the benthic δ¹⁸O record from DSDP 552 [Shackleton et al., 1988]. The timing and slowness of the increase in standard deviation support the timing and slowness of the increase in mean ice volume in the Northern Hemisphere: Ice sheets have to exist before fluctuations in their size occur.

4.4. Climate Persistence

[16] The estimated persistence time, τ, of δ¹⁸O fluctuations is of the order of a few kiloyears (Table 3). Correlation analyses showed that τ does not depend on the sedimentation rate (Table 2; r = 0.06) or the temporal spacing (Table 2; r = −0.04), which means that the sampling or effects within the sediment cores such as bioturbation likely did not influence τ. Rather, the persistence is thought to be of climatic origin and reflect the time constant [Imbrie and Imbrie, 1980] of the fluctuating part of the Antarctic and the developing Arctic ice sheets. The values are slightly smaller for planktonic than for benthic records, which we explain by a stronger “dilution” of the planktonic signal by short-term atmospheric interactions of surface waters. The benthic τ estimates, averaging to 7.1 ± 0.4 kyr, should come closest to the time constant of Pliocene ice sheets. This value is significantly smaller than the 17 kyr for the larger late Pleistocene ice sheets [Imbrie and Imbrie, 1980], in agreement with recent sedimentation modeling [Lisiecki and Raymo, 2005].

4.5. NHG Dynamics

[17] Over the 2–4 Myr interval, more or less all of the analyzed records (Figures 1, 2, and 3) exhibit a long-term trend, which is well described by the ramp model. There were few short excursions (extremes) and midterm deviations from this trend. The earliest δ¹⁸O extreme of likely global significance within 2–4 Myr was at marine isotope stages (MIS) M2–MG2 (3295–3340 kyr [Lisiecki and Raymo, 2005]), detected in 11 of 25 high-resolution records as a glacial event exceeding three standard deviations (Figure 1). Stage M2 seems to reflect a stronger glaciation than stage MG2 [Lisiecki and Raymo, 2005]. The global signature and IRD evidence from several sites in the North Atlantic/Nordic Seas [Kleiven et al., 2002] indicate a significant northern ice volume proportion. A systematic, midterm deviation from the glaciation trend was the Pliocene climate optimum [Dowsett et al., 1994], documented in 12 of 25 globally distributed, high-resolution records as an interval of lowered δ¹⁸O from ∼3250 kyr to ∼3050 kyr (Figure 1). Given the duration and the size of the δ¹⁸O decrease in many records (e.g., 607 b), a significant deglaciation seems likely. As regards the dispute about the stability of the Pliocene Antarctic ice shield [Shackleton et al., 1995a], our results combined with the IRD evidence [Kleiven et al., 2002] suggest that this deglaciation could have been a Northern Hemisphere event, leaving the Antarctic shield intact. Other, more heterogeneously distributed, shorter-term cold and warm extremes accompanied the long-term NHG glaciation trend (Figure 1) until a critical threshold was passed at 2.7 Ma [Haug et al., 1999] and strong glaciations (often exceeding 3 σ in Figure 1) at MIS 96–98–100 (2433–2523 kyr [Lisiecki and Raymo, 2005]) occurred. This heterogeneity and the above finding of spatial variations in start, end, and δ¹⁸O amplitude indicate differences on short-term to midterm timescales (∼41–100 kyr) among regional temperature histories across the NHG (Table 3). This supports the idea that during a global climate change, regional subsystems with varying feedback mechanisms emerge in response [Ravelo et al., 2004].

[18] In the geographical distribution of the long-term NHG temperature change (Figure 5), we see following provisional grouping. Bottom water cooling was, within the error ranges, more or less globally homogenous. Only 6 of 29 benthic records reflect larger deviations (Figure 5a), which may be due to reduced data quality (low-resolution records) or regional influences such as changes in upwelling (e.g., at DSDP 397 site). It is unclear why the result from the benthic record at site ODP 925 (0.1°C warming) deviates significantly from the results from nearby sites ODP 927 (0.9°C cooling) and ODP 929 (0.7°C cooling). The remaining 23 records have an average BWT decrease of ∼0.7°C.

[19] SST changes (Figure 5b) are divided into three groups. First, high latitudes likely cooled strongly (704 p). Second, at lower latitudes, the surface waters in the western Pacific region warmed significantly, by ∼0.7°C on average. Third, low-latitude surface waters elsewhere had, to first order, no long-term temperature change, which supports to some degree the hypothesis of a temperature buffer [Jenkins, 1992]. However, this does not apply to the western Pacific (warming) or some other regions (e.g., Sites 851 and 722). A more detailed picture of NHG temperature trends will emerge when more high-resolution BWT and SST records for the 2–4 Myr interval become available.