The impact of topographically forced stationary waves on local ice-sheet climate

Abstract

A linear two-level atmospheric model is employed to study the influence of ice-sheet topography on atmospheric stationary waves. In particular, the stationary-wave-induced temperature anomaly is considered locally over a single ice-sheet topography, which is computed using the plastic approximation. It is found that stationary waves induce a local cooling which increases linearly with the ice volume for ice sheets of horizontal extents smaller than ∼1400 km. Beyond this horizontal scale, the dependence of stationary-wave-induced cooling on the ice volume becomes gradually weaker. For a certain ice-sheet size, and for small changes of the surface zonal wind, it is further shown that the strength of the local stationary-wave-induced cooling is proportional to the basic state meridional temperature gradient multiplied by the vertical stratification in the atmosphere. These results are of importance for the nature of the feedback between ice sheets and stationary waves, and may also serve as a basis for parameterizing this feedback in ice-sheet model simulations (e.g. through the Pleistocene glacial/interglacial cycles).

Introduction

At the Last Glacial Maximum (LGM; ∼18–20 ka BP) the Northern Hemisphere was occupied by two major ice sheets, both absent today: the Laurentide ice sheet over North America and the Eurasian ice sheet. The LGM extents of these ice sheets are well constrained by observational data (Clark and Mix, 2002; Dyke, 2004), whereas their evolution to their maximum extents remains unconstrained (Dyke and others, 2002; Kleman and others, 2002). In addition to the orbital variations of the Earth (Milankovitch, 1930; Berger, 1978), feedbacks between the ice sheets and their environment must be taken into account when considering the evolution of the ice sheets through the last glacial/interglacial cycle. A number of early studies have shown that the high albedo and the elevation of the ice sheets are important feedbacks that further support ice growth (e.g. Weertman, 1976; Källén and others, 1979; Oerlemans, 1980). Furthermore, the ice sheets constitute topographic barriers, which influence the mid-latitude atmospheric circulation (Kutzbach and Guetter, 1986; Cook and Held, 1988; Kageyama and Valdes, 2000; Justino and others, 2005; Abe-Ouchi and others, 2007). In turn, the changes in the circulation patterns may alter the temperature and precipitation fields over the ice sheets themselves (Oerlemans, 1979; Lindeman and Oerlemans, 1987; Hall and others, 1996; Roe and Lindzen, 2001a,b). Some studies even suggest that variations in the orbital parameters influence the atmospheric circulation, which may be an important mechanism supporting the onset of the Northern Hemisphere glaciations (Kageyama and others, 2004; Cubasch and others, 2006; Otieno and Bromwich, 2009). The presence of ice sheets may also shift the precipitation pattern to an increased dominance of upslope precipitation (Sanberg and Oerlemans, 1983; Roe, 2005; Van den Berg and others, 2008).

Several studies have successfully simulated the last glacial/interglacial cycle in reasonable accord with observational data (e.g. Pollard, 1982; Tarasov and Peltier, 1997, 2004; Paillard, 1998; Bonelli and others, 2009). However, uncertainties associated with, for example, ice dynamics, basal sliding and the atmospheric state still remain. Simulations of the last glacial/interglacial cycle using dynamical ice-sheet models forced by climatology, computed by Atmospheric General Circulation models (AGCMs), show the sensitivity to the atmospheric state (Abe-Ouchi and others, 2007; Charbit and others, 2007). Charbit and others (2007) found that the simulated ice volumes and spatial extents of the ice sheets were highly dependent on the AGCM used to force the ice-sheet model. Depending on the AGCM, the maximum difference in simulated ice volume in the Northern Hemisphere at the LGM was 30 × 10¹⁵ m³, which is more than half the estimated ice-equivalent sea-level reduction at the LGM.

While the influence of albedo and ice-elevation feedbacks on the surface mass balance (SMB) of an ice sheet are relatively straightforward, the feedbacks from atmospheric circulation may be more complex. Variations of the zonal mean wind have an impact on the amplitude and phase of the temperature anomalies induced by the topographically forced stationary waves due to the presence of an ice sheet (Hoskins and Karoly, 1981; Held, 1983). This may also have an impact on the SMB of ice sheets. An illuminating study (Roe and Lindzen, 2001a) examined how a single ice sheet on an idealized continent evolved from a regional-scale initial size to a continental-scale equilibrium size using a coupled stationary-wave ice-sheet model. This study suggests that features of the ice sheet at equilibrium are strongly shaped by temperature anomalies created by topographically forced stationary waves, which serve to increase its equilibrium extent. Roe and Lindzen primarily focused on how the interactions between the stationary waves and the ice sheet affected its equilibrium features. Given the complex dependence of the stationary-wave-induced temperature anomalies over the ice sheet on its spatial structure, however, the nature of the feedback between the stationary waves and ice-sheet topography may vary qualitatively throughout the evolution of the ice sheet.

