1. Introduction

Experiments with ice-sheet models show that their evolution is particularly sensitive to the parameterization of mass balance, particularly in the ablation zone (Reference OerlemansOerlemans, 1981; Reference HallHyde and Peltier, 1985). Accumulation depends primarily on synoptic scale motions, while ablation depends on the detailed properties of the boundary layer over the ice sheet. Meteorological studies of the mass balance associated with ice sheets have often focused on GCM simulations of climate (e.g. Reference Manabe and BroccolliManabe and Brocolli, 1985; Reference Kutzbach and GuetterKutzbach and Guetter, 1986; Reference Mitchell, Grahame and NeedhamMitchell and others, 1988). The resolution of these models places limitations on their ablility to reproduce boundary-layer characteristics, and their large demands on computer time mean they can only provide snapshots of a steady-state climate for a specific ice-sheet configuration. This is not necessarily the most appropriate approach to take in order to increase understanding of the interaction between mass balance or climate and ice-sheet evolution. On the other hand, models which allow the ice sheet to grow often incorporate highly parameterized climates, which are unable to respond to the evolving ice sheet. In order more fully to understand the way in which mass balance and ablation may change as the ice sheet evolves, a more comprehensive understanding of the boundary-layer processes over ice sheets is required. This paper uses a slab model in order to look at the way in which the atmospheric boundary layer may respond to a change in ice-sheet size and/or shape.

During the last decade or so, measurements over Antarctica (e.g. Reference GosinkGosink, 1982; Reference Kodama and IshikawaKodama and others, 1989) have given greater insight into the katabatic wind regime of the glacial slopes. These observations can be used, in conjunction with numerical models, to increase our understanding of the boundary-layer processes. In the future, this knowledge will be able to provide a basis for improvements to the parameterization of ablation in ice-sheet models. In the absence of any reliable data from elsewhere, the model presented here was developed using the observations from Antarctica in the parameterizations and boundary conditions.

2. The Model

