Introduction

The glacial geology of Finland offers a very demanding test of our theoretical understanding of ice-sheet flow. This is because many workers (Reference PunkariPunkari, 1980; Reference DongelmansDongelmans, 1995; Reference Kleman, Hättestrand, Borgström and StroevenKleman and others, 1997) have used this geology as evidence for distinct zones of fast ice flow (ice streams) within Scandinavian ice sheet (SIS) during the Younger Dryas period (approximately 11000 years before present or 11kyrBP). This is important because the factors normally associated with fast ice flow (bedrock troughs and/or deformable sediment) are largely absent in Finland. The origin of these geological features is therefore unknown and could potentially further our understanding of the process of ice-streaming.

The principal evidence comes from the regional-scale interpretations of satellite images and air photography. Large lineations (several km in length) often dominate such images and can be used to infer both former ice-flow direction and basal thermal regime. The former information can simply be inferred from the orientation of the streamlined landforms. Information on the latter stems from the assumption that streamlined glacial landforms can only be created where water is present at the bed and there is differential movement between the basal ice and the underlying bedrock (Reference Kleman, Hättestrand, Borgström and StroevenKleman and others, 1997). Areas with abundant lineations are therefore interpreted as zones of former fast ice flow (former ice streams) and areas where lineations are sparse are interpreted as former ice divides or ridges where flow was slow and/or the ice was frozen to its bed. Further information on the relative ages of the lineations, and hence the flow events which created them, can be obtained where lineations cut across one another. In this case, the younger lineation can be assumed to lie over the older one. This type of analysis has allowed glacial geomorphologists to develop detailed palaeoglaciological reconstructions of the former ice streams of the SIS and their interactions (Reference PunkariPunkari, 1984; Reference Kleman, Hättestrand, Borgström and StroevenKleman and others, 1997).

Many of these reconstructions share common elements. The most obvious of which are the fan-like features which trend from northwest to southeast over much of southern Finland. Similar features trend northward in northern Norway; southwards over the Lake Vännern region of Sweden; and eastwards over Lapland and the Kola Peninsula. All are marked by high concentrations of well-aligned lineations and are slightly divergent at their distal end. Their formation is universally attributed to fast ice flow, although Reference Kleman, Hättestrand, Borgström and StroevenKleman and others (1997) describe this as “surging” while Reference DongelmansDongelmans (1995) calls it “streaming”.

If these features are evidence of fast ice flow, then they have two main implications. First, they all occur on the hard-rock Baltic Shield. This runs contrary to the suggestion by Reference Clark and WalderClark and Walder (1994) that fast ice flow is limited to soft, Tertiary sediments. Reference Clark and WalderClark and Walder (1994) base their argument on the theoretical differences in subglacial meltwater drainage over hard rock and soft sediment, which they assume to be the key determinant of fast ice flow. The second implication is that the concentration of ice flow in topographic troughs is often cited as being another key determinant of ice-streaming (Reference GlasserGlasser, 1995). Many of the inferred ice streams mentioned above lie in areas of very low topographic relief (the exceptions being those in the fjord topography of northern Norway). In particular, the inferred ice streams in Finland are separated by regions which are interpreted as having been covered by slow-flowing, cold-based ice. Reference DongelmansDongelmans (1995) and Reference PunkariPunkari (1984) go on to suggest that some lineation patterns reflect the stagnation of one ice stream and later activation of a neighbouring ice stream. The significance of these features, therefore, is that they imply fast ice flow in areas where this is very hard to explain using existing glaciological theory.

This information is particularly interesting to ice-sheet modellers because it represents a coherent picture of ice-sheet dynamics which is relevant at both the regional spatial scales and millennial time-scales over which ice sheets operate. As such, it represents a unique opportunity of testing ice-sheet models and furthering our theoretical understanding of how ice masses flow. This contrasts with data on the present-day ice sheets, which can only reveal information on time-scales up to 200 years at most.

A further advantage in the use of former ice sheets to test ice-sheet models is that the geological record reflects processes occurring at the former ice-sheet bed. In contrast, very few studies have collected data from the bed of a present-day ice sheet because of the logistical challenge of penetrating several kilometers of overlying ice. Geophysical techniques, such as ice-penetrating radar and seismics, have partly addressed this issue; nevertheless data are available only for a very small fraction of the present-day ice sheets.