In this study, we examine the stationary-wave response to changes in ice-sheet topography in a linear, steady, quasigeostrophic two-level model. We have deliberately chosen a highly idealized atmospheric model and a plastic ice-sheet topography to allow for qualitative analyses and analytical solutions of the stationary-wave response over a broad range of parameters. The emphasis of this study is on the influence of ice-sheet topography on local stationary-wave-induced temperature anomalies, rather than on fully coupled interactions between stationary waves and ice sheets that was addressed by Roe and Lindzen (2001a). The present study, which does not encompass any surface mass-balance calculations for the ice sheets, aims to derive relations for the stationary-wave-induced temperature anomalies over ice sheets that can serve as a basis to parameterize the stationary-wave feedback in ice-sheet models. One central result is that the mean stationary-wave-induced cooling over the ice sheet increases linearly with ice volume for small to intermediate-sized ice sheets. In the context of a one-dimensional (1-D) ice-sheet model, Roe and Lindzen (2001b) proposed that the stationary-wave-induced temperature anomaly should be proportional to the maximum height of the ice sheet. Our results suggest that a near-linear dependence of the temperature anomaly on the ice volume is appropriate for a two-dimensional (2-D) ice sheet; for a 1-D ice sheet we obtain a near-linear dependence on the ice volume per unit width.

Model Formulation and Boundary Conditions

Ice-sheet topography

The ice-sheet profile is calculated under the assumption of perfect plasticity (Van der Veen, 1999). As a result of the plastic approximation, the ice sheet deforms instantly when the applied stress exceeds a critical value given by the yield stress, τ ₀. This implies that the stress within the ice sheet never exceeds the yield stress and that the ice thickness can be related to its horizontal extent, here measured as the half-length of the ice sheet, L, in both horizontal directions. The assumption of perfect plasticity, which is a first-order approximation of ice-sheet topography, has been successfully used in previous studies to obtain analytical expressions of ice-sheet/climate interactions (e.g. Weertman, 1976; Källén and others, 1979). Accounting for local isostatic depression, the maximum height of the ice-sheet surface, η ₀, and the maximum ice thickness, H ₀, are related to the horizontal extent as

(1)

(2)

Here σ ≡ 2 τ ₀/(ρ _i g) relates the height of the ice sheet to its length scale, and ζ ≡ ρ _i/(ρ _m − ρ _i) yields the bedrock depression, where ρ _i is the ice density, ρ _m the mantle density and g gravitational acceleration. For each gridpoint, the height of the ice-sheet surface is given by

(3)

where x _d and y _d represent the distance from the centre of the ice sheet to x and y, respectively, and

(4)

where p _d is centre distance. Finally, the total ice volume, V, is given by

(5)

where H _m and η _m represent the mean ice thickness and height of the ice surface, respectively. Integrating Equation (3) from 0 to L in both the x and y directions, it can be shown that (H _m, η _m) = 4/(H ₀, η ₀).

Stationary-wave model

The atmospheric response to ice-sheet topography is examined in a linear two-level quasi-geostrophic model on a β-plane channel. The model is schematically outlined in Figure 1. Note that, although the two-level model is highly idealized, it has been widely used in the past as it enables analytical solutions. The rationale for choosing a linear model with topographic forcing is partly based on the results obtained by Cook and Held (1988), who suggested that the winter stationary-wave pattern over North America at the LGM was primarily a linear response to topographic forcing. Unfortunately, the linear model cannot resolve interactions between ice sheets and transient waves acting on synoptic timescales (3–7 days). Topographically forced changes of these waves can influence both the precipitation and the temperature variability, two features that could have a significant impact on the ice-sheet mass balance. In Cartesian coordinates, using the rigid-lid approximation, the barotropic and baroclinic stream function anomalies, ψ _M and ψ _T, of the linear quasi-geostrophic two-level model take the form

(6)

(7)

where η(x, y) is the height to the ice-sheet surface and x and y are the Cartesian coordinates in the zonal and meridional directions, respectively. Equations (6) and (7) are identical to Holton’s (2004) equations (8.15) and (8.16), except for the surface-topography terms on the right-hand sides of our equations. These terms enter the equations through the lower boundary condition of the vertical velocity, w, which is nonzero in our case, due the non-flat surface. The remaining notations of Equations (6) and (7) are the Laplace operator, ∇², the geometric model depth, H, the Coriolis parameter, f ₀, and β ₀, which is the meridional gradient of f ₀. Both f ₀ and β ₀ are calculated at 50° N. The exact choice of this latitude does not change the conclusions of this study. The Rossby radius of deformation, L _d, is given by

(8)

where is the Brunt–Väisälä frequency (θ is the potential temperature). Hence the value of L _d increases with the vertical stability of the atmosphere. The subscripts M and T denote barotropic and baroclinic quantities, respectively, and are related to each model level as

(9)

(10)

where a is an arbitrary variable, and subscripts 1 and 2 denote the upper and lower model levels, respectively. Using this notation, U _M is the vertically averaged zonal mean wind, and U _T is the thermal wind (i.e. the vertical wind shear scaled by 2). Both U _M and _{U T} are spatially uniform and related to the corresponding stream function through the geostrophic approximation.

