Introduction
With more powerful computers in the future, it will become feasible to perform continuous ice-age climate simulations on 1000 to 100 000-year time-scales, using global atmospheric and oceanic general circulation models (GCMs) coupled to three-dimensional ice-sheet models. To date such simulations have only been performed with simplified energy-balance climate models and one- or two-dimensional ice-sheet models (e.g. Reference Neeman, Ohring and JosephNeeman and others, 1988, and references therein). These models are useful in identifying important ice-age mechanisms, but more powerful models will eventually be needed to test more detailed ice-age scenarios.
In these simulations the deep oceans change gradually on time-scales of hundreds of years, and the major ice sheets change gradually on time-scales of thousands of years. It is therefore computationally wasteful to run the atmospheric model over the entire simulation, since hundreds of years would go by with no appreciable change in the atmospheric climate. Instead, the complete model can be run intermittently for relatively short “synchronous” periods (T s), each followed by a much longer “asynchronous” period (T A) during which the ocean and /or ice sheets are run with atmospheric information extrapolated from the previous synchronous period(s). This will cause slight differences in the results from a fully synchronous run, but the differences will be negligible as long as T s is not too short and T A is not too long.
Several asynchronously coupled atmosphere-ocean GCM experiments have been performed since the 1970s, designed mainly to simulate the present climate (e.g. Reference Manabe, Bryan and SpelmanManabe and others, 1979; Reference Washington, Semtner, Meehl, Knight and MayerWashington and others, 1980). More recently simple energy-balance coupled models have been used to evaluate various asynchronous coupling schemes for the atmosphere and ocean (Dickinson, Reference Dickinson1981; Harvey, Reference Harvey1986; Reference Schneider and HarveySchneider and Harvey, 1986), by assessing the differences in results from exact synchronous solutions.
In a similar vein, the present work evaluates asynchronous coupling schemes between ice sheets and the atmosphere. In the following two sections ice-sheet models, atmospheric forcing and asynchronous coupling schemes are described, and resulting variations in ice-sheet extent versus time are presented, emphasizing the differences from the fully synchronous results. We identify the asynchronous scheme that yields the smallest errors for given values of T s and T A, and explore the ranges of T s and Τ A that keep the errors within acceptable bounds. The main conclusions and implications for computer resource requirements are outlined in the final section.
Model Description
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 0, the elevation of the snowline at the Arctic Ocean shoreline, given by
Fig. 1. Schematic Cross section of the ice sheet. h 0 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 0 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).
Types of asynchronous coupling
Our base-line coupling method is fully synchronous, with the “atmospheric state” h 0 computed from Equation (I) at every ice-sheet time step. This method is referred to as “synchronous” below, and differences from the synchronous results caused by other coupling methods are termed “errors”.
Fig. 3. As for Figure 2, except using the plastic ice-sheet model and smaller amplitude in the long-term sinusoidal forcing h orb.
In our simplest asynchronous coupling scheme, the atmospheric state h 0 throughout each asynchronous period T A is taken as its mean value during the previous synchronous period T s. This is referred to below as “non-extrapolated”.
In a slightly more refined scheme, h 0 during each asynchronous period T A is obtained by linearly extrapolating in time from its mean values in the previous two synchronous periods. Harvey (Reference Harvey1986) used a similar scheme for atmosphere-ocean coupling (his TSE), except that he extrapolated quadratically using the previous three synchronous values. Our linear scheme is referred to below as “extrapolated”.
In an effort to minimize the sampling errors due to h sto term in Equation (1), we experimented with a more complex scheme as follows. The mean values of h 0 during many (about 40) synchronous periods are treated as a time series containing a smoothly varying signal plus additive noise. This time series is least-squares fitted by a sum of (about 5) smoothly varying sinusoidal or polynomial functions, which is then extrapolated to obtain h 0 through the next asynchronous period T A. This method is referred to below as “least-squares”.
Results
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 0 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 Ast0 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
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 Aorb 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,
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
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.
