The ice-sheet model

The model follows those of Mahaffy (Reference Mahaffy1976) and Budd and Smith (Reference Budd and Smith1981), and is based on the continuity equation for ice thickness. Changes in ice thickness are tested through time using the equation

where Ζ is ice thickness, t is time, b is net mass balance, U _d is vertically averaged deformation velocity and U _s is sliding velocity. The flow section of the model is detailed in Payne and others (Reference Payne, Sugden and Clapperton1989) and is driven by the balance between basal shear stress and surface slope driving stress. Sliding is found from an empirical relationship between basal shear stresses and normal stress on the ice-sheet bed. Equation (1) is approximated over a regular rectangular grid using finite differences and can be used to simulate the time dependent behaviour of the ice sheet.

This paper is mostly concerned with the net mass-balance term (b) in Equation (1), which is defined as:

where Ρ is the precipitation rate, Ab is the ablation rate and Eff is the fraction of annual precipitation falling as snow (effectiveness). Annual precipitation rate provides the link between the ice-sheet model and the atmospheric moisture model and is discussed later.

The method used to estimate net mass balance first involves the specification of surface air temperatures. For each year of the model run the annual temperature cycle is approximated using a sinusoidal curve, the amplitude and mean value of which are estimated from latitude, surface elevation and a specified temperature change (Payne, Reference Payne1988). Monthly air temperatures are obtained from this curve and these are then used to estimate monthly effectiveness and ablation values (Table I). In the former case we use the relationship of Lauscher (Reference Lauscher1976), which is based on western European data. In the latter case, the monthly air temperature is converted to the monthly sum of day-degrees above 0°C and this is then used to calculate the ablation rate. This latter procedure is based on the West Greenland work of Braithwaite and Olesen (Reference Braithwaite and Olesen1984).

The ice-sheet model requires three sets of input data. The first is bedrock topography, for which a rectangular grid of 5 km by 5 km contains the maximum altitude in each grid element (Fig. 1). The second data set consists of similar grids containing annual precipitation rates at the surface of the ice sheet. These are obtained from the output of the atmospheric moisture model at 500 year intervals. The third data set is a time series of regional July air temperatures for sea level, which drives the ice-sheet model. These data are shown in the upper graph of Figure 2 in the form of differences from the present July temperature of 17°C. These temperature estimates are obtained from the ecology of fossil Coleoptera assemblages (Coope, Reference Coope, Wright and Moseley1975). They partially determine the modelled values of the ablation rate and the precipitation effectiveness, and are fed directly into the equation modelling the annual temperature cycle.

Table I.

The equation used to estimate monthly effectiveness and ablation values.

Fig. 1.

Contoured representation of Scottish bedrock topography, based on the highest point in each 5 km by 5 km grid square.

Fig. 2.

Ice-sheet model input and output. The upper graph shows air temperature and predicted mean precipitation rate; the lower curve shows predicted ice-sheet volume.

The atmospheric moisture model

A separate numerical model is used to assess the effect of falling North Atlantic sea-surface temperatures on the distribution of precipitation over Scotland. The model is based on the vertically integrated continuity equation for atmospheric moisture (Reference Sanberg and OerlemansSanberg and Oerlemans, 1983):

where W is the amount of water vapour in a vertical column extending from the ground surface to the “top” of the atmosphere; W _sat is the amount of moisture held in this column when the air is saturated; T* is the time scale over which evaporation occurs; S is the upwind surface slope: f ₀ and f₁ are empirically derived parameters that determine the relationship between surface slope, atmospheric moisture and precipitation; U is a specified wind velocity field; and D _w is the diffusivity of atmospheric moisture.

The sophistication of the model is intermediate between a purely empirical approach and a full global circulation model. The equation contains terms dealing with evaporation, precipitation, advection and diffusion. The first term in Equation (9) concerns evaporation. The saturated moisture content (W _sat) is determined from surface elevation at a point and from the equivalent sea-level air temperature at that point. It is the variation of sea-level air temperature which drives the atmospheric moisture model. Precipitation, represented by the second term in Equation (9), is particularly important because it is the link between the atmospheric moisture model and the ice-sheet model. Precipitation responds to changes in both atmospheric moisture and upwind surface slope. Upwind surface slope is specified and provides the basis for introducing the effects of both orographic enhancement and the “precipitation shadow”. Predicted precipitation distribution therefore varies in response to the changing topography of the modelled Loch Lomond ice sheet. The third and fourth terms in Equation (9) represent advection and diffusion respectively, and are both concerned with the horizontal movement of moisture through the atmosphere.

The atmospheric moisture model is run at two resolutions. A low-resolution flowline is used to model the flow of atmospheric moisture across the North Atlantic towards Scotland. The horizontal moisture fluxes over Scotland are then fed into the high-resolution two-dimensional rectangular grid of the ice-sheet model. In this way, gross North Atlantic moisture fluxes are found quickly using the low-resolution model and the pattern of precipitation over Scotland is found in detail using the 5 km by 5 km high-resolution model.

Fig. 3.

The transects along which data was collected for the atmospheric moisture model. Also shown are the positions of the North Atlantic polar front 10 000–8 000 a B.P. After Ruddiman and Mclntyre (1981).