Fig. 1. The two-level model configuration. The stream function anomalies, ψ ₁ and ψ ₂, correspond to the 250 and 750 hPa levels, respectively. The barotropic and baroclinic stream function anomalies, ψ _M and ψ _T, are given by ψ _M = (ψ ₁ + ψ ₂)/2 and ψ _T = (ψ ₁ − ψ ₂)/2. The temperature anomaly at 500 hPa, T′, is proportional to ψ _T, and the topographic forcing is represented by the height of the ice-sheet surface, η.

As the timescale of ice-sheet expansion is much longer than that of the atmosphere, we neglect the local time derivatives in Equations (6) and (7). The resulting steady-state barotropic and baroclinic equations are then given by

(11)

(12)

A typical feature of the linear quasi-geostrophic steady-state equations ((9) and (10)) is the occurrence of a resonance singularity at the critical stationary zonal wavenumber, which depends mainly on the zonal mean wind (Held, 1983). The singularity, which gives rise to an infinite stream function response, is removed by adding linear damping on vorticity with the damping timescale, τ _F, to the lower model equation. We have also included radiative damping in the baroclinic equation (10), which relaxes the topographically forced temperature anomalies towards a zonal mean state, with timescale τ _T.

Through the hydrostatic equation, the temperature anomaly at the 500 hPa level, T′, is given by:

(13)

where R is the gas constant for dry air. The background meridional temperature gradient, , is given by the thermal wind relation as

(14)

We assume that U _M and U _T are related to each other as U _T = 0.4U _M, which is in rough agreement with the annual mean present-day climate at mid-latitudes in the Northern Hemisphere (Peixoto and Oort, 1992). The equations above are solved using standard fast Fourier transforms on a domain which is periodic in the x direction (here interpreted as longitude). In the meridional direction, the boundary conditions ψ _M = 0 and ψ _T = 0 are applied at 90° N and 0°. We use 512 gridpoints in both horizontal directions.

Experimental design

The boundary condition of the two-level model is represented by a single ice sheet. Due to the periodic nature of the atmospheric model in the zonal direction, the east–west location of the ice sheet is irrelevant. The northern margin of the ice sheet is chosen to be fixed at y _N = 75° N, implying that the position of the southern margin is determined by the ice-sheet half-length, L.

To calculate atmospheric temperature anomalies of reasonable amplitude in the two-level model, we need to assign a proper atmospheric basic state, which is determined by the zonal mean wind, U _M, the vertical wind shear, U _T, and the deformation radius, L _d. Because the basic state in a glacial climate is not precisely known, we need to find combinations of U _M, U _T and L _d that seem physically reasonable. To proceed, we consider the following two atmospheric basic states:

• Present-day: U _M = 9 m s⁻¹; U _T = 3:6 m s⁻¹; L _d = 700 km.

• Glacial: U _M = 12 m s⁻¹; U _T = 4:8 m s⁻¹; L _d = 810 km.

Roughly, the present-day state represents annual mean present-day conditions in mid- to high latitudes in the Northern Hemisphere. The annual variations of the basic state parameters at these latitudes are quite small; the jet stream is weaker in the summer but compensated by a northward shift (Peixoto and Oort, 1992). Because most ablation occurs in the summer months, the stationary-wave influence on the local ice-sheet climate should be interpreted as the summer response. The glacial state is calculated using a stronger background meridional temperature gradient. This is a rough estimate of the atmospheric conditions during periods with extensive ice sheets in the Northern Hemisphere, such as at the LGM (e.g. Cook and Held, 1988). The value of L _d for both states is calculated under the assumption that the atmospheric state, on longer timescales, is near the critical shear of baroclinic instability, which is given by 2U _T = β ₀ L _d ⁻² (Stone, 1978). The maximum amplitude of the stationary-wave-induced temperature over a continental-scale ice sheet of similar extent to the Laurentide ice sheet at the LGM is −4.7°C for the present-day basic state and −5.7°C for the glacial basic state, which is the same order of magnitude as simulated by Roe and Lindzen (2001a).

Due to the sparse vertical resolution of the two-level model, the temperature anomalies computed at the 500 hPa level (∼6 km height) should be interpreted as temperature anomalies at the ice-sheet surface, which reaches a maximum elevation of 3.7 km in this study. Due to the uncertainties associated with the meridional distribution of the local temperature anomaly in channel models, we use the ice-sheet area-averaged temperature anomaly, , to monitor the net stationary-wave influence on the temperature over the ice sheet. The numerical values of the parameters used in this study are given in Table 1.

Table 1. Numerical values of the parameters used

The Scale-Dependent Interaction Between Ice Sheets and Stationary Waves

Here, we analyse qualitatively how temperature and snowfall anomalies forced by the ice sheet can affect its mass balance and evolution. In particular, we examine the scale dependence of effects that act either to change the ice volume or to reorganize the ice-sheet structure. Conceptually, the rate of change of the surface height of an ice sheet can be expressed as

(15)

where T _s is the surface temperature. Here the mass-balance function, f, contains all processes that locally affect the height of the ice sheet, including accumulation, ablation and ice dynamics. Note that for a perfectly plastic ice sheet, the shape and the volume are independent of the local mass balance. However, in this conceptual analysis we allow for departures from the perfect plastic behaviour that enforces a parabolic ice-sheet profile. To qualitatively analyse the interaction between the stationary waves and the ice sheet, we assume the stationary-wave-induced temperature anomaly, T′, and the topographically forced vertical velocity, w, give rise to small perturbations on the surface mass balance of an ice sheet that is essentially in equilibrium, i.e. f ≈ 0. Hence, we crudely approximate Equation (15) as

