Ice-sheet and atmospheric forcing

The models used here have found standard usage in many ice-age studies; for brevity only a general description is given below, with the model components shown schematically in Figure I. For most of our results we used the perfectly plastic ice-sheet model (e.g. Weertman, Reference Weertman1976), in which the ice-sheet profile versus latitude is constrained to be parabolic. The bedrock depression below the ice is assumed to be isostatic, and the northern ice-sheet tip is fixed at the Arctic Ocean shoreline. The ice sheet varies in size according to the surface accumulation minus ablation on its southern half, always keeping the parabolic profile. As in Weertman (Reference Weertman1976), the atmospheric forcing is simply a prescribed pattern of net annual accumulation minus ablation, which is shifted vertically to represent long-term orbital perturbations, stochastic inter-annual weather variability and ice-sheet albedo feedback. The forcing is shown in the figures below in terms of h ₀, the elevation of the snowline at the Arctic Ocean shoreline, given by

Fig. 1.

Schematic Cross section of the ice sheet. h ₀ is the elevation of the snowline at the Arctic Ocean shoreline.

h _orb is a long-term sinusoidal forcing analogous to the effects of orbital perturbations on the atmospheric climate, and its form for each run can be seen in the figures below. The second term h _alb represents the albedo feedback of the ice sheet, and lowers the snowline in proportion to the current ice-sheet extent.

The third term h _sto in Equation (1) represents the year-to-year variability due to inter-annual weather variations, and is important for the present study since it introduces the possiblity of sampling errors in asynchronous coupling schemes (Harvey, Reference Harvey1986; Hasselmann, Reference Hasselmann and Schlesinger1988). In reality and in GCMs, the inter-annual variability of storm tracks, temperature and other atmospheric fields causes the net mass balance on ice surfaces to fluctuate from year to year, with amplitudes on the order of one-tenth of the glacial-interglacial changes. If in an asynchronous run the synchronous period T _s is set too short, one or two abnormal weather years can skew the extrapolated forcing used throughout the next asynchronous period. h _sto is assigned a random value annually during synchronous periods, distributed normally with zero mean and standard deviation of 100 m and with no correlation from year to year. The amplitude of 100 m is derived from observed summer temperature variability in Oort and Rasmusson (Reference Oort and Rasmusson1971) and winter snow-cover variability in Wiesnet and Matson (Reference Wiesnet and Matson1979), but the most appropriate value is uncertain as discussed below (cf. Oerlemans, Reference Oerlemans1979).

Fig. 2.

(a) Solid curve: ice-sheet extent (distance from northern to southern tip) vs time using the diffusive ice-sheet model with synchronous coupling. Dotted curve: atmospheric forcing h ₀ as defined in Equation (I), (b-d) Differences from Figure 2a using: (b) non-extrapolated asynchronous coupling; (c) extrapolated asynchronous coupling; (d) least-squares asynchronous coupling.

If the southern half of the ice sheet lies either entirely below or entirely above the snowline, the assumption of a fixed parabolic profile is no longer valid (Weertman, Reference Weertman1976; Birchfield, Reference Birchfield1977). We encountered this situation in only one type of run, in which the ice sheet vanishes periodically; as a nascent ice sheet starts to grow it remains entirely above the snowline for a few thousand years, and vice versa as a shrinking ice sheet vanishes. This caused spurious effects in our results, so for that type of run we used a more general diffusive ice-sheet model. (For all other types of runs the results obtained from the two models agreed well.) The diffusive ice-sheet model has been used with variations in many ice-age studies (e.g. Reference Birchfield and GrumbineBirchfield and Grumbine, 1985; Reference Hyde and PeltierHyde and Peltier, 1987). It is based on a vertically integrated ice-flow law, and predicts ice thickness as a general function of latitude. We take the bedrock depression to be isostatic, constrain the northern ice-sheet tip at the Arctic Ocean shoreline, and use the same atmospheric forcing pattern as Birchfield and Grumbine (Reference Birchfield and Grumbine1985).

Examples

Figures 2 and 3 show typical results for sinusoidal forcing h _orb, with T _s = 10 year, T _A = 990 year. The amplitude of h _orb in Figure 2 is large enough to make the ice sheet vanish periodically, which necessitated the use of the diffusive ice-sheet model as explained above. All the error curves in Figure 2 show large spikes when the ice sheet is decaying rapidly (ice extents <∼500 km), because distortions in the forcing applied for one 990 year asynchronous period can have a relatively large effect when the ice is thin. When the ice sheet is large its relaxation time is on the order of 7000 year (Birchfield, Reference Birchfield1977), and the uncorrelated distortions of several consecutive 990 year asynchronous periods usually cancel out before the ice sheet reacts much.