The response time of the atmospheric moisture system is considerably shorter than that of the ice-sheet system. The atmospheric moisture model is therefore run to steady-state by iterating Equation 9 forward through time until the volume of atmospheric moisture becomes constant. This process is repeated at regular 500 year intervals using the appropriate North Atlantic sea-surface temperatures to determine the modelled saturation moisture content and Loch Lomond ice-sheet thicknesses to determine the orographic component of precipitation. The final steady-state precipitation distribution is then fed back into the ice-sheet model, which is run forward for a further 500 years. This procedure is used to reduce excessive computation time and yet does not greatly affect the behaviour of the overall model.

The primary input to the North Atlantic flowline component of the model are sea-surface temperature transects across the ocean. These are derived from the maps of Ruddiman and Mclntyre (Reference Ruddiman and McIntyre1981) and Mclntyre and others (Reference McIntyre1976) based on the palaeoecology of deep sea sediments. The data set encompasses the migration of the North Atlantic polar front from a position between Newfoundland and Iceland to its maximum southerly extent of 45°N 10 000 years B.P. and its subsequent retreat (Fig. 3). Sea-surface temperatures, which fall to a minimum of 11°C below present-day values, are converted to the sea-level air temperatures needed by the model using the regression equation:

where T _a is sea-level air temperature and T _w is sea-surface temperature, both in °C.

The main inputs to the high resolution Scottish component of the model are the horizontal moisture fluxes from the North Atlantic model and the Scottish topography data used in the ice-sheet model. This latter has two components: the unchanging bedrock topography and the fluctuating form of the modelled ice sheet.

To investigate the effects of changing wind directions during the Loch Lomond stadial two alternative data sets are used to represent westerly and south-westerly airflow conditions. These alternative wind directions have two main consequences in the model. First, the fetch of the airflow alters its geometry relative to the North Atlantic polar front (Fig. 3). Secondly, the distribution of upwind ground-surface slopes obtained from the Scottish topography alters.

Testing the atmospheric moisture model

The atmospheric moisture model was run to steady-state using present-day North Atlantic sea-level air temperatures and a westerly wind direction. The purpose was to test the reliability of the model. The predicted precipitation distribution is shown in Figure 4. It compares favourably with the actual present-day pattern. Both have peak precipitation rates of 2400 mm a⁻¹ and typical sea-level values of 750–1000 mm a⁻¹. The high rates to the west are associated with ground-surface slopes of 0.005– 0.007, while the reduced rates in the east are caused by moisture depletion in the precipitation shadow of the uplands. Atmospheric moisture content (W) in the east is typically one-third of its value in the west. The similarity of predicted and actual patterns gives some credence to the low-resolution flowline model.

Fig. 4.

Contoured representation of (a) present-day precipitation as modelled with a westerly wind direction compared with (b) the actual precipitation pattern.

The Loch Lomond ice sheet

The coupled model was run for the period 12 500– 7 500 a B.P. using the July air temperature input shown in Figure 2 (upper graph) and a westerly wind direction. Predicted ice-sheet volumes are shown in the lower graph of Figure 2, while the predicted precipitation rates are shown as bars in the upper graph of Figure 2. The predicted precipitation pattern for 10 000 a B.P. is shown in Figure 5. Precipitation rates are typically half of their present-day values. However, high rates still occur to the west, where general atmospheric moisture reduction is countered by increased orographic precipitation caused by the elevated Loch Lomond ice-sheet surface.

Fig. 5.

Contoured representation of precipitation at 10 000 a B.P. as modelled with a westerly wind direction.

Fig. 6.

Predicted ice extent and thickness at the maximum of the Loch Lomond glaciation 10 000 a B.P.

The ice sheet reaches a maximum at 10 000 a B.P. and its dimensions are shown in Figure 6. The maximum represents a fleeting transition between phases of rapid growth and decay. These rapid changes in volume reflect both changing air temperatures and the positive feedback between ice-surface elevation and net mass balance.

There is generally good agreement between the modelled maximum ice-sheet extent and reconstructions based on geomorphological evidence. The ice sheet builds up in the western Highlands and is centred over the topographic basin of Rannoch Moor. Comparison with field-based reconstructions such as that of Sissons (Reference Sissons1976) shows that the fit is best in the south and worst in the northwest Highlands, where the model overpredicts the extent of ice.

Alternative wind directions during the Loch Lomond stadial

AU of the results described above were obtained using a model incorporating westerly airflow. This section explores the results of model simulations using a south-westerly air flow.

The pattern of precipitation 10 000 a B.P. is similar to that for a westerly wind direction (Fig. 5). However, there are two differences. First, precipitation rates are generally higher owing to the increased atmospheric moisture fluxes advected from latitudes south of the North Atlantic polar front (Fig. 3). Secondly, the area of peak precipitation is further south, in the southern Grampians and away from the north-west Highlands and towards the eastern Grampians.

These changes in the precipitation pattern have two effects on the development of the ice sheet. First, the generally higher precipitation rates cause a faster rate of growth and a larger maximum ice sheet. Secondly, the southern Grampians play a more important role in ice-sheet growth.

This latter effect is interesting because it implies that the southern Grampians are more favourable for ice-sheet growth than the north-west Highlands when the airflow is south-westerly. One of the highlighted mismatches between the model using a westerly air flow and geomorphological reconstructions of the Loch Lomond ice-sheet maximum was the overprediction of ice in the north-west Highlands. This mismatch argues in favour of a south-westerly airflow during the Loch Lomond glaciation.