(16)

where p ≡ − ∂f/∂T _s and q ≡ ∂f/∂w are positive constants converting the temperature anomaly and the vertical velocity to actual mass balance. The rationale for Equation (15) is that we assume the temperature anomalies primarily influence the surface mass balance via ablation (Roe and Lindzen, 2001a,b) and that the accumulation increases with vertical velocity, resulting in enhanced upslope precipitation (Sanberg and Oerlemans, 1983; Roe and Lindzen, 2001a). In reality, accumulation can also be influenced by the stationary-wave-induced temperature anomalies, an effect that is omitted here for the sake of simplicity. Note that the simplistic temperature dependence of the ablation in Equation (15) assumes some pre-existing ablation, which a small temperature perturbation may either enhance or reduce. Essentially, T′ > 0 implies ice-elevation decrease due to enhanced melting, whereas T′ < 0 implies ice-elevation increase due to reduced melting.

The isolated effect of our crude representation of topographically induced snowfall is to move the ice sheet upwind. Using the fact that w = U ₂(∂η/∂x) in Equation (15), we obtain the advection equation ∂η/∂t − qU ₂(∂η/∂x) = 0, describing a shape-preserving upwind translation of the ice sheet. Note that this simple linear accumulation formula has no net effect on the ice-sheet mass balance. In reality, the topographically induced snowfall is a nonlinear function of the vertical velocity and tends generally to enhance the accumulation over an ice sheet (e.g. Sanberg and Oerlemans, 1983; Roe and Lindzen, 2001a). We can crudely estimate the implied upwind translation speed, qU ₂, resulting from the accumulation formula of Roe and Lindzen (2001a, fig. 4). By linearizing this formula around w = 0 and taking the surface temperature to be at freezing point, we obtain q ≈ 1:5 × 10⁻⁶. For a surface wind speed of 5 m s⁻¹, the upwind translation speed is found to be ∼250 m a⁻¹. Thus, this crude consideration suggests that topographically induced snowfall can induce a tendency for the ice sheets to move upwind with a speed of a few hundred m a⁻¹.

To examine the effect on the ice sheet of the simple ablation representation of Equation (16), we set q = 0, and consider a single Fourier component of the ice sheet and the forced response

(17)

where k and l are the zonal and meridional wavenumbers, respectively. The time, t, is associated with the growth and decay of the ice sheet, which is much longer than the timescale of the atmosphere. Inserting Equation (17) into Equations (11) and (12), multiplying by i/k and using Equation (13), we obtain after some algebra (see Appendix):

(18)

where μ ≡ 4f ₀ ²(U _M − U _T)/RH. The function χ(k, l) relates the Fourier component of the ice sheet to the corresponding temperature anomaly. The analytical expression for this function is given in the Appendix. Substituting Equation (18) into Equation (16) yields

(19)

which has solutions of the form

(20)

As p and μ are positive constants, Fourier components for which Re(χ) > 0 experience a positive feedback from the ablation anomalies; in the absence of stabilizing feedbacks they would grow exponentially. Wavenumbers for which Re(χ) < 0, conversely, tend to decay. For a hypothetical infinite ice sheet characterized by a single Fourier component, Re(χ) > 0 would correspond to a situation with cold/warm temperature anomalies at the ice-sheet crests/troughs. The associated pattern of decreased and increased ablation will act to increase the amplitude of the sinusoidal ice-sheet wave.

The quantity pμ Im(χ) can be interpreted as the frequency of the ice-sheet Fourier components, implying that the associated phase velocity in the zonal direction is given by

(21)

and the zonal group velocity is given by

(22)

For a hypothetical sinusoidal ice sheet, Im(χ) < 0 implies negative temperature anomalies and hence reduced ablation on the eastern (leeward) side of the crests; on the western (windward) side of the crests the opposite is true. This ablation pattern will act to advance the ice sheet eastward, i.e. inducing an eastward phase velocity. Hence, the phase shift between the ice sheet and the temperature anomalies, determined by Im(χ), affects the spatial evolution of the ice sheet but not its volume.

The scale-dependent interactions between the ice sheet and the stationary-wave-induced temperature anomalies are illustrated in Figure 2, which shows the real and the imaginary part of χ(k, l). As a response to the ablation anomalies, low, as well as high, ice-sheet wavenumbers tend be amplified, whereas intermediate wavenumbers tend to be attenuated (Fig. 2a). The positive interval of the low wavenumbers terminates at the wavenumber of the stationary equivalent barotropic Rossby wave (say K _s; grey lines in Fig. 2a), beyond which an interval with negative wave-numbers is encountered. The stationary Rossby waves tend to dominate the atmospheric response (Charney and Eliassen, 1949; Held, 1983). When U _T = 0, the wavenumber of the stationary Rossby wave is K _s ≡ (k _s ² + l _s ²)^1/2 = (β ₀ U _M ⁻¹)^1/2; for typical atmospheric values of U _T and L _d, K _s is slightly smaller. Note that it is for length scales comparable to the stationary wave that the ice sheet experiences the strongest growth due to the stationary-wave-induced ablation. The shift from negative temperature anomalies to positive anomalies at k _s is essentially a barotropic phenomenon; an analogous transition from negative to positive stream-function response occurs in a barotropic model. The shift back to negative temperature anomalies (i.e. Re(χ) > 0) at higher wavenumbers occurs when zonal advection of relative vorticity begins to dominate.