The major disadvantage of former ice sheets as tests of model behaviour is that the geomorphological record is laid down time-transgressively, which complicates the interpretation of these features dramatically. This initial study will, therefore, concern itself purely with gross features such as the presence or absence of flow concentrations and not the details of their evolution over time.

The present paper is an initial attempt at using this palaeoglaciological record to test a model of ice-sheet flow. Despite the strong case made above for using former ice sheets to test models, there has been little exploration of this in the literature. The few studies that have attempted to model former ice sheets have aimed at the accurate simulation of former ice margins identified in the geological record (Reference Boulton, Payne and DuplessyBoulton and Payne, 1994; Reference Huybrechts and T’siobbelHuybrechts and T’siobbel, 1995). This then allows information on former climate to be inferred.

The focus of the present study is somewhat different. It aims to use the record of former ice-stream dynamics discussed above to test whether a model can generate the types of ice streams thought to have existed over much of Scandinavia during the last glaciation. Previous work (Reference Payne and DongelmansPayne and Dongelmans, 1997; Reference PaynePayne, 1999) has shown that bands of rapid ice flow (oriented parallel to regional ice-surface gradients) can arise spontaneously within an ice sheet. The origin of these ice streams is the interaction between the evolving ice temperature and flow fields. In particular, the “creep instability” feedback identified by Reference Clarke, Nitsan and PatersonClarke and others (1977) appears to control the onset of fast ice flow within these models. This instability links temperature and flow in two ways. First, ice temperature is known to play a crucial role in determining ice viscosity (Reference PatersonPaterson, 1981) and therefore the rate at which ice deforms. This effect is thought to have an exponential form so that ice near pressure melting point deforms approximately 1000 times faster than ice at 230 K (Reference Paterson and BuddPaterson and Budd, 1982). The second link between flow and temperature is that the flow of ice generates heat by the process of dissipation. Faster-flowing ice therefore generates much more heat than slower-flowing ice. These two effects are linked so that initially small perturbations in flow could develop into ice streams. Reference PaynePayne (1999) uses these feedbacks to investigate the dynamics of ice-streaming in West Antarctica, while Reference Payne and DongelmansPayne and Dongelmans (1997) find they can also lead to the formation of ice streams in the absence of topographic variation. This study forms a logical progression in that it attempts to simulate the ice streams of the former SIS, which (as noted above) cannot readily be explained by existing theories of rapid ice flow.

Previous thermomechanical modelling work on the SIS has produced patterns of ice flow, temperature and thickness which are approximately concentric (Reference Boulton, Payne and DuplessyBoulton and Payne, 1994; Reference Huybrechts and T’siobbelHuybrechts and T’siobbel, 1995). The only exception to this is where the topography of the coastal mountains and fjords produce local concentrations of flow. As such, these results are incompatible with interpretations based on the area’s glacial geomorphology.

Model Description

This section summarizes the thermomechanical ice-sheet model. Reference Payne and DongelmansPayne and Dongelmans (1997) provide a more detailed discussion. The model is time (t) dependent and is three-dimensional (with horizontal dimensions x and y, and vertical dimension z which is positive upwards from sea level).

The evolution of ice thickness (H) is given by the continuity equation

where ∇. is the two-dimensional horizontal divergence operator, b is the local snow accumulation/ablation rate, and U is the vertically averaged, horizontal velocity vector.

Ice deformation is assumed to be driven solely by horizontal shear stresses (τ _xz and τ _yz ) which, at the spatial scales over which ice sheets operate, can be approximated by

where s is ice-surface elevation, ρ is ice density and g is acceleration due to gravity.

Assuming a non-linear flow law for ice (Reference GlenGlen, 1955), one can find horizontal velocity as

where h is bedrock elevation, n is usually taken as 3 and the parameter A introduces the temperature dependence of ice deformation. Similar expressions are used to find τ _yz and u_y (z), and vertical velocity is found from the horizontal velocity field using the incompressibility condition. Sliding is not included in the work reported in this paper, so the only mechanism of flow is internal deformation.

The temperature dependence of ice deformation is incorporated using an Arrhenius relation, choosing constants based on work by Reference Paterson and BuddPaterson and Budd (1982):

