Integrated Summer Insolation Forcing and 40,000-Year Glacial Cycles: The Perspective from an Ice-Sheet/Energy-Balance Model

Although the origins of the 40,000-year glacial cycles during the early Pleistocene are readily attributed to changes in Earth’s obliquity (also having a 40,000-year period), the lack of ice-volume variability at precession periods (20,000 years) is difficult to reconcile with most parameterizations of the insolation forcing. It was recently proposed that precession’s influence on glaciation is muted because variations in the intensity of summer insolation are counterbalanced by changes in the duration of the summertime, but no climate model has yet been shown to generate obliquity period glacial cycles in response to the seasonal insolation forcing. Here we present a coupled ice-sheet/energy-balance model that reproduces the seasonal cycle and, when run over long time periods, generates glacial variability in response to changes in Earth’s orbital configuration. The model is forced by the full seasonal cycle in insolation, and its response can be understood within the context of the integrated summer insolation forcing. The simple fact that obliquity’s period is roughly twice as long as that of precession results in a larger amplitude glacial response to obliquity. However, for the model to generate almost exclusively obliquity period glacial variability, two other conditions must be met. First, the ice sheet’s ablation zone must reside poleward of ! 60 ! N because insolation intensity is more sensitive to changes in Earth’s obliquity at high latitudes. Second, the ablation season must be long enough for precession’s opposing influences on summer and fall insolation intensity to counterbalance one another. These conditions are consistent with a warm climate and a thin ice sheet, where the latter is simulated as a response to subglacial sediment deformation. If a colder climate is prescribed, or in the absence of basal motion, ice sheets tend to be larger and undergo greater precession period variability, in keeping with proxy observations of late Pleistocene glaciation. and (2008), Integrated summer insolation forcing and 40,000-year glacial cycles: The perspective from an ice-sheet/energy-balance model, Paleoceanography , 23 , PA1208, doi:10.1029/2007PA001463.