However, spikes of about half the size of those in Figure 2 occur using synchronous coupling alone, just by choosing different realizations of the random number sequence used for h _sto. (When the ice sheet is large, the differences between synchronous realizations are negligible.) Hence only about half of the spike amplitudes in Figure 2 is attributable to asynchronous error, the rest being due to unpredictability inherent in the model.

The errors due to the non-extrapolated scheme in Figures 2b and 3b are systematically positive when the ice sheet is shrinking and negative when it is growing. They are mainly due to the long-term sinusoidal forcing effectively being phase-shifted by an amount T _A, simply because the “true” value of h ₀ computed during T _s is used without change for the next T _A years.

The extrapolated scheme drastically reduces the phase-lag errors (to ∼15 km, not shown), and the remaining errors in Figures 2c and 3c are due mainly to sampling errors related to the random component A_st0 in Equation (1), as discussed further below.

In Figures 2d and 3d, the least-squares method produces noticeable improvement, at least after the first -40 000 year when information from the previous 40 synchronous periods has been accumulated. However, the method works well only when the form of the long-term forcing function is known and smooth as with the sinusoids here; for step-function forcing (not shown), the results are considerably worse than those using the other schemes (cf. Guest, Reference Guest1961, p. 265). The practicality of this method for future GCM runs is also questionable, since it requires the fitting of ∼40 sets of global data. For these reasons we put aside the least-squares method for now, and concentrate on the more viable extrapolated method.

Effect of T _s and Τ _A on stochastic forcing

In other runs not shown here, we have found that the main points deduced above hold generally for various long-term sinusoidal, step-function and real orbital forcing functions. In the remainder of this paper we examine how the errors in the extrapolated method depend on the choices of T _s and T _A. These errors are basically caused by inadequate sampling of the inter-annual variability h _sto in Equation (1) during each synchronous period. In synchronous runs this term is felt by the ice sheet as white noise, and causes red-noise fluctuations in our model ice sheet of about 3 km (cf. Oerlemans, Reference Oerlemans1979; Oerlemans and Van der Veen, Reference Oerlemans and van der Veen1984; it also causes a reduction in mean size of ∼7 km due to the ice-sheet curvature, as seen in Figure 5a). With asynchronous coupling, the forcing during T _A is based on averages taken over the previous two synchronous periods T _s. These averages depart randomly from the true mean climate by an amount

the inter-annual variability, since the standard deviation of the mean of Ν random samples is ] the individual standard deviation. Furthermore, they are applied not just over one year as in synchronous coupling, but over the entire period T ^S + T _A.

Fig. 4.

Errors in ice-sheet extent versus (T _s + T _A)/T _S using extrapolated asynchronous coupling. Each point represents a run of 200 000 year with various types of long-term forcing A_orb and different values of T _s, all with T _s + T _A= 1000 year. The upper and lower shaded regions show maximum (absolute) errors and root-mean-square errors respectively. Points using the same long-term forcing are connected by straight lines, with labels A, B, C, and D indicating the forcings used in Figure 2, Figure 3, Figure 3 with doubled frequency, and Figure 5, respectively. The sporadic small ice-sheet spikes such as in Figure 2c have been excluded from the “A” lines (see text).

For large ice sheets, this can be modeled as linearized departures from equilibrium due to a step-wise discontinuous random forcing;

where l(t) are small perturbations in ice-sheet extent, τ is the ice-sheet time-scale (∼7000 year), and ∂l _eq/∂h is the dependence of equilibrium ice-sheet size on snowline elevation (∼4000). The forcing h′ (t) is composed of linear segments over time intervals t _c, with random start and end values that depend on h _sto and the coupling method. By defining a discrete sequence l _n = l(nt _c) and solving Equation (2) over one interval t _c, the sequence l _n is seen to be a first-order Markov process, and for t _c ≪ τ its variance is (Priestley, Reference Priestley1981; cf. Reference Saltzman, Sutera, Hansen, Berger, Imbrie, Hays, Kukla and SaltzmanSaltzman and others, 1984):

For synchronous coupling,