The concepts of the model are shown in Figure 1, The boundary layer is treated as a layer of cold air flowing over a dome-shaped ice sheet, with a surface profile of slope angle α within a stably stratified environment. The boundary layer is subject to sustained cooling, surface stress, entrainment of the air above and pressure gradient and Coriolis forces. Only mean-layer quantities are evaluated, deviations from well-mixed flow being incorporated into a set of profile parameters, as used by Reference Manins and SawfordManins and Sawford (1979). The s-axis follows the surface of the ice in the direction of maximum slope, parallel to the direction of flow component u. The y-axis is directed perpendicularly to the left, in the direction of flow component υ, parallel to the contours, and the n-axis is orientated normal to the s-y plane. The air is treated as a Boussinesq fluid, so that variations in density are ignored, except for those producing a buoyancy force when combined with gravity (Reference Manins and SawfordManins and Sawford, 1979); this approximation is reasonable for heights of a few hundred metres. The upper boundary of the model is constrained by a prescribed potential temperature gradient (dθ/dz) and wind speed (V _g). The temperature gradient implies a Brunt-Vaisalla buoyancy frequency (N ² which is the velocity of propagation of small-amplitude gravity waves. In a stable layer, this is defined as

(1)

where g is the acceleration due to gravity and θ_br is the potential temperature of the “ambient” air and is defined in Figure 1. The wind above the boundary layer (Vg) is assumed to be geostrophic and exclusively in the y-direction. Surface cooling is prescribed according to a constant net-radiative divergence in the boundary layer (∂R _n/∂n). The air within the boundary layer is acted on by a buoyancy force due to the potential temperature deficit (Δθ), in the ambient temperature field (θ _br), and flows down-slope under the influence of this force, combined with the over-riding pressure-gradient force and the Coriolis force. Friction acts at the surface with the ice sheet (surface stresses in the x-and y-directions are τ_x and τ_y , respectively) and entrainment is permitted with the warmer air above the boundary layer, producing additional drag.

Fig. 1. The idealized boundary layer of the slab model.

2.1. The governing equations

The model is based on the following set of equations for the steady-state conservation of momentum, mass and heat, using the nomenclature given in the previous section and taking (p-p _a) as the pressure depth of the boundary layer and ∂p/∂s as the pressure gradient of the ambient air. Momentum

(2)

(3)

Hydrostatic balance

(4)

Heat

(5)

Mass

(6)

The wind above the boundary layer is assumed to be geostrophic and exclusively in the y-direction, so that the first term on the righthand side of Equation (3) allows for the velocity gradient above the boundary layer. It is assumed that this is a linear profile so that (∂V _g/∂n) ≈ (V _g/h). The second term on the righthand side of Equation (2) can be replaced by the pressure-gradient force of the free atmosphere and the effect of the thermal wind between the free atmosphere and the boundary layer, so that

(7)

Neglecting vertical accelerations with respect to the pressure variation, the hydrostatic equation can be rewritten as follows:

(8)

This can be integrated with respect to n and substituted into Equation (2) using Equation (7).

2.1.1. The profile factors

The vertical profiles of and θ are given in Figure 2. The integration of Equations (2)–(6) assumes well-mixed flow with average mean layer-velocity components U and V and, following Reference Ellison and TurnerEllison and Turner(1959) and Reference Manins and SawfordManins and Sawford (1979), profile factors are introduced to scale the integrals in order to account for deviations from well-mixed flow. These profile factors (S₁,S ₂ and S ₃) are defined by the following equations:

(9)

(10)

(11)

where h is the depth of the boundary layer, scaled using the profile factors in Equations (9)–(11), H is the height assumed to be unaffected by katabatic flow and W _e is the entrainment velocity. If there was no entrainment, h would be at the same height as H (personal communication from J. P. Gosink, 1991), and there would be a discontinuity at the top of the boundary layer. For well-mixed katabatic flow, S ₁ = S ₂ = S ₃ = 1. The profile factors account for the distribution of the density relative to the buoyancy profiles, and arise because integrating each according to average characteristics of the layer is not strictly equivalent to summing the integrals of each of the sub-layers within the boundary layer.

Fig. 2. Schematic of the vertical profiles of the wind velocity components U and V, potential temperature and vertical heat flux of the boundary layer (solid lines) and profiles used in the model (dashed lines).

S₁ accounts for the fact that both Δθ and θ in the equation describing the pressure distribution associated with the thermal wind (Equation (9)) are functions of n. 52 in the integration of the buoyancy Equation (10) allows for the deviations in the value of θ throughout the layer. However, variations in the ambient density are likely to be relatively small compared to the other terms, and on this basis the value of S ₂ is generally close to unity. S ₃ is associated with the vertical velocity within the boundary layer. Rather than defining a mean value for the layer, the value at H is used in accordance with entrainment theory (Reference Manins and SawfordManins and Sawford, 1979). Entrainment of air from above the boundary layer into the katabatic flow will induce a vertical velocity, W _e, in the ambient air which will ultimately result in a velocity parallel to the isotherms in a stably stratified environment (Reference Manins and SawfordManins and Sawford, 1979). For horizontal isotherms, this small horizontal velocity, U _H is related to W _e by the following equation

(12)

This should be taken into account in the integration of the flow speed U, as follows:

(13)

Unless H is much larger than h U_H << so that the second term on the righthand side of Equation (13) can be neglected, except in the integration of the equation governing the conservation of heat (Equation (5)).

Using these profile factors, and using Equations (8)–(11), Equations (2)–(5) are integrated in the n-direction to height H, and in the y-direction around the circumference of the ice sheet y _i. Replacing (∂θ/∂n) and (∂θ/∂s) with (∂θ/∂z) according to Figure 1, parameterizing surface stress according to −τf/ρ = C _d U ² and neglecting second-order derivatives, gives the following set of equations describing the evolution of the mean layer-velocity components (U and V), the depth of the layer (h) and the temperature deficit (Δθ):

(14)

(15)

(16)

(17)

where R _div is the net radiative divergence across the depth of the boundary layer and f is the Coriolis parameter.

The model is valid for shooting flow which assumes a critical Richardson number (see e.g. Reference BallBall, 1956; Reference Lalaurette and AndréLalaurette and André, 1985), defined by Equation (18) to be less than unity:

(18)

This occurs where the slope is large enough so that the layer velocity exceeds the velocity of small-amplitude gravity waves (N in Equation (1)). Over East Antarctica, this occurs at the top of the glacial slope, around 700 km from the centre of the ice sheet (approximately 300 km from the margin). Initial values of U, V, h and Δθ are therefore chosen to be appropriate to the top of the glacial slope in order to be consistent with this assumption.

2.2. Model parameters and boundary conditions

The profile factors S ₁, S ₂ and S ₃ have been given values of 0.5, 0.9 and 1.0, respectively. These are the same as the values used by Reference GosinkGosink (1989). Further discussion of these values can be found in Reference Ellison and TurnerEllison and Turner (1959), Reference Manins and SawfordManins and Sawford (1979), Reference Lalaurette and AndréLalaurette and André (1985) and Reference HallHall (1992). The cooling rate (R _div/ρC _p) is taken to be 0.03 K m s⁻¹, equivalent to a heat loss from the boundary layer of around 30 W m², and is the same as the value used by Reference GosinkGosink (1989). This value compares well with the observed net longwave flux from the surface of the Antarctic ice sheet (Reference SchwerdtfegerSchwerdtfeger, 1984), which should be close to the net flux divergence, since the net radiative flux at the top of the boundary layer is small.

The drag coefficient in the x-direction, C _du is given a value of 0.0015. This is consistent with previous studies, e.g. Reference Hyde and PeltierInoue (1989), Reference GosinkGosink (1989). A larger value is used in the y-direction and Cdv is taken to be 0.04. This is to take account of other processes which may retard the flow in this direction, which are not included explicitly in the model, namely:

If C _dv is too small, the Coriolis force begins to dominate the evolution of the U component of flow, and the drainage flow becomes supercritical. The geostrophic wind is assumed to be exclusively in the y-direction. This is approximately true of the flow over Antarctica (e.g. Reference Lettau and SchwerdtfegerLettau and Schwerdtfeger, 1967; Reference Schwerdtfeger and MahrtSchwerdtfeger and Mahrt, 1968; Reference Schwerdtfeger and OrvigSchwerdtfeger, 1970), where it is generally observed to be between 2 and 5 ms⁻¹, but these values may be exceeded over Greenland and past ice sheets. In the work presented in this paper a value of 3 m s⁻¹ is used. Higher values of the geostrophic wind produce a deeper and faster boundary layer due to the additional down-slope force (Reference HallHall, 1992). If the geostrophic wind is too small, the flow becomes subcritical where the slope is small because the down-slope buoyancy force is not sufficient alone to maintain drainage flow.

3. Entrainment

Entrainment is parameterized according to the formulation used by Reference Ellison and TurnerEllison and Turner (1959), which was developed following a series of laboratory experiments. It has been used by several authors to describe the entrainment process into the boundary layer over ice (e.g. Reference Manins and SawfordManins and Sawford, 1979; Reference GosinkGosink, 1989). The equation for entrainment is given by

(19)

where E is the entrainment coefficient and A_c and A_k are both constants, with values of 0.0004 and 0.02, respectively, these are the same as those used by Reference GosinkGosink (1989), S₁ is the profile factor described earlier and Ri is the Richardson number for the layer given by

(20)

4. Model Results

Figure 3 shows the evolution of the model boundary layer over two different profiles. The solid line represents a surface profile of Terre Adélie, Antarctica (Reference DrewryDrewry, 1983), and the dotted line is an idealized profile of the form described in Equation (21) with surface slope angle tan α = dH _i/dx

(21)

The model was initialized at the top of the slope (300 km from the margin) with a temperature deficit Δθ of 12 K, consistent with observations (e.g. Reference Parish and WaightParish and Waight, 1987). Initial values of U, V and h, were chosen iteratively to eliminate any abrupt changes in U, V, h or Δθ at the top of the slope. On this basis, for the Terre Adélie profile U ₁ = 10.4 m s⁻¹, V ₁ = 3 m s⁻¹ and h ₁ = 245 m. For the idealized profile U ₁ 12.0 m s⁻¹, V ₁ = 3 m s⁻¹ and h ₁ = 235 m. For the Terre Adélie profile, the resultant wind speed increases from 10.8 m s⁻¹ to around 16 ms⁻¹ at the coast, the largest acceleration occurring in the region of rapidly increasing slopes in the coastal strip. There is a tendency for the wind to be orientated more nearly perpendicular to the contours as the coast is approached, so that at the margin the U and V components of the wind speed are 15.48 and 3.73 m s⁻¹, respectively. The temperature deficit decreases from 12K to around 4 K. Both these features show a generally good agreement with the observations and numerical model results of Reference Parish and WaightParish and Waight (1987). From the observations of Reference Phillpot and ZillmanPhillpot and Zillman (1970), the temperature deficit is between 10 and 15 K on the glacial slopes, decreasing to between 2 and 5 Κ at the coast and the wind-speed observations of Reference Mather and MillerMather and Miller (1967) show an increase from around 10 m s⁻¹, orientated at 45° to the slope at Pionerskaya, to around 18 ms⁻¹, orientated at between 5° and 15° at the coastal sites of Port Martin and Cape Denison. The layer depth of the slab model is slightly less than that which is observed, around 300 m, rather than increasing from 300 to 400 m.

Fig. 3. The evolution of the boundary layer according to the model using a profile taken from Terre Adélie inAntarctica (solid line) and an idealized profile described in the text (dotted line).

The idealized profile shows a reasonable agreement with these results, although the absence of the detailed slope structure, especially in the upper regions of the slope, causes the wind speed to be too high, and the acceleration in the lowest 100 km is less. The broad features of boundary-layer development, namely an increasing wind speed becoming orientated closer to the line of the maximum slope, the slight increase in boundary-layer depth and decrease in temperature deficit are reproduced, giving confidence in the model results. Therefore, the use of idealized profiles in the following section seems reasonable.

4.1. The effect of ice-sheet shape and size

Model runs have been carried out using ice sheets of different size and shape to investigate the boundary-layer evolution as an ice sheet grows. The ice-sheet profile in these cases is always of the form

(22)

where H _i is the height of the ice above mean sea level, h _d is the height of ice at the divide, x is the horizontal distance from the ice divide and s _p is the span of the ice sheet. In the first group of experiments Fig. 4, the ice sheet is allowed to maintain a similar shape, but successively smaller spans are used; for profile 1, h _d = 2.4 km and s _p = 500 km, for profile 2, h _d = 1.75 km and s _p = 360 km, and for profile 3, h _d = 1.25 km and s _p = 215 km. A second group of experiments, shown in Figure 5, was carried out where the maximum height of the ice was changed in each case, while the span s _p was kept constant at 500 km; in profile 4, h _d = 1.75 km, and in profile 5, h _d = 1.25 km. At the top of the slope (100 km from the centre), the profiles are almost identical and the initial boundary-layer depth at the top of the slope, 100 km from the centre, is taken to be 110 m, Δθ ₁ is given a value of 6 K and V ₁ a value of 2 ms⁻¹ in each case. Initial values of U have been chosen to allow smooth boundary-layer development; for profiles 1 to 5, these are 5.4, 5.5, 6.0, 5.0 and 4.6 ms⁻¹, respectively. The results in Figure 4 show that, as the ice sheet shrinks, the temperature deficit, boundary-layer depth and wind speed all change more rapidly, but conditions at the margin of the ice sheet are little changed. The second group of experiments Fig. 5 indicates that a flatter ice sheet causes the layer to slow down, there is less entrainment and the increase in the radiative time-scale relative to the advective time-scale tends to produce a colder boundary layer dominated by radiative rather than turbulent processes. The difference in the temperature deficit between the profiles exceeds 1 K.

Fig. 4. The evolution of the boundary layer for three different idealized profiles with similar shape but different spans (see text for details).

Fig. 5. The evolution of the boundary layer for three different idealized profiles, having the same span but with different slopes (see text for details). Profile 1 is the same as that in Figure 4.

The ablation of ice on an ice sheet depends on the energy budget at the surface as well as the temperature; therefore any model which attempts to assess the potential ablation of ice should include an energy-balance model, but these results show that the evolution of the boundary layer is also important. Investigations by Reference Braithwaite and OlesenBraithwaite and Olesen (1990) suggest that the Greenland ice sheet is becoming steeper as it retreats. If this is the case, the results presented here suggest that entrainment of warm air into the boundary layer contributes significantly to the existing ablation.

5. Conclusions

The experiments presented in this paper show that the evolution of the boundary layer depends on an interplay between ice-sheet slope, entrainment and radiation budget. In particular over the flatter ice sheets, where flow is slower, entrainment is reduced and the radiation balance begins to dominate the heat budget of the boundary layer. In the case where the ice sheet is steeper, the flow is much faster and the additional entrainment introduces a large amount of heat to the surface layer, which would be significant for the ablation of ice. Thus, steeper profiles tend to produce warmer boundary layers. Boundary layers over ice sheets with the same shape but different spans appear to differ only in the horizontal distance over which the change takes place. Of more significance in this case is the temperature deficit, which is likely to be reduced on a smaller ice sheet where the opportunity for radiative cooling over the central plateau is less.

6. Discussion

Experiments with ice-sheet models show the evolution of ice to be sensitive to ablation. Ablation is generally modelled using either energy-balance considerations or more simple degree-day models. This paper has shown that the detailed boundary-layer properties, especially entrainment, are also important, and that an evolving boundary layer over an evolving ice sheet is likely to produce feed-back effects, so that, as the ice profile becomes steeper, there is more ablation. Entrainment is greatest near the edge of the ice sheet, therefore the effect is likely to be more marked at low altitude, steepening the profile further.

The understanding of boundary-layer processes as it stands at present is inadequate. More observations are required, particularly those which enable an assessment of the entrainment of mass and heat into the boundary layer. The role of the geostrophic wind as a significant driving force to the drainage flow must also be considered. Work by Reference Shinn and BarronShinn and Barron (1989) and Reference Kutzbach and WrightKutzbach and Wright (1985) has demonstrated the changing strength of the glacial anticyclone and the mean meridional flow. Experiments with this model confirm the fact that the strength of the geostrophic wind is important for the maintainance of drainage flow, particularly where the slope is small, as has been found in similar studies over Antarctica (e.g. Reference LettauLettau, 1966; Reference Parish and WaightParish and Waight, 1987).

A detailed investigation of the interaction between the glacial anticyclone and the upper mean meridional flow, however, is beyond the scope of a single-layer model such as the one presented here, but should be a consideration of future work. This work has highlighted the need to take account of the interplay between ice-sheet shape and boundary-layer development, particularly entrainment, in estimating ablation of ice. Modelling projects in the future therefore should incorporate these boundary-layer processes. In this respect, the model is of use to both glaciologists and climatologists. It would be easily combined with a snow-melt model by incorporating an energy balance at the surface, and its simplicity means it could provide a parameterization of these important processes in GCMs which cannot explicitly reproduce these processes as they occur on a scale smaller than that of a typical model grid.