Skip to main content Accessibility help


  • Access



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Asynchronous Coupling of Ice-Sheet and Atmospheric Forcing Models
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Asynchronous Coupling of Ice-Sheet and Atmospheric Forcing Models
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Asynchronous Coupling of Ice-Sheet and Atmospheric Forcing Models
        Available formats
Export citation


Asynchronous coupling schemes between ice sheet and atmospheric forcing models are evaluated for use in long-term ice-age simulations. In these schemes the ice sheet and atmospheric forcing are run together for short synchronous periods (T s), alternating with longer asynchronous periods (T A) during which the ice sheet is run with atmospheric information extrapolated from the previous synchronous period(s). Two simple ice-sheet models are used that predict ice thickness as a function of latitude, and the atmosphere is represented by a prescribed pattern of net annual accumulation minus ablation. The pattern is shifted vertically to represent long-term orbital variations, stochastic inter-annual weather variability and ice-sheet albedo feedback.

Several asynchronous schemes are evaluated by comparing results with those from fully synchronous runs. The best overall results are obtained using a scheme in which the forcing during each asynchronous period is linearly extrapolated from its means in the previous two synchronous periods. Differences from the synchronous results are caused primarily by poor sampling of the stochastic forcing component, which exaggerates the stochastic ice-sheet fluctuations. We examine how these errors depend on T s and T A, and outline implications for GCM ice-age simulations.


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. Neeman 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. Manabe and others, 1979; Washington 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, 1981; Harvey, 1986; Schneider 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, 1976), 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 (1976), 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, 1986; Hasselmann, 1988). 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 (1971) and winter snow-cover variability in Wiesnet and Matson (1979), but the most appropriate value is uncertain as discussed below (cf. Oerlemans, 1979).

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, 1976; Birchfield, 1977). 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. Birchfield and Grumbine, 1985; Hyde 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 (1985).

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 (1986) 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”.



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, 1977), 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, 1961, 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, 1979; Oerlemans and Van der Veen, 1984; 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 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, 1981; cf. Saltzman 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.


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, 1979). 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, 1982). 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 (1988) 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.


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


Birchfield, G.E 1977. A study of the stability of a model continental ice sheet subject to periodic variations in heat input. J. Geophys. Res., 82(31), 49094913.
Birchfield, G.E and Grumbine, R.W. 1985. “Slow” physics of large continental ice sheets and underlying bedrock and its relation to the Pleistocene ice ages. J. Geophys. Res., 90(B13), 11,29411,302.
Dickinson, R.E 1981. Convergence rate and stability of ocean-atmosphere coupling schemes with a zero-dimensional climate model. J. Atmos. Sci., 38(10), 21122120.
Guest, P.G 1961. Numerical methods of curve fitting. Cambridge, Cambridge University Press.
Harvey, L.D.D 1986. Computational efficiency and accuracy of methods for asynchronously coupling atmosphere-ocean climate models. Part II. Testing with a seasonal cycle. J. Phys. Oceanogr., 16(1), 1124.
Hasselmann, Κ 1988. Some problems in the numerical simulation of climate variability using high-resolution coupled methods. In Schlesinger, M.E, ed. Physically-based modelling and simulation of climate and climatic change. Part 1. Dordrecht, etc. Kluwer Academic Publishing Company, 583614.
Hyde, W.T and Peltier, W.R. 1987. Sensitivity experiments with a model of the ice age cycle: the response to Milankovitch forcing. J. Atmos. Sci., 44(10), 13511374.
Manabe, S, Bryan, Κ and Spelman, M.J. 1979. A global ocean-atmosphere climate model with seasonal variation for future studies of climate sensitivity. Dyn. Atmos. Oceans, 3, 393426.
Neeman, B.U, Ohring, G and Joseph, J.H. 1988. The Milankovitch theory and climate sensitivity 2. Interaction between the Northern Hemispheric ice sheets and the climate system. J. Geophys. Res., 93(D9), 11,17511,191.
Nicolis, C 1982. Stochastic aspects of climatic transitions — response to a periodic forcing. Tellus, 34(1), 19.
Oerlemans, J 1979. A model of a stochastically driven ice sheet with planetary wave feedback. Tellus, 31(6), 469477.
Oerlemans, J and van der Veen, C.J. 1984. Ice sheets and climate. Dordrecht, etc., D. Reidel Publishing Company.
Oort, A.H and Rasmusson, E.M. 1971. Atmospheric circulation statistics. NOAA Prof. Pap. 5.
Priestley, M.B 1981. Spectral analysis and time series. Vol. 1. New York, Academic Press.
Saltzman, B, Sutera, A, and Hansen, A.R. 1984. Earth-orbital eccentricity variations and climatic change. In Berger, Α, Imbrie, J, Hays, J, Kukla, G and Saltzman, B eds. Milankovitch and climate. Part 2. Dordrecht, etc., D. Reidel Publishing Company, 615636.
Schneider, S.H and Harvey, L.D.D. 1986. Computational efficiency and accuracy of methods for asynchronously coupling atmosphere-ocean climate models. Part I. Testing with a mean annual model. J. Phys. Oceanogr. 16(1), 310.
Washington, W.M, Semtner, A.J jr, Meehl, G.A, Knight, D.J and Mayer, T.A. 1980. A general circulation experiment with a coupled atmosphere, ocean and sea ice model. J. Phys. Oceanogr., 10(12), 18871908.
Weertman, J 1976. Milankovitch solar radiation variations and ice age ice sheet sizes. Nature, 261(5555), 1720.
Wiesnet, D.R and Matson, M 1979. The satellite-derived Northern Hemisphere snowcover record for the winter of 1977–78. Mon. Weather Rev., 107(7), 928933.