and t _sto ≡ 1 year to represent observed inter-annual variability, which yields . (The variance of l(t) can be shown to be almost the same as that of l _n, for t _c ≪ τ For extrapolated asynchronous coupling, the forcing variance is reduced by a factor t _sto/T _s and t _c becomes T _s + T _A, so Equation (3) becomes

For example, with T _s = 10 year and T _A = 990 year. Equation (4) yields

Our model results are generally consistent with this analysis, as shown for varying (T _s + T _A)/T _S in Figure 4 and for fixed (T _s + T _A)/T _S in Figure 5. In the lower envelope of Figure 4, the rms errors for a variety of long-term forcing functions are about the predicted size and have a roughly square-root dependence on (T _s + T _Α)/T _s. For a fixed value of(T _s + T _Α)/T _s Figure 5 shows that the magnitude of the errors is fairly independent of the individual T _s and T _A values. However, the value of T _A does affect the spectrum of the ice-sheet response by eliminating periods shorter than ∼T _A, as evident in Figure 5b-d.

Uncertainties

We have varied most model parameter values over plausible ranges and used other atmospheric forcing patterns, and found only slight (∼20%) changes in the asynchronous errors reported here. However, there are other sources of uncertainty as outlined below.

The most appropriate amplitude of the inter-annual forcing h _sto in Equation (1) is not well constrained by the present data. We have used 100 m, but values as high as 200 m are plausible (cf. Oerlemans, Reference Oerlemans1979). From Equation (4) and from other runs not shown, the size of the asynchronous errors with extrapolated coupling is proportional to the amplitude of h _sto, so this source of uncertainty translates to a factor of about 2 in our quantitative results.

The errors are also proportional to ∂l _eq/∂h in Equation (4). This value is relatively small for stable ice sheets as in this study, but it is conceivable that the equilibrium size of real Pleistocene ice sheets is extremely sensitive to snowline changes, so the stabilizing feedback in Equation (2) becomes very weak allowing the stochastic forcing to be an important “random-walk” factor in ice-age evolution (e.g. Nicolis, Reference Nicolis1982). We have not tuned the model to investigate this regime.

We have implicitly assumed that stochastic ice-sheet forcing is not important for the real ice ages, but is a nuisance to be filtered out during synchronous periods T _s. As Hasselmann (Reference Hasselmann and Schlesinger1988) points out, if it really is important its statistics should be measured during each synchronous period so that a synthetic stochastic component can be superimposed during the asynchronous periods. In future GCM runs, that (and the spin-up of the upper oceans) would probably require longer synchronous periods than 10 years.

Concluding Remarks

Errors in asynchronous coupling of simple ice-sheet and atmospheric forcing models are due to (i) the phase-lagging of the long-term “orbital” forcing by ∼T _A, and (ii) the amplification and distortion of the stochastic forcing component.

The phase-lag errors can be reduced to negligible levels by using the extrapolated scheme with T _A on the order of 1000 year or less. The remaining errors are then due to the spuriously altered stochastic forcing, which amplifies the resulting ice-sheet fluctuations by a factor of about

compared to synchronous coupling, and concentrates them at periods of T _A and above.

Fig. 5.

(a) Solid curve: ice-sheet extent using the plastic model and synchronous coupling, with stochastic forcing alone (constant h _orb). Dashed curve: without stochastic forcing, (b-d) Differences from the solid curve in Figure 5a using extrapolated asynchronous coupling with (b) T _s = 1 year, T _s + Τ _A = 100 year; (c) T _s = 10 year, T _s + T _A = 1000 year; (d) T _s = 30 year, T _s + T _A = 3000 year.

Figure 4 shows the general dependence of error size on (T _s + T _A)/T _S, using the extrapolated scheme. The large error spikes such as in Figure 2c have been excluded from Figure 4 on the grounds that they are sporadic, occur only for decaying ice sheets less than ∼500 km in extent, and are partly due to inherent model unpredictability. We consider that non-sporadic errors of ∼150 km or less will be acceptable for future GCM ice-age simulations, since that is an order of magnitude smaller than the full glacial-interglacial range and is a fraction of the typical horizontal resolution in atmospheric GCMs. To keep asynchronous errors below about 150 km, Figure 4 shows that (T _S + T _A)/T _S must be ∼100 or less.

The ratio (T _S + T _A)/T _S is also the approximate increase in efficiency afforded by asynchronous coupling in long-term GCM ice-age simulations, since most of the computer costs will be incurred by the atmospheric and upper oceanic components running only during synchronous periods. Hence the future choice of (T _S + T _A)/T _S will involve a trade-off between efficiency and accuracy. Coupled atmospheric and oceanic GCMs currently use about 10 supercomputer (CRAY X-MP/48) CPU hours per simulated year, which means that a fully synchronous simulation of the last 150 000 years would require 170 years of CPU time. If asynchronous coupling were used with (T _S + T _A)/T _S = 100, “only” 1.7 years of CPU time would be needed. That is still prohibitive, but should become feasible in the 1990s when available computer power increases by one more order of magnitude.

Acknowledgement

The National Center for Atmospheric Research is sponsored by the U.S. National Science Foundation.