Fig. 2. (a) Real and (b) imaginary parts of the transfer function, χ (Equation (17)), for the following atmospheric basic states: ‘present day’ (solid curve; U _M = 9 m s⁻¹, U _T = 3:6 m s⁻¹ and L _d = 700 km) and ‘glacial’ (dashed curve; U _M = 12 m s⁻¹, U _T = 4:8 m s⁻¹ and L _d = 810 km) as a function of zonal wavenumber, k (meridional wavenumber is set to l= 1). The stationary Rossby wavenumbers, k _s, for the present-day and glacial basic states are indicated by the grey vertical solid and dashed lines, respectively. The values of χ are set nondimensional by the scale factor U _M L _d ⁻². Re(χ) < 0 corresponds to a positive stationary wave-induced temperature anomaly over the ice sheet, and Im(χ) < 0 implies a positive phase speed, i.e. eastward ice-sheet propagation.

The phase shift of the stationary-wave-induced temperature anomalies, determined by Im(χ), results in an eastward downstream advancement of the phase for all wavenumbers (Equation (21); Fig. 2). This results from enhanced melting due to positive anomalies on the upstream side and reduced melting due to negative anomalies on the downstream side of the ice-sheet crests, which is the typical barotropic stationary-wave response over a topographic barrier (Held, 1983). Further, Figure 2 reveals that the group velocity can be downstream as well as upstream depending on the wavenumber. The strongest growing wavenumbers, encountered near the wavenumber of the stationary Rossby wave, have a downstream group velocity. This tendency of downwind advancement due to ablation will compete with the effect due to topographically induced snowfall, which acts to advance the ice sheet upwind.

Earlier, we estimated that topographically induced snowfall can advance an ice sheet upwind at a speed of a few hundred m a⁻¹ and further noted that this speed is independent of the wavenumber. We now attempt a rough estimate of the downwind phase speed of the ice sheet caused by the ablation anomalies due to the stationary waves. For this purpose we use the ablation formula presented by Roe and Lindzen (2001b), which yields p = 1.2 m a⁻¹ °C⁻¹. Substituting this value into Equation (20) indicates phase speeds of the ice-sheet wavenumbers between 5000 and 15 000 m a⁻¹ for the largest scales (k < 6), and approaching zero for smaller scales. This crude analysis suggests that the ablation-induced phase speeds of long waves, which are directed downwind, can be an order of magnitude greater than the upwind phase speed that is induced by upslope snowfall. This very crude estimate of the ablation-induced phase speeds is probably on the high side, because it assumes some background ablation everywhere over the ice sheet all year round. Thus, if, for instance, ablation occurs during a quarter of the year, the phase speeds would be reduced by 75%.

The Influence of Ice Volume on the Local Stationary-Wave-Induced Temperature Anomaly

In the previous section, we found that different scales of the stationary waves, described by the function χ, act to either enhance or reduce ice-sheet ablation. In this section, we consider how the shape and volume of the ice sheet influence the local stationary-wave-induced temperature response. The plastic nature of the ice sheets is exploited in the analysis.

Figure 3 shows the stationary-wave-induced temperature anomaly, averaged over the ice sheet as a function of the ice volume. The temperature in Figure 3 is nondimensionalized by the multiplication of . In addition to the present-day and glacial basic states, we include two simulations with different values of L _d (the grey curves in the figure). Because the deformation radius, L _d, is associated with the vertical stratification in the atmosphere (Equation (8)), these additional simulations serve to illustrate the sensitivity to changes of the lapse rate in the glacial atmosphere, an issue that is still uncertain for the LGM (Abe-Ouchi and others, 2007). The larger value of L _d corresponds to a strongly stratified atmosphere, implying a smaller lapse rate. Equivalently, the smaller value of L _d correponds to a larger lapse rate. The scaled stationarywave-induced temperature exhibits nearly the same response for all basic states in Figure 3. This means that the scaled temperature anomaly is proportional to T_yL _d, provided the surface zonal wind, which is proportional to the temperature anomaly (Equation (18)), remains roughly constant.

Fig. 3. The scaled area-mean stationary-wave-induced temperature anomaly, , as a function of ice volume, V, for two atmospheric basic states: ‘present-day’ (solid curve; U _M = 9 m s⁻¹, U _T = 3.6 m s⁻¹ and L _d = 700 km) and ‘glacial’ (dashed curve; U _M = 12 m s⁻¹, U _T = 4.8 m s⁻¹ and L _d = 810 km). The grey curves represent the same basic state as ‘glacial’, with the exception that L _d = 600 km (solid grey curve) and L _d = 1000 km (dashed grey curve). These states are included to illustrate the sensitivity of the stationary temperature anomaly to changes in the atmospheric lapse rate.