where T* is absolute temperature corrected for the dependence of melting point on pressure, f is used to tune the model, a is a constant of proportionality, Q is the activation energy for ice creep, and R is the universal gas constant. The values used in this study for these, and all other constants, are given in Table 1.

The ice temperature (T) evolves according to

where k is ice conductivity and c is its specific heat capacity. These and the other thermal parameters used in the model are assumed to be independent of ice temperature. The terms in Equation (5) represent vertical diffusion (horizontal diffusion is assumed negligible following Reference PatersonPaterson (1981)), horizontal advection, vertical advection and dissipational heat generation, respectively.

Two boundary conditions are required for the integration of Equation (5) forward through time. The boundary condition at the ice-sheet base depends on whether the basal ice is at pressure melting point or not:

where G is the geothermal heat flux (assumed constant). The first condition applies to basal ice maintained at pressure melting point, when the heat generated at the bed (geothermal and sliding) exceeds the heat flux away from the bed into the ice. The second condition applies when the basal ice is frozen. The upper boundary condition is a specified ice-sheet surface temperature (T _a), the details of which are discussed below.

When basal temperatures reach pressure melting point, a melt rate is calculated from the difference between the heat supplied to the ice-sheet base (the geothermal heat flux) and the heat diffused up through the ice mass. The latter depends on the basal gradient of the temperatures found from the evolution equation (Equation (5)).

Isostatic depression and rebound are included in a rudimentary fashion using

where Θ is a time constant controlling the rate of deflection, ρ_a is the density of the asthenosphere and h ₀ is the isostatically relaxed topography (taken as the present-day one).

The model summarized above is solved using finite differences on a regular horizontal and a stretched vertical grid using mostly implicit methods. All numerical experiments described here were performed using a 20 km horizontal grid- and time-steps of two years.

Experimental Design

Before the model described in the previous section can be numerically integrated forward through time, atmospheric and topographic boundary conditions must be specified as well as initial conditions.

