1.Introduction

Spectral analysis of δ¹⁸O records and the application of other statistical techniques have shown that significant climatic transitions have occurred during the Pleistocene epoch (e.g. Reference Pestiaux and BergerPestiaux and Berger, 1984; Reference Ruddiman, Shackleton and McIntyreRuddiman and others, 1986; Reference MaaschMaasch, 1988; Reference Ruddiman and RaymoRuddiman and Raymo, 1988). Maasch (Reference Maasch1988) concluded from the statistical analysis of several oxygen-isotope records and also sea-surface temperature time series distributed worldwide that there was an abrupt global transition in the δ¹⁸O records around −900 ka. Prior to this mid-Pleistocene transition, the ice sheets were of moderate size. After that time their amplitudes increased, leading to lowered sea-surface temperatures and a global climate cooling. He showed that there was an abrupt transition in the mean; however, no firm conclusions could be made about the rapidity of the transition in the variance. A speculation as to why this might be is the existence of a gradual increase of the amplitude of the 100 ka period after the mid-Pleistocene transition, that reaches a plateau around 450 ka before present.

In this paper we present the results of ice-volume simulations for the entire Pleistocene that were obtained using a one-dimensional time-dependent ice-sheet model. The orbitally-forced one-dimensional ice-sheet model is briefly described in Section 2. The incorporation of a generalized feedback mechanism (Reference PollardPollard, 1983a, b) leads to complete deglaciations throughout the entire late Pleistocene which are in fact independent of the time of initiation of the integrations; that is to say that if the integration was started for example at −2000 or −900 ka, one would obtain complete deglaciations over the entire time spans considered. This is in disagreement with observations. A late-Pleistocene simulation of this type is described in Section 3. Multiple-window harmonic analysis (Reference ThomsonThomson, 1982) and techniques for detection of a jump in the mean and and variance have been applied to the ODP 677 record (Reference Peltier and ShackletonPeltier and Shackleton, 1989) from the Panama Basin which is of especially high resolution. This statistical information is employed to redesign the ice model so as to achieve an improved simulation using procedures described in section 4.

2. Ice-Age Model Formulation

The ice-age model includes the effects of the cryosphere, the asthenosphere, and the atmosphere. A detailed description of it has been provided in Deblonde and Peltier (Reference Peltier and Shackleton1989). For present purposes it will suffice to begin with a very brief summary of its main characteristics.

The non-Newtonian flow of the ice sheet is assumed to be governed by a power law rheology. The evolution equation for the ice thickness (H), which follows from conservation of mass, is highly non-linear and is solved in spherical coordinates for an axisymmetric spherical geometry. This evolution equation for H is simply:

in which θ is the co-latitude, r is the radius of the Earth, A is the net mass-balance function, λ = 1.825 × 10^–11 m^–3 s¹, and h is the ice height above a given reference level (e.g. Reference Birchfield, Weertman and LundeBirchfield and others, 1981).

The simulation of the bedrock depression h′ (h′ = H − h) takes into account the vertical variations of density and viscosity of the Earthȁs interior. The isostatic adjustment theory developed by Peltier (e.g. 1982) shows that h′ is given by the following convolution integral:

in which ϕ is the longitude and u_r is the Green function for radial displacement. Since this model is axisymmetric, Equation (2) with the addition of initial topography assumed to be in isostatic balance reduces to:

with the Love numbers q_l given as follows:

ρ _a is the average density of the Earth (5515 kg/m³), ρ_i is the ice density, and

is the instantaneous elastic contribution to the response, is the amplitude of the mode of relaxation j at degree l and the are the inverse relaxation time constants. The relaxation strength is measured by the ratio of over The five strongest modes of relaxation are kept in the system (Reference PeltierPeltier, 1982). A lithospheric thickness of 120 km and a lower mantle viscosity of 2 x 10²¹ Pa s are also assumed. The net mass-balance function is prescribed (Reference OerlemansOerlemans, 1980) as follows:

where a = 0.81 × 10^–3 s^–1 and b = 0.3 × 10^–6 m^–1 s^–1. E is the equilibrium line as is defined as follows:

S is the slope of the equilibrium line and θ₀ is a fixed co-latitude. For the simulations to be described based upon this model, we shall assume S = 0.1 × 10^–2 and θ₀ = 20°. E₀ and k are given constants. δQ_S is the caloric summer insolation anomaly defined according to Milankovitch and has been computed after Berger (1978a, b).

3. Late-Pleistocene Ice-Volume Simulations (−700 to 0 ka)