To further examine the temperature anomaly changes with ice volume in Figure 3, we consider the ice-sheet Fourier components (Equation (18)). To begin, it is illustrative to consider a general ice-sheet formulation, defined as

(23)

Here, η ₀ is the maximum height of the ice-sheet surface. The extent of the ice sheet is constrained by its half-length, L, and its shape is given by η _*. In the zonal direction, the Fourier coefficients of Equation (23) are given by

(24)

where R _e is the circumference of the Earth at a certain latitude and k_n = 2πn/R _e. Substituting x/L by r and rearranging

(25)

where a ≡ k_nL and

(26)

Equation (25) states that the Fourier coefficients of the 1-D ice-sheet topography (and hence the Fourier coefficients of the temperature anomalies) are proportional to the horizontal extent in addition to the maximum height. The 2-D Fourier representation of the ice sheet takes the form

(27)

where l_n is the meridional wavenumber, b ≡ l_nL, and G(b) is proportional to F(b). Relating the mean height of the ice sheet to its maximum height as η _m = γη ₀, where 0 < γ < 1, Equation (27) can be rewritten as

(28)

Thus, for any ice-sheet representation, the Fourier coefficients of the ice-sheet topography are proportional to the ice volume multiplied by the functions F(a) and G(b), which are determined by the shape of the ice sheet. Using the first-order plastic approximation of ice-sheet topography, the general formulation of the ice-sheet Fourier coefficients can be substituted by γ = 4(1 + ζ)/9 and , implying that:

(29)

Using Equation (18), we can relate the Fourier components of the plastic ice sheet to the corresponding temperature anomalies:

(30)

where F(a) is now given by

(31)

This function can be written using the Fresnel integrals, C and S, as

(32)

where F(0) = 4/3. Function F is shown in Figure 4. For small arguments, the functions F and G decrease only slowly, implying a regime where the lowest Fourier coefficient of the ice sheet is roughly proportional to V. Specifically, F decreases by only 3% in the interval 0 < a < 0.5. Thus for a < 0.5, a linear relationship between the first few Fourier coefficients of the ice sheet and ice volume is a good approximation. At 50° latitude this condition applies when L < 2000/n km or equivalently when

(33)

For the parameters given in Table 1 (σ = 10 m, ζ = 0.4), there is a linear relationship between and V when V < 40 × 10⁶/n ^5/2 km³. This is illustrated in Figure 5: for wavenumber n = 2, increases linearly with V until V = 7.1 × 10⁶ km³, which is equivalent to L = 1000 km. For increasing zonal wavenumbers, the necessary condition for linearity (Equation (32)) is strongly constrained by the n ^−5/2 dependence. As a consequence, for n = 3 and n = 4,the linear relationship breaks down at about V = 2.6 × 10⁶ km³ (L = 700 km) and V = 1.2 × 10⁶ km³ (L = 500 km), respectively.

Fig. 4. Solution of the Fresnel integral defined in Equation (31). Parameter a is given by a = kL, where k is the zonal wavenumber and L is the ice-sheet half-length. When a < 0:5, the value of F(a) changes ∼3%, implying a near-linear relationship between the Fourier amplitudes of ice-sheet topography and ice volume.

Fig. 5. The amplitude of the Fourier coefficients of ice-sheet topography, , as a function of ice-sheet volume, V, for zonal wavenumber k = 1 (solid curve), k = 2 (dashed curve), k = 3 (dashed-dotted curve) and k = 4 (dotted curve). All solutions have been calculated using meridional wavenumber l = 1.

The properties of the spectral plastic ice-sheet representation show that if the response is dominated by a few low wavenumbers, as suggested in Figure 2, then the temperature anomaly should initially increase linearly with ice volume. However, as the ice sheet grows, the increase declines and eventually stops, owing to higher values of a and b. This pattern is analytically described for each Fourier component in Equation (29) and, because the full solution in a linear model is given by the superposition of the individual Fourier components, it is also evident for all basic states in Figure 3. Notably, the present-day basic state temperature anomaly decreases at a slower rate than the other basic states. This occurs because the strongest response is encountered at higher wavenumbers in the present-day case due to a lower value of the zonal mean wind.

The Stationary-Wave Influence on the Temperature Over an Equatorwardexpanding Ice Sheet

In the previous section we showed that if the ice sheet grows, the stationary-wave-induced cooling becomes stronger. In addition, the surface temperature over the ice sheet depends on the latitude as well as the height of the ice sheet. Simple 1-D models have demonstrated that the general temperature decrease with altitude and latitude may give rise to two equilibrium solutions: a small unstable ice sheet and a larger stable one (Weertman, 1976). Here we follow Roe and Lindzen (2001b) and consider also the effect of stationary-wave-induced temperature anomalies on the mean temperature over the ice sheet. For an anchored northern margin, we write , as

(34)

Our representation of the mean surface temperature over the ice sheet contains the temperature at the northern margin, T _N, the surface background meridional temperature gradient, , the lapse rate, Γ, and the area-mean stationary-wave-induced temperature, . Both and Γ are negative, implying that for a southward-expanding ice sheet there is a warming due to the meridional temperature gradient and a cooling due to the atmospheric lapse rate. To examine the different processes contributing to the mean surface-temperature change, we differentiate Equation (33) with respect to V. Using Equation (4) for a plastic ice sheet, we obtain