The gridded bedrock topography was constructed using the Terrain Base global, 5 minute dataset which contains data both above and below sea level (Row and Hastings; ftp://ftp.ngdc.noaa.gov/Solid_Earth/Topography/tbase.5min). The map projection chosen for this study was the Lambert Conic Conformal with standard latitudes of 52 and 72° N, and centred on 60° N and 15° E.

The present-day topography is used as an initial condition in the model and is assumed to be fully relaxed isostatically. Figure 1 shows this topography.

Fig. 1. The present-day topography of the Fennoscandian land mass (contours are solid above sea level and dashed below). The units of the x and y axis are kilometers.

The mean annual air temperature required as an upper boundary condition for Equation (5) is based on sea-level air temperature T _a* from a global circulation model simulation of the 11 kyr BP climate (Kutzbach and others; ftp://ftp.ngdc.noaa.gov/paleo/gmcoutput/ccml).

Where μ is the environmental lapse rate.

The parameterization of the net mass balance is secondary to the main aim of this study, the main criteria being that the scheme employed should allow the growth of an ice sheet that is comparable to the ice sheet that developed over Scandinavia during the Last Glacial Maximum.

A simple, but very effective parameterization is that used by Reference OerlemansOerlemans (1981). The local snow accumulation/ablation rate b is parameterized in terms of height above an equilibrium plane altitude, s_* , by

where, following Oerlemans, a = 0.732 × 10^-3m yr^-1 and d = 0.268x10^-6m^-1yr^-1. This gives a parabolic form for mass balance that reaches a limiting value of 0.5 m yr^-1 at 1500 m above the equilibrium plane altitude (EPA).

It remains to specify the EPA and this is prescribed as a function of geographical position (latitude, ?; longitude, λ) in the following manner:

Thus EPA depends on a meridional gradient (α) centred on some reference latitude (ϕ ₀), on a zonal gradient (γ) and on a reference elevation (E ₀). Physically, the meridional gradient accounts for variations in solar insolation, while the zonal gradient accounts for the higher precipitation rates found in the western part of Europe. The climate at the start of glaciation is assumed to be similar to that of the present-day.

While there is a physical basis for the parameterization employed here, the choice of values is governed by our need to have an ice sheet that will first develop and then grow to a size and extent comparable to the SIS. For our purposes, values of E ₀ = 890 m, α =70 m/° N, ϕ ₀ = 60° N and γ = 25 m/° E are used.

The model does not include the complicated physics of grounding, ice-shelf flow and iceberg calving. This is principally because the focus of this study is on thermomechanical interactions within the grounded ice mass. A simple means of incorporating these effects is to artificially lose all ice that flows past a 300 m bathymetric limit. This effectively constrains the evolving ice sheet to the Fennoscandian land mass. We therefore deliberately ignore the possibility of ice growth in Svalbard which could grow to cover the Barents Sea area. This choice was made to reduce model computation because the main aim of the work is to understand the features generated in Scandinavia.

Results

The model experiment described in this paper started with a zero ice initial condition and lasted for 150 kyr. Figure 2 shows the evolution of ice-sheet volume and area during this run. The mass balance forcing (Equations (9) and (10)) was constant during the time although there will be feedback between increased ice-surface elevation and the local snow-accumulation rate (6). The first 25 kyr of the experiment is characterized by slow growth (both in area and volume) as the ice sheet expands out of the Scandinavian mountains. A 75 kyr period of rapid growth then follows, during which the modelled ice sheet reaches a maximum extent similar to that reconstructed from the area’s glacial geomorphology (Reference Boulton, Smith, Jones and NewsomeBoulton and others, 1985). This maximum then remains in equilibrium over the course of the last 50 kyr of the experiment. The proportion of the ice sheet that experiences basal pressure melt lags overall growth slightly and reaches equilibrium at approximately 40 %.

Fig. 2. Time series of modelled (a) volume and (b) area for the experiment described in the text. The dashed curve in the area plot refers to area at pressure melting point, while the solid curve is total area.

Various features of the 150 kyr, equilibrium ice sheet are shown in Figure 3. The predicted ice-surface elevation distribution is dominated by a central dome over southern Finland, accompanied by ridges running toward the Jotunheimen and Finnmark areas. By 150 kyr, basal ice temperature is in equilibrium with this pattern of ice-surface elevation. Its distribution shows cold ice under the domes and ridges described above, which is in accord with previous studies (Reference Boulton, Payne and DuplessyBoulton and Payne, 1994; Reference Huybrechts and T’siobbelHuybrechts and T’siobbel, 1995). However, the basal temperature distribution away from these areas is dominated by a series of wedges of warm ice extending for distances up to 600 km inward of the margin. These wedges are separated by the colder ice, more typical of the dome and ridge areas. The wedges are also evident in the distributions of both the vertically integrated flow coefficient (A in Equation (3), which is linked to ice temperature via Equation (4), and the vertically averaged horizontal ice velocity.

Fig. 3. Model output at 150 kyr for (a) ice-surface elevation in m, (b) basal ice temperature in K, (c) vertically integrated ice-flow coefficient A in 10^-16 × s^-1 Pa^-3 and (d) vertically averaged horizontal velocity in myr^-1. The arrows in (d) are scaled to the maximum velocity of 300 myr^-1. The present-day coastline is shown to aid interpretation.

Figure 4 shows the distribution of vertically averaged horizontal velocity at various times during the growth of the ice sheet. At 30 kyr, the ice sheet is still restricted to upland areas. There is little pattern to the distribution of flow, other than an increase away from the divide and towards the steeper topographic slopes of the coast. By 50 kyr, the ice sheet has expanded through southern Finland. Flow along the Norwegian coast has become concentrated in several locations, however flow along the eastern margin is still relatively uniform. The ice sheet has expanded into Poland by 80 kyr and the eastern margin has now developed wedges of fast ice flow separated by relatively slow-flowing areas. Comparison with the pattern at 100 kyr indicates that once a fast-flow feature has formed, it is carried forward with the advancing position and does not change its relative position along the margin significantly. The pattern of flow at 100 kyr is very similar to that at 150 kyr (Fig. 3), indicating that it has reached a steady-state.

Fig. 4. Model output for vertically averaged horizontal velocity in myr^-1 at (a) 30, (b) 50, (c) 80 and (d) 100 kyr into the experiment. The arrows are scaled to a maximum velocity of 300 myr^-1.

A second experiment was performed using exactly the same design outlined above, except that a global circulation model simulation of the 21kyrBP (Last Glacial Maximum) climate was used as the upper boundary condition for Equation (5). Mean annual air temperatures in this dataset are typically 15 K colder than those of the 11kyr BP dataset (used previously).

The equilibrium ice sheet produced in this second experiment is both cooler and larger than that of the initial experiment. The increase in size arises from a coupling between cooler, slower-flowing ice, increased ice-surface elevations (thicker ice) and enhanced accumulation (Equation (9)). The same development of warm ice wedges occurs, however individual wedges tend to extend further back into the ice sheet and are narrower than those shown in Figure 3. A similar effect is noted by Payne and others (submitted).

Discussion

The preliminary results described above indicate that localized creep instability could potentially account for the fanlike landform assemblages reported by Reference DongelmansDongelmans (1995) and Reference Kleman, Hättestrand, Borgström and StroevenKleman and others (1997). The formation of these flow concentrations or ice streams within the model occurs as a consequence of ice-flow—temperature coupling. Initially, small differences in ice flow can become magnified because fast-flowing ice generates more heat via dissipation; becomes warmer and less viscous and, consequently, flows even faster. The surface of the ice sheet becomes drawn down into these faster-flowing areas. This leads to a further concentration of ice flow into these areas because ice always flows locally downslope. In contrast, the developing ridge areas are robbed of ice and experience the above feedback operating in reverse. Ridge ice therefore becomes progressively colder and flows more slowly than the neighbouring stream ice. This thermomechanical mechanism does not require the presence of either topographic troughs or soft sediment to account for local concentrations of ice flow. Indeed, the present model can generate this behaviour without the inclusion of any form of basal slip.

The widths of the flow concentrations generated by the model vary from 100–200 km, which is in good agreement with the widths of the fan-like landform assemblages identified by Reference DongelmansDongelmans (1995). The distal divergence of these landform assemblages is also reproduced by the model.

The widths cited above mean that the flow concentrations are resolved by the model fairly well (typically by 5–10 grid-cells across the feature). However, it is clear the numerical grid does play some role in controlling the precise form of the features. In particular, the ice streams flowing directly northwards over the Kola Peninsula and northern Sweden appear to be controlled by the numerical grid.

The flow concentrations have a wedge-like form. This contrasts with the features observed by Reference Payne and DongelmansPayne and Dongelmans (1997), which had lower width-to-length ratios and were more linear. Two main factors could account for this difference. First, the present experiment has a free, ablating margin which is convex in plan form. This contrasts with the ice margin in the Reference Payne and DongelmansPayne and Dongelmans (1997) experiments, which coincided with the straight edge of the model domain and was not ablating. Ice flow near the margin in the present experiments is therefore compressive as opposed to extensional, and its precise direction is also less-constrained. Second, the present experiments do not incorporate sliding, whereas those of Reference Payne and DongelmansPayne and Dongelmans (1997) did. The onset of sliding represents a discontinuity in the model, whereby adjacent gridpoints can have very different velocities. This strengthens the feedback between ice flow and temperature, and results in more tightly defined concentrations of flow.

Another contrast with the work of Reference Payne and DongelmansPayne and Dongelmans (1997) is that the flow concentrations in some of their experiments oscillated. The flow concentrations in the present experiments are steady-state features. The cause of this difference could again lie with the inclusion of sliding. The discontinuity associated with the onset of sliding (outlined above) leads to dramatic increases in ice flow, which has the effect of capturing flow from adjoining flow concentrations. The more gradual nature of flow concentration in the present study appears to lessen this competition.

This work is not an attempt to accurately model the overall evolution of the SIS during the Weichselian. We have deliberately kept the simulation of climate and ice-sheet mass balance simple. The parameterization we employ does, however, generate an ice mass which is similar in size and shape to that reconstructed from the area’s glacial geomorphology. In addition, the time required for this ice sheet to grow to its maximum seems realistic and would certainly fit within a typical 100 kyr glacial cycle.

Further work should investigate the stability of individual ice streams during the growth and decay of the ice sheet, and the effect of the various parameters employed within the model on overall ice-stream formation. Once this has been done, more detailed ablation and accumulation parameterizations could be employed to generate a more realistic simulation of SIS evolution during the Weichselian.