The major ice retreats that are obtained with both a locally-damped and a thin-channel asthenospheric bedrock depression model are no longer so prominent when the bedrock despression is modelled by the viscoelastic model described above (Reference Peltier and ShackletonDeblonde and Peltier, 1989). These differences are due to (1) a different dependence of the relaxation time on wavelength and (2) the variation of the relaxation strength with wavelength in the viscoelastic model. However, the addition of a generalized melt-water parameterization mechanism which depends on the time rate of change of the ice volume itself (Reference PollardPollard, 1983a, b) leads to complete deglaciations but not to a maximum in harmonic amplitude at the 100 kyear period as was the case for the thin asthenospheric channel model (Reference PollardPollard, 1982, Reference Pollard1982a, b). A time series of the ice volume for two different strengths of the feedback mechanism are shown in Figure 1. The linear correlation coefficient of the synthetic ice-volume record with SPECMAP (Reference ImbrieImbrie and others, 1984) is 0.57 (Fig. 1; solid line).

4. Pleistocene Ice-Volume Simulations (−1900 to 0 ka)

The same statistical techniques for the detection of the jump in the mean and variance as employed in Maasch (Reference Maasch1988) have been applied to the planktonic B¹⁸O record ODP 677. The δ¹⁸O series in shown in Figure 2a. For this series, it is quite obvious from visual inspection that it has distinct properties before and after about −900 ka. Figure 2b is a time series of the running mean with a segment length of 300 ka. Figure 2c is a curve of the changing standard deviation with the same segment length. The signal-to-noise ratio for the detection of a jump in the mean is shown in Figure 2d. Two distinct peaks (that must be greater than 1.0 to be significant at the 99% level) are observed around −950 and 600 ka. The signal-to-noise ratio has been computed for segment lengths in the range 300-450k; and then averaged. Figure 2e illustrates the Mann-Kendall rank statistic (Reference SneyersSneyers, 1975) for the forward and retrograde time series. An intersection of the curves occurs within the 95% confidence level around −300 ka. Thus, the ODP 677 time series has a jump in the mean around −900 ka that is statistically significant. Figure 2f illustrates values computed Trom the 99% confidence intervals of the variance. If D ⁺ > 0 and D⁻ > 0, then there is an increase in variance. These values have also been averaged for a range of 300-450k: to ensure that the results are independent of segment length. This last test does not point to a sharp transition in variance but rather hints toward a gradual increase of the variance and hence a possible gradual increase of the amplitude of the 100 ka cycle. The ratio of the standard deviation after and before −900 ka is 1.75 and this is statistically significant (F-test applied to a ratio of variances).

Fig. 1. Total ice cross-sectional area for the ice model with orbital forcing at 55°N. The increases in height of the equilibrium line used in the GMP mechanism are at 1000 m (solid line) and 500 m (dahsed line). The function describing the increase in the height of the equilibrium line (E+) is as follows: it is zero until a critical melt-water rate of M_l is reached, the E⁺ increases linearly until the melt-water rate (M) reaches a value of M₂ (such that M₂ > M_t) where E* is set equal to a chosen value E+ ₂ for higher values of the melt-water rate, E⁺ is kept fixed at the value of E⁺ ₂. For all cases considered in this paper, M₁ = 12 cm year^–1 and M₂ = 20 cm year^–1 unless specified otherwise.E⁺ ₂ varies between 500 and 1000 m.

Fig. 2. Statistical analysis of site ODP 677 δ¹⁸O record, (a) Normalized time series, (b) Running mean of the time series shown in (a) with segment length of 300 ka. (c) Changing standard deviation of the time series shown in (a) computed for a segment length of 300 ka. (d) Signal-over-noise ratio averaged over a range of segment lengths of 300-450 ka. (The confidence level is 99%.) (e) Mann-Kendall rank statistic, U(Di), for forward (solid) and retrograde (dashed) time series, (f) Test for jump in variance, D⁺ (solid line) and D⁻ (dashed line), averaged over the same segment lengths as in (d).

This statistical input has been used to develop a simulation of ice volume over the entire range of Pleistocene time. The one-dimensional ice model is forced externally by the orbital variations. If one assumes that orbital reconstructions are correct in their prediction that the nature of the forcing has not changed significantly within the Pleistocene time then a mid-Pleistocene transition is clearly not obtainable with such a model. However, it is nevertheless possible to evaluate how well such a model would be able to reproduce this transition given further information as input to it. Additional degrees of freedom may be added to the system by varying the height of the equilibrium line at θ₀(E₀ in Equation (4)), the sensitivity to the orbital forcing (k in Equation (4)) and the strength of the generalized melt-water parameterization (GMP). E_Q is one of the variables that controls the maximum ice volume. Raising the equilibrium line generates a warmer climate and therefore smaller overall maxima in ice mass. k controls the standard deviation of the ice volume. Using a lower value of k results in a lower value for the standard deviation of the ice volume. The question as to how fast E₀ and k might change can be assessed by comparing model output with the δ¹⁸O time series ODP 677.