(35)

where

and

are positive constants containing the ice-sheet height-tolength parameter, σ, and the bedrock sinking parameter, ζ. For an equatorward-expanding ice sheet, the background meridional temperature gradient, (defined negative), leads to warming over the ice sheet. This warming is countered by the lapse rate and the stationary waves, both inducing a cooling over the ice sheet. The relative importance of the processes acting to change the mean temperature depends on the ice volume. The effect of the elevation (hereafter referred to as the ice-elevation effect) on the mean temperature is strongest for small ice sheets and decays rapidly with increasing ice volume. The near-linear dependence of the stationary-wave temperature on ice volume implies a slower decay as the ice volume increases.

Figure 6 shows the derivative of the mean temperature with respect to V, with and without the effect of the stationary waves. In the absence of stationary waves, the temperature initially decreases despite the fact that the ice sheet expands equatorward. The cooling for small ice sheets is related to the ice-elevation effect. The associated reduction of the ablation can destabilize a small ice sheet, prompting it to expand, as pointed out by Weertman (1976). Eventually, at V ∼ 0.3 × 10⁶ km³, the effect of the background meridional temperature gradient begins to dominate, inducing a net temperature increase over the expanding ice sheet.

Fig. 6. The derivative of the mean surface temperature over the ice sheet with respect to ice volume (i.e. the solution of Equation (34)) in the absence of stationary waves (grey curve) and with stationary waves (black curves). The black solid curve represents the case for the stationarywave-induced temperature anomaly computed in the two-level model with the ‘glacial’ basic state (U _M = 12 m s⁻¹, U _T = 4.8 m s⁻¹ and L _d = 810 km). The black dashed curve is calculated using the stationary-wave-induced temperature representation of Roe and Lindzen (2001b). The surface meridional temperature gradient is set to −7°C (1000 km)⁻¹ and the atmospheric lapse rate, Γ, to −6.5°C km⁻¹.

The stationary waves induce a mean cooling over the expanding ice sheet (the glacial basic state: solid curve in Fig. 3). This is especially true for a range of intermediate-sized ice sheets, for which the stationary-wave effect is relatively strong compared with the other processes. There is a regime between V ∼ 1 × 10⁶ km³ and V = 10 × 10⁶ km³ (equivalent to L = 470 and 1200 km, respectively), where the cooling due to stationary waves along with the ice-elevation almost completely cancel out the effect of the meridional temperature gradient on the mean temperature.

For comparison, Figure 6 also illustrates the stationary-wave parameterization suggested by Roe and Lindzen (2001b) (dashed curve). They assumed that the local stationary-wave-induced temperature anomaly is proportional to the maximum height of the ice sheet, i.e. it depends on the ice volume in the same way as the ice-elevation effect. The strength of the stationary-waveinduced cooling in their parameterization is given by −5°C/2 km (equivalent to their CTT case). Compared with the case without any stationary waves, the inclusion of the stationary-wave representation due to Roe and Lindzen enhances the cooling over the equatorward-expanding ice sheet by up to an order of magnitude for larger ice sheets. However, this representation of the stationary-wave-induced temperature anomaly decays faster with ice volume than our results obtained from the two-level model (solid curve in Fig. 6). This implies that the representation of Roe and Lindzen shifts into the warming regime for smaller ice sheets.

Figure 6 suggests that stationary waves can strongly modify the mean temperature over an equatorwardexpanding ice sheet. Potentially, this can result in a strong stationary-wave/ablation feedback, especially for small to intermediate sizes of ice sheets. The strength of this feedback, however, is not entirely determined by the area-mean temperature anomaly. Because ablation occurs at the ice-sheet margin, the stationary-wave feedback depends also on the distribution of the temperature anomaly over the ice sheet. To further examine this issue, one should use a spherical geometry instead of the channel geometry used here.

Summary and Conclusion

In this study, interactions between the steady atmospheric stationary-wave-induced temperature anomalies and ice-sheet topography have been considered in a linear quasigeostrophic two-level model on a β-plane channel. We emphasize that the idealized two-level model neglects several atmospheric processes of importance for the interaction between stationary waves and ice sheets. For example, we have not considered interactions between the topography and the mean flow, a feature that could also influence the local atmospheric wave response (Charney and DeVore, 1979). To obtain the qualitative behaviour of the stationary-wave-induced temperature response, we chose to monitor the stationary-wave feedback using the mean temperature anomaly, averaged over the whole ice-sheet area rather than the ice-sheet margin. However, the local temperature response in a two-level model is quite symmetric around the centre of the ice sheet. Therefore, the qualitative analysis also applies at other parts of the ice sheet (e.g. the southern margin). To study the spatial distribution of the temperature anomaly over an ice sheet in more detail, it is appropriate to use spherical geometry rather than the channel geometry used here. Furthermore, to allow for physical interpretations and analytical solutions, we have used an idealized representation of the ice-sheet topography, based on the plastic approximation (Van der Veen, 1999).