Introduction
[2] That the $40,000-year glacial and temperature variability during the early Pleistocene [e.g., Shackleton and Opdyke, 1976;Pisias and Moore, 1981;Shackleton and Hall, 1984;Raymo et al., 1989;Ruddiman et al., 1989;Imbrie et al., 1993;Tiedemann et al., 1994;Venz et al., 1999;Liu and Herbert, 2004;Huybers, 2007;Lawrence et al., 2006] are in some way related to the 40,000-year changes in Earth's obliquity seems almost certain, but such a description is incomplete. The caloric summer half year [e.g., Milankovitch, 1941] and the insolation intensity on any summer day [e.g., Hays et al., 1976;Imbrie et al., 1992] varies primarily with the precession of the equinoxes, or at $20,000-year periods (assuming that the eccentricity of the Earth's orbit is not anomalously small). Why, then, is there not more precession period variability during the early Pleistocene?
[3] Kukla [1968] as well as Muller and MacDonald [2000] speculate that the driver of glaciation may be winter insolation at 65°N, sensitive primarily to obliquity because it dictates whether 65°N experiences polar night. However, while this yields an insolation quantity having the correct timescale, it is difficult to rationalize how winter insolation variability could affect glaciation, both because no ablation is expected in the winter and the amplitude of the variability is much smaller than for summer. Another quantity sensitive primarily to obliquity is the meridional insolation gradient [Berger, 1978;Young and Bradley, 1984;Johnson, 1991;Raymo and Nisancioglu, 2003]. Increasing obliquity causes an annual average redistribution of insolation from latitudes equatorward of 44°N to latitudes poleward of this hinge point, whereas precession causes insolation to change more uniformly with latitude. Raymo and Nisancioglu [2003] proposed that insolation gradients would modulate the poleward transport of moisture and heat and thereby determine glacial mass balance, but model simulations suggest that the obliquity-induced variations in the insolation gradient have a minor influence on glaciation [Jackson and Broccoli, 2003;Nisancioglu, 2004]. Another suggestion [Ashkenazy and Tziperman, 2004] is that early Pleistocene glacial cycles are nonlinear oscillations with an intrinsic period close to 40 Ka, and that these become phase-locked to the obliquity forcing [Gildor and Tziperman, 2000;Tziperman et al., 2006].
[4] A more recent hypothesis calls upon how glacial variability is recorded to mute the precession signal. Raymo et al. [2006] suggest that precession period variations in glaciation are antisymmetric between the hemispheres, but whose effects are averaged in the oceans and thus absent from records of marine calcite d 18 O. To explain why precession variability is not present in ice-rafted debris records from the North Atlantic [e.g., Shackleton and Hall, 1984;Ruddiman et al., 1989] and Nordic Seas [e.g., Jansen et al., 2000], which covary with d 18 O at predominantly 40-Ka periods, Raymo et al. [2006] call upon delivery of ice-rafted debris to occur only during obliquity induced excursions in sea level and not during precession induced ablation events. Note, however, that North Atlantic sea surface , intermediate [McIntyre et al., 1999], and deep water temperatures [Dwyer et al., 1995] all appear to vary at 40-Ka periods. Also, ice volume appears to have varied at 40,000-year periods prior to the glaciation of Antarctica [Zachos et al., 2001], indicating that at some point in Earth's history a mechanism existed for causing Northern Hemisphere ice sheets to vary primarily at obliquity periods.
[5] Last, it has been suggested by one of the authors [Huybers, 2006] that the absence of precession period variability during the early Pleistocene occurs because the Earth's distance from the sun is inversely proportional to its angular velocity, i.e., Kepler's second law. Annual ablation is expected to be a function of both the solar radiation intensity (related to Earth's distance from the sun) and the duration of the ablation season (related to Earth's angular velocity). These two influences are both taken into account by integrating insolation intensity over the time variable duration of the summer, a quantity referred to as the summer energy. In this case, summer is defined as the period during which daily average insolation intensity exceeds a specified threshold. The summer energy is generally insensitive to precession because, for example, just when perihelion occurs at summer solstice, summertime is shortest. That is, Kepler's second law dictates that duration and intensity counterbalance one another. Here we elaborate on what controls the glacial sensitivity to changes in obliquity and precession using a coupled ice-sheet/energy-balance model.
[6] Numerous studies have used numerical models to explore glacial variability, but for such a model to address the summer-energy hypothesis requires the representation of both the full seasonal cycle and the secular variations in precession, obliquity, and eccentricity, thus spanning 5 orders of magnitude in timescales. We are aware of only four previous studies which include the full seasonal cycles and attempt to simulate 40-Ka glacial cycles. Berger and Loutre [1997] developed a model having zonally averaged ocean and land components coupled to a parabolic ice sheet. It generated glacial variability of roughly equal magnitude at obliquity and precession periods during the early Pleistocene as well as considerable energy spread across a continuum of timescales. Nisancioglu [2004] coupled a zonally averaged atmospheric energy balance model, including a representation of moisture and heat transport, to a parabolic ice sheet model and obtained somewhat greater precession than obliquity period variability. Jackson and Broccoli [2003] coupled an atmospheric GCM to a slab ocean and ran it to equilibrium under various orbital configurations. As their model did not include an ice sheet, they analyzed variations in precipitation and positive degree days and, again, find somewhat more precession than obliquity period variability. Finally, DeConto and Pollard [2003aPollard [ , 2003b employed a coupled atmospheric GCM and ice sheet model to represent Antarctica. While their focus was on initiation of glaciation in Antarctica, the orbital period fluctuations in ice volume appear (from visual inspection of their figures) to be almost exclusively at obliquity periods. It should be noted, however, that to speed up the model integration, idealized variations in orbital parameters were used wherein obliquity and eccentricity variations were approximated to have periods which are integer multiples of the precession cycle, making it difficult to interpret sensitivities to the individual orbital parameters.
[7] Raymo et al. [2006] point out that no climate model has simulated the 40,000-year glacial cycles as a response to the full insolation forcing: a fair challenge to any hypothesis seeking to explain the record of early Pleistocene glacial variability as directly indicating variations in ice volume.
Here we present a coupled ice-sheet/energy-balance model which calculates ablation as a function of the surface energy balance and is driven by the full seasonal cycles of insolation, sensitive to both precession and obliquity. The model generates 40-Ka glacial cycles consistent with the early Pleistocene observations and the predictions of the summer energy. We first describe some features of the summer energy and then discuss results from the model.

Summer Energy
[8] The concept of positive degree days (PDDs) provides a simple relationship between temperature and glacial ablation [Braithwaite and Zhang, 2000], Here T d is the daily temperature and d d is zero if T d < 0°C and one otherwise. The summer energy is defined in close analogy with PDDs, as the sum of the insolation intensity exceeding a threshold, where F d is the daily average insolation intensity, and b d is one when F d is above a threshold, t; otherwise, b d is zero. Note that the full magnitude of the insolation intensity, not only the portion above the threshold, is summed under the assumption that most incident radiative energy will eventually lead to ablation once the freezing point is obtained. In order to provide a conceptual framework by which to interpret the model results, which are presented in subsequent sections, it is useful to first explore how summer energy depends on latitude and the threshold. Huybers [2006] only discussed summer energy at 65°N with a threshold of 275 W/m 2 .
[9] The variance of the summer energy, computed between 2 and 1 Ma, ranges from zero to (3000 Giga-Joules/ m 2 ) 2 , depending on the threshold and latitude (Figure 1a). Variance in J is greatest when the threshold is just below the maximum insolation received at each latitude, partly because insolation intensity is large but also because the time rate of change in insolation is small; shifts in insolation intensity then cause relatively large changes in the number of days exceeding the threshold. In the tropics and midlatitudes there is negligible summer-energy variability at low thresholds because intensity exceeds t all year, and only changes in the annual average intensity influence J. Nearly all the summer-energy variance is partitioned between the obliquity and precession bands (Figures 1b and 1c). Obliquity dominates at most threshold values below about 200 W/m 2 because insolation is integrated over nearly the whole of the year and precession has no influence on annually integrated insolation [Rubincam, 1994]. An exception is at 44°N and nearby latitudes, where shifts in Earth's obliquity have negligible influence on annual insolation, leaving only the small contribution from changes in eccentricity. At latitudes above 60°N, presumably the most relevant for determining the extent of glaciation, J is primarily a function of obliquity up to thresholds of 350 W/m 2 .
[10] The switch from obliquity to precession period variability at thresholds above 350 W/m 2 can be understood by considering the change in the structure of the annual cycle resulting from the precession of the equinoxes. When perihelion (i.e., Earth's closest approach to the Sun) occurs in Northern Hemisphere summer, the insolation intensity increases between May and August but decreases between September to November ( Figure 2). This fall deficit owes to the fact that summer perihelion is associated with Earth being unusually far along its orbital path. Low threshold values tend to intersect the insolation intensity curve in the fall, and therefore summer perihelion causes a reduction in the number of days during which insolation intensity exceeds t. Thus the increased intensity associated with summer perihelion is counterbalanced by a decrease in exposure time, giving a summer energy nearly independent Figure 1. Summer-energy variability as a function of latitude and threshold. (a) Variance of the summer energy contoured in units of log 10 (Giga-Joules/m 2 ) 2 . (b) Fraction of summer-energy variance occurring at the obliquity band (1/41 ± 1/150 Ka À1 ) and (c) fraction of variance at the precession band (1/21 ± 1/ 150 Ka À1 ). Together, the obliquity and precession bands account for over 90% of the summer-energy variance, with the remainder associated with eccentricity and the first overtone of precession, i.e., $11 Ka. Only the Northern Hemisphere is shown; values are symmetric for the Southern Hemisphere. Representative time series of summer energy, as well as software to calculate insolation and summer energy, can be downloaded at http://www.ncdc.noaa.gov/paleo/pubs/huybers2006b. For the lower threshold, the summer energy is insensitive to changes in precession (the intensity is counterbalanced by changes in duration), while the higher threshold gives a summer energy sensitive primarily to precession. of precession. Higher threshold values, however, intersect the insolation intensity curve between July and August, in which case summer perihelion causes an increase in insolation intensity as well as duration and a corresponding increase in summer energy. Note that, following convention, the year is defined with respect to Northern Hemisphere Spring Equinox.
[11] There exists a conspicuous valley bisecting the precession period variability in Figure 1c. The lobe of precession variability associated with higher thresholds is intensity controlled; that is, maxima in summer energy are associated with perihelion occurring during summer. Conversely, the lobe of precession variability at low thresholds is durationally controlled; that is, summer energy is largest when aphelion (Earth's furthest excursion from the Sun) occurs during summer. The intervening valley is where summer-energy maxima occur for both summer aphelion and summer perihelion, giving a period of half a precession cycle, or $11 Ka.
[12] The latitude-threshold space in which an ice sheet's ablation zone exists helps determine the ratio of obliquity to precession period variability. When ablation occurs at thresholds below $350 W/m 2 and at latitudes above 60°N, glaciers are expected to be most sensitive to changes in Earth's obliquity. Moving the ablation zone equatorward or requiring a higher ablation threshold would increase the precession period variability. These inference are further explored in the following sections using a coupled ice-sheet/ energy-balance model.

Model Description
[13] The model formulation is similar to that presented by Pollard [1978], asynchronously coupling an energy-balance model (EBM) and an ice sheet (see Figure A1 in Appendix A1). The EBM is zonally averaged and spans the equator to the pole at one degree resolution. Atmospheric heat transport is parameterized using the formulation of Stone [1972]. Ablation is computed as a function of surface energy balance, as opposed to the positive degree day formulation, which would largely presuppose the validity of the summerenergy concept. The hydrological cycle is treated very simply, ignoring latent heat fluxes, and assuming that a constant amount of ice accumulates each day for which the atmospheric surface temperature is below freezing; otherwise precipitation falls as rain and is assumed to run off.
The EBM state appears to be an unique function of orbital configuration and ice sheet topography and converges to a stable seasonal cycle within a few years. The annual ice mass balance is computed at each grid box from the equilibrated EBM and then used to force an ice sheet.
[14] The ice-sheet model is also zonally averaged and the domain spans 30°N to 85°N, outside of which calving is assumed to impose zero ice-thickness boundary conditions. Ice deformation is calculated from horizontal stress according to Glen's flow law. Early Pleistocene ice sheets may have been as spatially expansive as during the late Pleistocene [Boellstorff, 1978;Fisher et al., 1985;Joyce et al., 1993] (but see Barendregt and Irving [1998]), while ice-volume variations appear to have been smaller [e.g., Pillans et al., 1998] and as indicated by d 18 O, Figures 3c and 3e). A possible resolution is that a deformable bed underlaid the Laurentide during the early Pleistocene and led to faster flowing, more expansive, and thinner ice sheets [Fisher et al., 1985;Beget, 1986;Alley, 1991;Clark and Pollard, 1998]. In this study we assume the ice sheet rests atop a layer of deformable sediment, represented using the model of Jenson et al. [1996]. The presence of deformable sediment is, however, not necessarily required as a thin ice sheet could also result from basal sliding in response to the presence of melt water [Jenson et al., 1996].
[15] A detailed description of the energy balance and ice sheet components of the model, along with a list of the specified parameters, is provided in Appendices A1 and A2. Parameter values are selected to give reasonable agreement with the modern climate system. At certain points it will be useful to demonstrate model behavior under different parameterizations, involving changes in surface temperature and the deformability of sediment; all these deviations from the ''standard'' parameterizations are explicitly stated.

The 40-Ka Glacial Variability
[16] While the model is seasonally ice-free under modern orbital conditions, orbitally induced shifts in the insolation at the top of the atmosphere are sufficient to initiate glaciation. Eighty-two percent of the variance in ice volume is within the obliquity band (1/41 ± 1/150 Ka À1 ) and only 12% at the precession band (1/21.1 ± 1/150 Ka À1 ) (see Figure 3). The remainder of the ice-volume variance is spread throughout a continuum of timescales. The partition of variance in the composite d 18 O records of Huybers [2007] for the same time period is 54% of the variance in the obliquity band and 2% in the precession band. A similar partition is found in the composite d 18 O record of Lisiecki and Raymo [2007]. The lower fractions of variance in the d 18 O orbital bands relative to the model orbital bands is anticipated from some combination of age-model error [e.g., Huybers and Wunsch, 2004], observational noise, the fact that Earth's climate variability comprises myriad stochastic phenomena not included in the model, that d 18 O is also sensitive to changes in temperature, and that approximately 80-Ka-long glacial cycles are found near 1.6 and 1.2 Ma in the d 18 O record, associated with a concentration of variability near 1/80 Ka. The bottom line is that the model results are consistent with the observation of negligible precession period variability during the early Pleistocene ( Figure 3c). There are more nuanced discrepancies between the d 18 O record and the model results, and these are addressed in more detail in sections 4 and 5.
[17] The size of the model ice sheet is in rough agreement with estimates for 50-m amplitude sea level variations during the early Pleistocene [e.g., Pillans et al., 1998]. At its maximum extent the model ice sheet spans 60°N to 85°N with an average thickness of 1500 m ( Figure 4). Although longitude is not represented, if we assume the ice sheet spanned North America, $3000 km at high latitudes, the ice sheet's volume corresponds to $40 m of sea level. During maximum glacial extent the annual average temperatures at high latitudes decrease by approximately 20°C from an interglacial to glacial state, primarily owing to changes in albedo and surface height. Decreased temperatures lead to greater accumulation as more of the precipitation falls as snow, and is balanced by increased ablation at more equatorward and warmer latitudes. Ablation occurs almost exclusively at the leading edge of the ice sheet as this region receives relatively more insolation, is associated with the largest flux of heat from the atmosphere, and is at a low elevation.
[18] Glaciation is controlled in the model by orbital variations in two senses. First, model runs starting from different initial conditions all converge within a few glacial cycles. Second, in the absence of orbital or stochastic variations the model would reach a steady state with either a permanent ice sheet or seasonally ice-free conditions. In the presence of substantial noise and orbital variations, the model produces fluctuations in ice volume similar to when it is forced only by orbital variations. For example, when a noise term is introduced into the ice mass balance having a standard deviation of 2 m (twice the precipitation rate), the model continues to undergo 40-Ka glacial cycles. In this case, the noise is uncorrelated at the 1°grid spacing and 5year time step of the ice-sheet model. The ice sheet apparently acts as a very effective smoother in space (redistributing ice-mass anomalies) and time (being largely insensitive to high-frequency anomalies) so that the lowfrequency systematic changes induced by insolation still dominate the result. This damping of high-frequency variability may explain why the early Pleistocene glacial cycles appear to occur so regularly at 40-Ka intervals, even in the presence of the inevitable stochastic variability.
[19] This description of the model behavior still leaves the question of why the model generates primarily 40-Ka glacial variability? As a partial answer, consider the insolation threshold at which ablation is induced in the model.  shows an estimate of this threshold, obtained by sampling the minimum insolation at which ablation occurs as a function of latitude and elevation. Ablation always occurs at latitudes above 60°N and initiates at insolation intensities below 320 W/m 2 . These are just the conditions under which summer energy was shown to vary primarily at the obliquity period (see section 2). Thus, in so much as summer energy and the model behave similarly, obliquity period variability is expected.

Sensitivity Experiments
[20] Three sensitivity studies are used to illustrate the differences between how the model responds to obliquity and precession. Negative values occur where ice has retreated from an isostatically depressed bed. Also shown are (d) the annual average temperature,°C, (e) accumulation, meters per year, (f) ablation, meters per year, and (g) total heat flux into the surface including short wave, long wave, and sensible heating, W/m 2 . Vertical lines indicate maxima in obliquity.

Thin Static Ice Sheet
[21] In the first sensitivity experiment a 10-m-thick ice field is specified to extend from 60°N to 85°N. Increasing obliquity from a minimum (22.1°) to maximum (24.5°) results in an increase in insolation of 25 W/m 2 at high latitudes during the summer and a decrease in mass balance of À0.43 m/a, when spatially averaged over the ice field (Figures 6a and 6b). The change in mass balance comes from increased ablation (À0.41 m/a) and a small decrease in accumulation (À0.02 m/a). (Changes in accumulation would probably play a more prominent role if precipitation rates were made to depend on the climate [e.g., Gildor and Tziperman, 2000].) When perihelion is shifted from winter to summer solstice and the long-term eccentricity average of 0.0275 is specified, insolation intensity increases in the spring and early summer by as much as 80 W/m 2 but decreases by a corresponding amount in the later summer and fall (as also described in section 2). Likewise, mass balance decreases in the early part of the ablation season but increases in the later part of the ablation season, leading to a low sensitivity of À0.18 m/a, or less than half that associated with obliquity (Figures 6d and 6e).

Thick Static Ice Sheet
[22] In the second sensitivity study, the ice-sheet elevation is increased to 1500 m, resulting in a temperature decrease of 10°C at the ice-sheet surface. The sensitivity to obliquity decreases by a factor of 2 to À0.22 m/a, but the sensitivity to precession only decreases by roughly 20%, to À0.14 m/a (Figures 6c and 6f). The sensitivity to precession is less influenced by the colder temperatures because no ablation occurs during the fall and the counterbalancing effect is lost. The change in sensitivity between the thin and thick ice sheet runs parallels the changes anticipated in the summer energy between specifying a low and high threshold value.

Dynamic Ice Sheet
[23] An alternative method to compute sensitivity is to average the mass balance associated with minima in obliquity and subtract this from the average maximum-obliquity mass balance. Averages are taken from a model run spanning 2 to 1 Ma. The mass-balance sensitivity to obliquity, averaged between 65°N and 85°N, is À0.41 m/a, and for precession the sensitivity is À0.16 m/a. The ratio between obliquity and precession variance is obtained by squaring the ratio of the sensitivities(À0.41/À0.16) 2 = 6.6, in good agreement with the spectral distribution of energy obtained from the model run (0.82/0.12 = 6.8, Figure 3b).
[24] The sensitivity to obliquity is larger in the dynamic than the static ice sheet runs. This suggests that the time dependence of the model also plays a role in determining the ratio of obliquity to precession related variance, a topic taken up in more detail in the next section.

Consequences of a Nonequilibrium Response to Astronomical Forcing
[25] A coherence analysis shows that model ice-volume variations are almost perfectly coherent with the changes in obliquity as well as precession, lagging these orbital variations by approximately 90°. This is consistent with observational studies of the marine d 18 O record as well as dynamical expectations [Imbrie and Imbrie, 1980;Huybers, 2006;Roe, 2006], and suggests that insolation does not control the magnitude of ice volume, but instead its rate of change. A simplification of this relationship is The two terms on the right-hand side represent the influence of obliquity and precession, having frequencies of w obl = 1/41 Ka À1 and w prec = 1/22 Ka À1 , respectively. The difference in phasing as well as the frequency and amplitude modulations of the orbital forcings have been ignored; their inclusion would not alter the point being made here. Integrating gives so that ice volume lags behind the cosine forcing terms by 90°, consistent with observations. The amplitude of the ice-volume variability depends on the frequency of the forcing, and ice volume is almost twice as sensitive to obliquity as to precession, w prec /w obl = 1.9.
[26] Is this simple representation of the ice sheet's orbital response reasonable? The spatial and temporal structure of the obliquity and precession insolation anomalies are distinct. Obliquity induced changes in insolation contain Figure 5. Insolation threshold for ablation to occur in the model. The threshold is diagnosed as the minimum diurnal average insolation intensity for which ablation occurs and is contoured as a function of latitude and elevation of the ice sheet. Absence of contour lines indicates that no ablation occurs in this region. Note that the model insolation thresholds are within the domain for which summer energy is obliquity dominated, i.e., t < 350 W/m 2 and a latitude above 60°N (see Figure 1b). low-frequency components (e.g., 1/41 Ka À1 ), whereas precession variations only alter the structure of the seasonal cycle [Rubincam, 1994]. To more fully explore the influence of different periods of forcing upon amplitude of the ice-volume response, a series of model runs are performed using idealized, sinusoidal variations in orbital parameters ranging in period from 10 Ka to 100 Ka ( Figure 7). As anticipated from equation (4), the amplitude of ice-volume variations is found to generally be proportional to the period, although the relationship breaks down at periods longer than 50 Ka. The breakdown seems to arise because of the physical limits in the ice-sheet extent. For example, under constant low obliquity (corresponding to an infinite period) the ice sheet eventually grows to a steady state limit reaching to 50°N. In a somewhat cooler climate ice-sheet volume would be larger, on average, but the scaling with respect to forcing frequency behaves similarly. Thus both the model observations indicate that the amplitude of the orbital response depends on the forcing period.
[27] The predominance of the obliquity period variability can apparently be explained as a result of both counterbalancing of the precession effect and a larger obliquity response because of its lower frequency. It is possible to quantify the relative importance of these mechanisms. The obliquity-induced changes in insolation intensity at high latitudes are roughly a third those caused by precession. However, for the case of a static ice sheet, model mass balance was found to be more sensitive to obliquity, suggesting that the influence of precession is reduced by about a factor of 4 by seasonal counterbalancing. Obliquity also has a period roughly twice as long as precession suggesting another factor of 2. Multiplying the relative influences and squaring gives a ratio between the expected variance of obliquity and precession, (1/3 Â 4 Â 2) 2 % 7, , and 6c, respectively, but for changing the location of perihelion from summer to winter solstice. All calculations are with an eccentricity of 0.0275 and, unless otherwise stated, an obliquity of 23.34°and perihelion occurring at vernal equinox. Note that for the thin ice sheet the mass-balance anomalies associated with precession during summer and fall are opposite and of nearly equal magnitude, whereas for the thick and cooler ice sheet the high-latitude anomaly is more uniform. consistent with the 7:1 partition of variance between obliquity and precession found in the model run.

Effect of a Cooling Climate and Implications for the Late Pleistocene
[28] The climatic response to orbital variations itself depends on the background state of the climate. Cooling the climate is expected to be analogous to increasing the threshold for ablation in the summer energy parameterization and thus increase the relative fraction of precession period variability. To check this inference, a cooler surface temperature is induced in the model by decreasing the emission level for outgoing longwave radiation. This decreases the effective thickness of the atmosphere with respect to longwave radiation and is similar to decreasing greenhouse gas concentrations [see, e.g., Goody, 1995]. The average surface temperature in the model decreases by $1°C for each 100 m the emission level is lowered, both because of increased transmission of longwave radiation to space and because of feedbacks involving an increase in the area and thickness of the ice sheet. Precession variability increases monotonically with cooling surface conditions. For example, when the emission level is lowered from 7 km to 6.5 km, the mean temperature over the ice sheet cools by 5°C (partially suppressing ablation during the late summer and fall), and 50% of the ice volume variance is concentrated at the precession periods (see Figure 8).
[29] The increase of precession period glacial variability with a cooling climate owes to three interrelated effects. First, the direct cooling of the ablation zone restricts ablation to a shorter portion of the summer so that the counterbalancing between summer and fall insolation changes is lost (see sections 2 and 3.2). Second, the cooler climate leads to a thicker ice sheet, which is associated with further cooling of the surface. Finally, the cooler climate also permits ice sheets to reach farther equatorward, and while this would normally be associated with a warmer climate, the insolation at lower latitudes is more sensitive to variations in precession. Conversely, raising the emission level to 7.1 km further suppresses precession period variability, though if the surface temperatures are made still warmer, glaciation is suppressed all together.
[30] Another model run is presented to illustrate how changes in basal conditions can also evoke precession period variability. The deformable till at the base of the glacier is set to zero thickness so that no basal motion occurs. The average slope associated with the leading edge of the ice sheet increases from 0.0003 m/m to 0.02 m/m (significant ablation at the margin sustains this very steep slope). Although the model lacks sufficient resolution to accurately model the ablation zone, such an effect is in the Figure 7. Glacial amplitude as a function of forcing frequency. Circles indicate the ice-volume amplitude resulting from sinusoidal changes in obliquity as a function of frequency. Crosses are similar but for the precession of the equinoxes. The actual obliquity and precession frequencies are 0.025 Ka À1 and 0.05 Ka À1 , respectively. Solid lines indicate the expected change if amplitude is inversely proportional to frequency. correct direction, with the much steeper ascent leading to a smaller ablation zone, and permitting the ice sheet to move equatorward to 50°N. The ice sheet is now almost exclusively sensitive to changes in precession, which accounts for 85% of the ice-volume variance. The loss of ablation zone area also makes the ice sheet less sensitive to insolation variations; changes in precession cause only 0.1% changes in net ice volume. Thus, in the absence of basal  Figures 3a and 3b, showing the time variability of ice volume and its periodogram, respectively. The dashed line is obliquity with its sign reversed and mean and variance scaled. The shaded regions in the periodogram indicate the (left) obliquity and (right) precession bands. The other plots follow but with the emission height changed from 7 Km to (a, b) 7.1 Km, (e, f) 6.5 Km, and (g, h) 6 Km. (i, j) Finally, the emission level is set to 7 Km and basal sediment thickness is set to zero, in order to prevent basal motion. In this case, ice volume is primarily sensitive to precession, as indicated by the dashed line (again with reversed sign and scaled mean and variance). A cooler climate or thicker ice sheet leads to greater precessional influence on ice volume. Note that the y axes have different limits, and that for Figures 8g and 8i, ice volume variability is a small fraction of the total. motion, the model ice sheet is largely insensitive to orbitally induced shifts in insolation. This is consistent with coupled ice-sheet/energy-balance models generally yielding meager variations in ice volume in direct response to insolation shifts [e.g., Oerlemans, 1984], except when unrealistically large ice deformability is specified [e.g., Weertman, 1976;Pollard, 1978]. Furthermore, the apparent need for a thin ice sheet in order to obtain predominantly obliquity period variability may explain why earlier modeling studies [Berger and Loutre, 1997;Nisancioglu, 2004] which simulated the full seasonal cycle and a dynamic ice sheet, but not basal sliding, obtained substantial precession period variability.
[31] The model generates greater ice volume and more precession period variability in response to climate cooling. Thus the progression of glacial variability over the last few million years toward greater ice volume and greater precession period variability [e.g., Huybers, 2007] can be understood as the glacial response to a cooling climate [e.g., Shackleton and Hall, 1984;Raymo, 1994;Ravelo et al., 2004]. The model also generates larger ice sheets and a greater fraction of precession period variability when the ice sheet's basal motion is stopped [e.g., Clark and Pollard, 1998]. The model does not permit differentiating between the two scenarios, although they may be related as a cooler climate would tend to suppress basal melting. The model fails, however, to generate some of the key features of late Pleistocene glacial cycles, for example, not sufficiently deglaciating and not generating $100-Ka variability.
[32] There are also some shortcomings in the model's simulation of early Pleistocene glacial variability. Using d 18 O records as a proxy for ice volume, Hagelberg et al. [1991] and Lisiecki and Raymo [2007] have shown that the climate increasingly resides in a glacial state toward the present, and Ashkenazy and Tziperman [2004] and Huybers [2007] showed that the asymmetry between rapid deglaciation and slow reglaciation is present during the early Pleistocene and increases toward the present. The model glacial cycles, however, are symmetric, appearing most similar to the glacial cycles near $2 Ma.
[33] Another shortcoming is that the maximum equatorward extent of the model ice sheet is only 60°N whereas ice sheets seem to have reached into Iowa and Nebraska (40°N, [Boellstorff, 1978;Barendregt and Irving, 1998]) and there are indications that the Laurentide repeatedly entered the Mississippi River drainage basin during the early Pleistocene [Joyce et al., 1993]. Perhaps the glacial cycles generated by this model represent a simple prototype, but other mechanisms associated with trends toward increasing asymmetry, larger amplitude and longer period glacial cycles and episodes of more southerly glacial extent must also be present. Indeed, a more complex and episodic relationship with insolation forcing seems necessary if one is to explain the obliquity cycle skipping and $100-Ka glacial cycles of the late Pleistocene.

Summary and Conclusions
[34] Insolation integrated over the time-variable summer period, termed the summer energy, provides a framework by which to understand the partition of glacial variability between obliquity and precession periods. This partition is a function of the latitude and ablation threshold. A higher threshold indicates that greater insolation intensity is required to cause ablation. Summer energy predicts that an ice sheet will be most sensitive to obliquity when its ablation zone exists at high latitudes and undergoes ablation at low insolation thresholds. Obliquity control occurs because it has a larger influence on insolation at higher latitudes and the low thresholds permit counterbalancing between the opposing precession-induced insolation anomalies in summer and fall [Huybers, 2006].
[35] An ice-sheet/energy-balance model is constructed to explore the relationship between insolation and ice volume in more detail. The model approximates the modern seasonal cycle in temperature, albedo, and energy fluxes, and when integrated over long time periods, variations in Earth's orbital configuration lead to 40-Ka glacial cycles extending southward as far as 60°N. Ablation is calculated according to a surface energy balance and is shown to initiate at an insolation threshold depending on latitude and elevation. The insolation threshold always resides within the parameter range in which summer energy predicts primarily obliquity period variability. Furthermore, sensitivity to precession is muted because ablation anomalies are counterbalanced between summer and fall, as anticipated from the summer energy.
[36] Another reason why model glaciation is controlled by obliquity owes to its longer period relative to precession. The insolation forcing determines the rate of change of ice volume, not its magnitude, suggesting that the amplitude of the response will be proportional to the forcing period, and as demonstrated by model simulations. Thus ice volume is about twice as sensitive to a 40-Ka obliquity periods as compared to a 20-Ka precession period.
[37] The focus is on the origins of the early Pleistocene 40-Ka glacial cycles, but we hope that study of these seemingly simpler variations will also facilitate understanding of the late Pleistocene glacial cycles. When the longwave emission level of the atmosphere is lowered, temperature decreases, the ice sheet becomes thicker, and it extends farther south, all of which tends to increase precession period variability. Likewise, in the absence of basal deformation, the ice sheet becomes thicker, extends farther south, and is relatively more sensitive to changes in precession. Greater precession period variability and larger ice sheets during the late Pleistocene can thus readily be explained as the response to either a cooler climate or less basal deformation.
[38] The model runs conducted with a cooler climate or with basal sliding switched off also indicate a deficiency in the present model, whereby larger ice sheets grow but never completely deglaciate. Furthermore, the 40-Ka variations produced by the model are symmetric, indicating a mostly linear response to Milankovitch forcing, consistent with the early 40-Ka glacial cycles, but not with the asymmetry of the later glacial cycles [Ashkenazy and Tziperman, 2004]. The model also fails to produce trends toward greater amplitude or to skip one or two obliquity forcing cycles prior to deglaciation [Huybers and Wunsch, 2005;Huybers, 2006], all of which indicates that at least one mechanism is missing.
[39] A number of mechanisms have been suggested to account for the deglaciation of large ice sheets, and which could also introduce asymmetries as well as skipping of obliquity cycles. Oerlemans [1980] showed that if a long time constant is prescribed for isostatic adjustment (10 Ka), rapid deglaciations can be made to occur once an ice sheet has sufficiently depressed its underlying bed, and that this can generate self-sustained $100-Ka glacial cycles. Pollard [1982] speculated that sea-water incursion at the base of an ice sheet or that calving associated with proglacial lakes would increase deglaciation. Using a thermo-mechanical ice sheet model, Oerlemans [1984] suggested that episodic production of basal melt water could lead to rapid collapse of an ice sheet. Peltier and Marshall [1995] called upon dirty snow to lower ice albedo. Gildor and Tziperman [2000] argue sea-ice growth during a glaciated state would diminish accumulation and eventually starve an ice sheet. There is no shortage of mechanisms for generating complete deglaciations, and much of the future work lies in distinguishing between these plausible hypotheses.

Appendix A: Details of the Model
[40] Our aim is to construct a model capable of representing insolation's influence on glaciation. Presently, owing to computational constraints, atmospheric GCMs can usually only take of a series of snapshots of the climate state through a glacial cycle [e.g., Jackson and Broccoli, 2003] or must utilize idealized variations in Earth's orbit [e.g., DeConto and Pollard, 2003a]. Thus, following a long tradition [e.g., Pollard, 1978;North et al., 1981;Shell and Somerville, 2005], we couple an energy balance model (EBM) and an ice-sheet model (see Figure A1).

A1. Energy Balance Model
[41] The EBM explicitly represents temperature at three levels: the middle atmosphere (T a ), surface (T s ), and subsurface (T ss ). Temperature at each level responds to fluxes of energy, C ss @T ss @t ¼ ÀF ss ; ðA3Þ Figure A1. Schematic of the energy fluxes. Levels from top to bottom are the upper and middle atmosphere, atmospheric surface layer, ground/ice surface, and subsurface. Arrows indicate locations at which radiative, diffusive, or turbulent heat fluxes are absorbed or reflected. Note that the atmosphere radiates upward only at the upper atmospheric level. The model has 1°resolution in latitude. Surface and subsurface boxes are represented as either ground or ice and, in this case, an ice sheet extends equatorward to 55°. For the sake of visual clarity, the y axis is not drawn to scale.
where C is the heat capacity associated with each layer, t is time, S is the solar (shortwave) heating, and I is the infrared (longwave) heating. F s is the sensible heat flux from the surface to the atmosphere, and F ss is the heat flux from the subsurface into the surface. Latent heat fluxes and transport of heat within the ground and ice are not represented. There is no ocean, except as implied by the ice sheet calving off at continental margins and the prescription of an ice albedo in keeping with sea ice at the highest latitudes. Values for each constant are given in Table A1.
[42] The incoming solar radiation is variously absorbed and reflected by the atmosphere and surface, depending on the reflectivity (R), absorptivity (A), and transmissivity (T) of the atmosphere and albedo (a) of the surface. Separate albedos are specified for ice (a i ) and bare ground (a g ). The solar energy absorbed by the atmosphere is where the denominator accounts for absorption of radiation reflected by the surface multiple times. The solar radiation absorbed by the surface is, A, R, and T are influenced by clouds and water vapor and are expected to change with a changing climate but here, for simplicity, are parameterized as fixed constants. The sum, A + R + T, must equal 1. Incoming solar radiation, S, is calculated using the approach of Berger [1978] and the orbital solution of Berger and Loutre [1992].
[43] The vertical profile of temperature in the atmosphere is approximated as having a constant moist adiabatic lapse rate, G m . The alternative of specifying a spatially and, possibly, temporally varying lapse rate would require further assumptions regarding the hydrological cycle which we choose to circumvent. The lapse rate is used to calculate the temperature at the atmospheric surface layer, T as = T a + G m H as , where H as is the distance from the middle atmosphere to the surface and which varies according to icesheet thickness and surface depression.
[44] The surface exchanges longwave radiation with the atmosphere according to their respective temperatures, Here a is the atmospheric emissivity, the surface emissivity is assumed to be one, and s is the Stefan-Boltzmann constant. Atmospheric longwave radiation fluxes are given by The final term on the right is the emission temperature of the upper atmosphere, computed as T ul = T a À GH ul , where H ul is the distance from the middle to upper atmosphere. Note that the a in the first term on the r.h.s. of equation (A7) accounts for the fact that some surface longwave radiation is transmitted directly to space. [45] The sensible heat flux, F s , is parameterized using a bulk coefficient, F s = K s (T s À T as ). Likewise, the heat flux from the subsurface into the surface is parameterized as having a linear dependence on the temperature gradient, F ss = K ss (T s À T ss ).
[46] The eddy heat flux divergence depends on the meridional temperature gradient in the middle atmosphere [Stone, 1972], where j.j indicates the absolute value, f is latitude, and K a is tuned to give a reasonable meridional heat flux.
[47] One meter of precipitation falls annually in the model at each grid point and is equally distributed throughout the year. Precipitation collects as ice when the temperature of the atmospheric surface layer is below freezing, and is assumed to run off otherwise. The simplifying assumption is made that evaporation and precipitation occur at the same latitude so that meridional latent heat fluxes are ignored. Ablation occurs when ice is present, the surface is at the melting point, and additional heat is fluxed into it. The amount of ablation is equal to the heat flux per square meter divided by the latent heat of ice per meter cubed. In this manner, ice tends to build up at high latitudes during the winter and then ablate during the spring and summer. Precipitation which accumulates as ice is assumed to give up its excess heat to the atmosphere, so that ice accumulates on the ground with the same temperature as the atmospheric surface layer. A regression between ablation and positive degree days (see equation (1)) calculated from the atmospheric boundary layer temperature yields a slope of 3 mm/ (°C day), in good agreement with generally accepted values [e.g., Braithwaite and Zhang, 2000;Lefebre et al., 2002].

A2. Ice Sheet
[48] The ice-sheet component of the model is zonally averaged and a function of meridional distance, y, and height, z. It utilizes a common shallow-ice approximation [e.g., van der Veen, 1999], assuming that deformation occurs only as a result of horizontal shear stress and that stress and strain are related by Glen's law. The ice is assumed isothermal and incompressible, and the evolution of its thickness is governed by the continuity equation, Here h is the thickness of the ice sheet, B represents the surface mass balance, and u is the vertically averaged velocity. The horizontal shear stress is approximated as a function of ice thickness and surface slope, where r i is ice density and H is the elevation associated with the base of the ice sheet. It follows that the horizontal velocity is where n is the exponent in Glen's flow law and A governs the deformability of the ice. A is known to depend on temperature, fabrics within the ice, and impurities [e.g., Paterson, 1994] but here, for simplicity, is taken as a constant, consistent with ice at a temperature of 270K. The final term, u b , represents the velocity at the base of the ice sheet owing to basal sliding or motion of subglacial sediment. Integrating u with respect to z and inserting into the continuity equation gives an expression for the temporal evolution of the ice sheet, [49] Basal velocity is prescribed according to the sediment deformation model of Jenson et al. [1996] and as used by Clark and Pollard [1998] (their equation (3)). The velocity at the ice-sediment interface is where the vertical bars indicate the absolute value and ''min'' is the minimum of the two quantities. Here h s is the thickness of the sediment at the base of the glacier. The value a sed equals r i gh@(h + H)/@{y} and is the shear stress imparted from the ice into the sediment, and b sed equals (r s À r w )g tan(f s ) and is the rate of increase of shear strength with depth in the sediment. The shear-zone thickness is given by the ratio a sed /b sed , and in the model of Clark and Pollard [1998] is typically between 1 and 10 m. We follow Jenson et al. [1996] and Clark and Pollard [1998] in assigning values for the sediment rheology. The angle of internal deformation is specified as f s = 22°, the reference deformation rate is D o = 7.9 Â 10 À7 s À1 , the exponent m equals 1.25, and the reference viscosity is set m o = 3 Â 10 9 Pa s. The appropriate value of m o is uncertain, and other experiments were conducted increasing and decreasing this value by an order of magnitude, and which yielded results qualitatively similar to those reported here. See Table A2 for a list of ice sheet parameter values.
[50] The scouring of continental regolith is ignored, and a constant sediment thickness is assumed. Furthermore, the sediment is assumed to always be liquid saturated and thus readily deformed. Tracking the temperature and heat transport of the ice sheet itself would permit a more realistic representation of basal melting and hydrology, though in this meridional ice sheet the difficult subject of water drainage would, at best, still be crudely represented.
[51] The adjustment of the bed height to the overlying ice load is modeled as a simple local relaxation toward isostatic equilibrium [e.g., Clark and Pollard, 1998], where H is the height of the bed above sea level, H eq is the equilibrium height of the bed given ice-free conditions (taken as 0 m), h is the height of the ice above the bed, r i is the ice density (910 kg/m 3 ), r b is the bedrock density (3370 kg/m 3 ), and T b is a time constant set to 5 Ka [following Peltier and Marshall, 1995;Clark and Pollard, 1998]. Ice flowing to 85°N is assumed to calve into the ocean and melt.
[52] The ice sheet influences the EBM in a variety of ways: the presence of ice changes the surface albedo in the EBM, the height of the ice sheet influences surface temperature according to the moist adiabatic lapse rate, and the latent heat in the ice can cause large surface-atmosphere temperature contrasts. Likewise, the EBM feeds back onto the ice sheet through determining the surface mass balance. The EBM and ice-sheet components are asynchronously coupled. The EBM is stepped forward at 15-min intervals for 5 years using the Euler method, and the glacial mass balance is calculated from the difference in accumulation and ablation averaged over the final year of the EBM run. The ice model is then integrated forward for 500 years using a semi-implicit Crank-Nicholson formulation at 5-year time steps, after which the energy-balance model is again run to equilibrium with the new ice distribution. Integrations using different time stepping, up to 1-day time steps in the EBM run for 3-to 25-year intervals and 1-to 5-year time steps in the ice-sheet model run for intervals of 50 to 500, yield essentially identical results indicating that the time stepping and coupling techniques are stable. Model runs are started at 2.1 Ma. Runs starting from different initial conditions all converge by 2 Ma.

A3. Model Climate
[53] Under the modern orbital configuration, no ice accumulates and the model equilibrates on a timescale of years. Figure A2 shows the seasonal cycles in various quantities contoured as a function of latitude. Surface temperatures are too warm in the tropics at 40°C, presumably owing to the model lacking evaporation, latent heat transport, and a Hadley circulation. Note that we are not attempting a realistic simulation, but rather to diagnose and isolate mechanisms responsible for the 40-Ka glacial cycles within a model containing only the most pertinent physics. Extratropical temperatures, which are by far the more important for glaciation in this model, are more in keeping with observations, decreasing to an annual average of À15°C at 85°N.
[54] Accumulation occurs throughout the winter in high latitudes, with a longer accumulation season at the highest latitudes. In unglaciated conditions, ablation occurs primarily in the spring and, at high latitudes, extends into the summer, keeping surface cool later into the year. The midtroposphere temperature undergoes a more symmetric seasonal cycles than at the surface, inducing large heat fluxes into the surface during the ablation season. The subsurface temperature adjusts more slowly, having a seasonal cycles which is lagged and attenuated relative to the surface, and tending to warm the surface in the winter and cool it in the summer.
[55] Instrumental records of temperature show that daily average insolation intensity has an almost perfectly linear relationship with zonally averaged surface air temperatures when they are lagged relative to the insolation time series by 30 days [Huybers, 2006]. The model shows a similar relationship, where the cross correlation between insolation and temperature lagged by 30 days is 0.99 or greater at all extratropical latitudes. Note, however, that when parts of the model become glaciated, surface ground temperatures cannot exceed the melting point, tending to suppress the atmospheric temperature and leading to a weaker correlation between insolation and lagged temperature of 0.9.
[56] It is also useful to compute some seasonally averaged quantities which can be directly compared against modern observations [e.g., Peixoto and Oort, 1992]. In the model the meridional transport of heat peaks at 4 Peta-watts at 50°N ( Figure A3b), whereas modern estimates of heat transport peak at 5 Peta-watts near 30°N [Trenberth and Caron, 2001]. The northerly peak in heat transport owes to absorption of heat by the seasonal snow cover, and the weaker heat transport at lower latitudes presumably owes to lack of latent and Hadley cell heat transports in the model.
[57] Two different albedos are specified in the model for ice (0.8) and bare land (0.3). The annual average albedo in the model is 0.3 at low latitudes but increases beginning in the midlatitudes to values of 0.8 at the permanent sea ice parameterized for latitudes above 85°N. The average albedo during winter is significantly higher in the mid and high latitudes owing to the presence of ice cover, and are commensurate with those estimated for Greenland and Antarctica. Estimates of zonal average albedo for the modern climate are similar to those produced by the model [Peixoto and Oort, 1992].
[58] The total absorbed solar radiation ( Figure A3d) is a function of the incident radiation; the albedo of the surface; and the atmospheric transmissivity, reflectivity, and absorptivity. The annual average absorbed solar radiation is 250 W/m 2 in the tropics and steadily decreases to about 100 W/m 2 in the Arctic. The emitted terrestrial radiation in the model ( Figure A3e) is somewhat less than 250 W/m 2 at the equator and somewhat more than 100 W/m 2 at the pole, but values are all within approximately 30 W/m 2 of those estimated for the modern climate [Peixoto and Oort, 1992]. The difference between the absorbed and emitted radiation ( Figure A3f) owes to local heat storage as well as the atmospheric transport of heat. In the arctic, the model atmospheric heat flux convergence is 30 W/m 2 on an annual average basis but can reach as high as 90 W/m 2 during polar night. Overall, we conclude that the EBM gives a Figure A2. Annual cycle in the energy-balance model using the modern orbital configuration. Each quantity is contoured as a function of day of the year (starting with 1 January) and latitude. (a) Insolation at the top of the atmosphere in W/m 2 . (b) Total heat flux across the land/atmosphere boundary, including sensible, shortwave, and longwave heat fluxes in W/m 2 . The densely spaced contour lines are in regions with snow and ice cover during the spring and summer and reach peak values of 220 W/m 2 . (c) Rate of ablation with contours spaced by 2 cm/day. At bottom are the temperatures for (d) the middle atmosphere, (e) surface, and (f) subsurface in degrees Celsius. Figure A3. Model seasonal and annual averages of (a) incoming solar radiation, (b) atmospheric energy transport, (c) albedo, (d) absorbed solar radiation, (e) emitted terrestrial radiation, and (f) net heating. Solid lines are for the annual average, while dashed lines indicate summer (JJA) and dash-dotted lines indicate winter (DFJ) averages. Values indicated in these plots can be directly compared with those estimated for the modern climate by Peixoto and Oort [1992, pp. 127 -129].
sufficiently reasonable representation of the seasonal cycle to permit exploration of the mechanisms connecting insolation to ice-volume fluctuations.