The gradual transition of the 100 ka amplitude is simulated by changing the strength of the sudden increase in equilbrium-line height that occurs when a critical value of ice volume time rate of change is obtained which sets off a wastage process. It is assumed that the GMP mechanism is absent before -900 ka and not as strong during a 400 ka transition period. The linear correlation coefficients for the cases discussed in this paper are listed in Table I. As expected, using a constant GMP during the entire Pleistocene leads to several unobserved complete deglaciations (case 3). Starting the GMP mechanism at −800 ka and gradually increasing its strength to −400 ka and then keeping it constant gives rise to a series that has a relatively close fit with the ODP 677 series. Starting the GMP mechanism at −900 ka instead of −800 ka leads to very similar results. The linear correlation coefficient of the simulated time series with SPECMAP for the time span −700 to Oka increased from 0.57 to 0.67 (case 1). The same coefficient for ODP 677 is 0.62. One should note that the correlation coefficient between SPECMAP and ODP 677 over the last 700 ka is 0.81, which rules out the expectation of obtaining a correlation coefficient of 1 between the model and individual observed time series. The values of k and E₀ giving rise to the highest correlation coefficients as well as to a similar ratio in standard deviation of the ice volume are 24 and 35 m (W m²)^–1, respectively and 700-550 m, respectively, before and after a transition period that is centred around −900 ka. The change in E₀ corresponds to a change of about 1°C using the current value of the atmospheric lapse rate (−6.5°Ckm^–1). The change in k of about 30% suggests the existence of a less

Fig. 3. Modelled ice-volume time series for case 1 (Table I) (solid line) and ODP 677 δ¹⁸O time series (dashed line).

sens itive climate system before −900 ka . As can be seen from the linear correlation coefficients (cases 1 and 2), it is not possible to deduce in a convincing manner that the transition in E₀ and k was rapid (i.e. less than 200 ka). The ratio of the standard deviation of case 1 (Fig. 3) before and after −900 ka is 2.1 . A harmonic analysis for the time span from −900 to 0ka for both case I and ODP 677 (Figs 4 and 5) shows that the modelled 100 ka amplitude is inadequate. However, if the harmonic analysis is confined to the last 400 ka, then the amplitude for both cases is a maximum at the 100 ka period. The harmonic analysis of the modelled time series before -900 ka (case 1 ) shows similarities in relative amplitude at the 41, 23, and 19 ka cycles provided the orbital forcing is taken at 55°N (Figs 6 and 7). The modelled maximum southern latitudinal extent of the ice sheets before −900 ka is 3° latitude further north. Using an orbital forcing at 65°N before −900 ka (case 4) since the ice sheets did not extend as far south does increase the amplitude of the 41 ka cycle and hence does not give rise to as good a fit with the observations.

Fig. 4. Multiple window harmonic analysis of the ice volume for case 1 for the time series covering −900 to 0 ka. The top figure shows the harmonic amplitudes. The bottom figure shows the F- variance ratio of 2 and 8 degrees of freedom. Peaks above the horizontal lines are significant above the 97.5% point. Five windows have been used and NW = 8 with N = 901. A linear regression line has been subtracted from the time series.

Fig. 5. Same as Figure 4 but for the ODP 677 time series covering −900 to 0 ka.

Fig. 6. Same as Figure 4 but for the modelled time series of ice-volume case 1 covering −1990 to −900 ka and N = 1091.

Fig. 7. Same as Figure 4 but for the ODP 677 time series covering −1990 to −900 ka and N 1091.

5. Conclusions

By introducing further degrees of freedom in the one-dimensional ice-sheet model that includes a GMP feedback mechanism, it is possible to simulate the general features of the mid-Pleistocene climate transition and the gradual increase of the amplitude of the 100 ka cycle after this climatic change. The reason for the rapidity of the mid-Pleistocene transition (probably > 200 ka) is still an open question.

According to the model, changes in E ₀ and k around −900 ka of 250 m and 11 m(W m^–2)^–1, respectively, and the introduction of the GMP mechanism around −900 ka are necessary to simulate the observed deep sea-core record changes. The changes to E₀ and k are necessary to simulate the reduction in ice-sheet volume before about −900 ka, and the changes in the strength of the GMP mechanism are ncessary to simulate the appearance and increase in 100 ka amplitude after the mid-Pleistocene transition. Possible physical mechanisms causing these changes include tectonic uplift (Reference Ruddiman and RaymoRuddiman and Raymo, 1988) and variations in the oceanic thermohaline circulation (e.g. Reference SachsSachs, 1976; Reference Ruddiman and McIntyreRuddiman and Mclntyre, 1981; Reference Manabe and BryanManabe and Bryan, 1985; Reference Boyle and KeigwinBoyle and Keigwin, 1987). Cyclic variations in carbon dioxide may be important for individual 100 ka cycles (Reference Shackleton and PisiasShackleton and Pisias, 1985; Reference Barnola, Raynaud, Korotkevich and LoriusBarnola and others, 1987). It is also conceivable that the long-term variations (≫100 ka) of carbon dioxide could also have changed over the Pleistocene epoch.