The Fourier analysis presented shows that ablation due to the stationary-wave-induced temperature anomalies acts to propagate the ice sheet downstream. This feature was noted by Lindeman and Oerlemans (1987), who investigated mass-balance perturbations on the LGM ice sheets using a two-level atmospheric model in combination with a statistically based mass-balance model. Further, the positive feedback on the ice sheet due to the ablation anomalies is most pronounced for length scales comparable to that of the stationary Rossby wave; a feature that can be inferred from the equilibrium ice-sheet shapes computed by Roe and Lindzen (2001a), with and without the effect of stationary-wave-induced temperature anomalies. On an infinite ice sheet with constant height, the ablation/stationary-wave feedback will act to generate sinusoidal perturbations of a wavelength close to that of the stationary wave, and these perturbations will propagate eastward (Fig. 2). This hypothetical scenario has similarities with the spatial evolution of sea-surface temperature anomalies that arise from interaction with stationary atmospheric waves (Nilsson, 2001). In the coupled atmosphere/ice-sheet model of Roe and Lindzen (2001a), the scale-selective stationary-wave feedbacks affect how a small initial ice sheet evolves to its equilibrium extent (see their fig. 13): initially, the ice sheet expands westward due to topographically forced accumulation, whereas stationary waves influence the later stages of the evolution, characterized by an eastward expansion of the ice sheet. Ice dynamics also play an important role in ice-sheet evolution in the simulations of Roe and Lindzen (2001a).

A central result of our study is that the area-mean stationary-wave-induced cooling over small to intermediate-sized ice sheets is directly proportional to their ice volume. The underlying reason is that a few low-wavenumber stationary waves dominate the response. As long as their wavelengths are large compared to the ice-sheet extent, the stationary-wave-induced temperature anomaly is proportional to the ice volume. However, the response also depends on the Fourier component of the ice sheet corresponding to the wavenumber of the stationary wave. The amplitude of this Fourier component increases with increasing ice volume. Using the first-order plastic approximation of ice-sheet topography, the stationary-wave-induced cooling increases linearly with ice volume as long as the dominant wavelength of the atmospheric response satisfies the condition λ > 4πL, where L is the half-length of the ice sheet. For an atmospheric response which is dominated by zonal wavenumber 3, we found that λ > 4π L applies for ice sheets with L < 700 km, or, equivalently V < 2.6 × 10⁶ km³ using the standard parameters in this study. For larger ice sheets, the dependence of the local stationary-wave-induced cooling on the ice volume becomes gradually weaker, implying that the linear relationship breaks down.

The proportionality between the stationary-wave-induced temperature and ice volume obtained here is in contrast with Roe and Lindzen (2001b), who parameterized the stationary-wave-induced temperature anomaly as being proportional to the maximum height of the ice-sheet surface. This representation of the stationary waves gives rise to a relatively strong impact on the mean temperature over small equatorward-expanding ice sheets, whereas the representation due to ice volume entails a relatively stronger influence on the mean temperature over larger ice sheets. More specifically, there is a range of intermediate-sized ice sheets (here between V = 1 × 10⁶ and V= 10 × 10⁶ km³) where the stationary waves contribute to a net temperature decrease over an equatorwardexpanding ice sheet.

The amplitude of the topographically forced stationarywave-induced temperature is a function of the size and shape of the ice sheet, as well as the atmospheric basic state. For small changes of the surface zonal wind, we find that the stationary-wave-induced temperature anomaly is proportional to , i.e. the background meridional temperature gradient times the Rossby radius of deformation. As L _d depends on the atmospheric lapse rate, this result yields an estimate of the impact of the background temperature state on the stationary-wave-induced temperature anomaly over the ice sheet. Unfortunately, changes of the atmospheric basic state through a glacial/interglacial cycle remain an unresolved issue (Abe-Ouchi and others, 2007; Charbit and others, 2007).

The present analysis only considered the stationary-wave interactions with a single isolated ice sheet. The Pleistocene glacials, however, are generally characterized by two continental-scale mid-latitude ice sheets in the Northern Hemisphere. Modelling studies suggest that at the LGM the stationary-wave response to the Laurentide ice sheet enhanced the baroclinic eddy activity downstream, which may lead to an intensification of storm tracks and associated precipitation over the Fennoscandian ice sheet (Kageyama and Valdes, 2000). Potentially, the Laurentide ice sheet could also influence the stationary-wave-induced ablation on Fennoscandia. However, the far-field stationary-wave response is not as strong as the local response (Held, 1983), suggesting that teleconnection effects are small compared to effects of stationary waves forced locally by the Fennoscandian ice sheet.

In agreement with Lindeman and Oerlemans (1987) and Roe and Lindzen (2001a), we found that the atmospheric flow response has a leading-order impact on the local climate over ice sheets, manifested by negative temperature anomalies of several degrees. The new important result is that the stationary-wave-induced cooling locally over the ice sheet is proportional to the ice volume and the product of the meridional temperature gradient multiplied by the Rossby radius of deformation. These intriguing new features deserve further attention in a more realistic framework that includes a more complete atmospheric model and a dynamical representation of the ice sheet.