1. Introduction
During the Weichselian glacial period, large parts of northern Eurasia have been covered by ice sheets. The Quaternary Environment of the Eurasian North (QUEEN) project has compiled the most up-to-date reconstruction of ice-sheet extent in Eurasia (Reference SvendsenSvendsen and others, 1999; Reference Thiede, Bauch, Hjort and MangerudThiede and others, 2001). This reconstruction differs from earlier compilations especially in the southern Kara Sea and the Taymyr Peninsula area, which are deduced to have stayed ice-free during the Last Glacial Maximum (LGM) around 20 kyr BP (Reference SvendsenSvendsen and others, 1999; Reference Polyak, Levitan, Gataullin, Khusid, Mikhailov and MukhinaPolyak and others, 2000; Reference Mangerud, V. and SvendsenMangerud and others, 2002). Details of this extent are still debated: Reference Alexanderson, Hjort, Moller, O. and PavlovAlexanderson and others (2001) conclude that there is evidence of LG M ice in northwest Taymyr, whereas the QUEEN results indicate arid and ice-free conditions. The QUEEN ice-sheet extent is depicted in Figure 1 together with the location of major palaeo-ice streams.
Fig. 1. Ice-sheet extent at LGM (solid line) according to QUEEN (Reference SvendsenSvendsen and others, 1999; Reference Thiede, Bauch, Hjort and MangerudThiede and others, 2001). Added to this are the approximate positions for the ice streams identified by Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others (2001): the Norwegian Channel ice stream (A), the Baltic Ice Stream (B), the Danish sub-lobe (B1), the Oder sub-lobe (B2), the Wisla sub-lobe (B3), the Kurshkiy sub-lobe (B4), the Riga sub-lobe (B5), the Finnish ice stream (D), the Karelian ice stream (E) and the Novgorod ice stream (F).
Proxies for local temperature and precipitation during the LGM are sparse. In addition to more global data from Greenland ice cores and deep-sea sediment records, some local proxies based on palaeobotanical evidence are available (Reference TarasovTarasov and others, 1999). These authors conclude that the temperature anomaly (LGM minus present) for the coldest month of the year at 18 kyr BP is largest in areas near Novgorod, in the range –29 to –20°C. From there, the temperature anomaly decreases to the east. In the Taymyr Peninsula and the Kara Sea areas it is –10 to –3°C. The temperature anomaly for the warmest month is approximately –3 to 0°C in the whole domain. The anomaly in precipitation rate is –450 to –200 mm a–1.
If these reconstructions are compared to Palaeoclimate Modelling Intercomparison Project (PMIP) results from the United Kingdom Meteorological Office (UKMO) (Reference Hewitt and MitchellHewitt and Mitchell, 1997), we can see that the winter temperature anomaly in southeast regions is in agreement with the results obtained by Reference TarasovTarasov and others (1999). However, the UKMO results show an anomaly of similar magnitude in the northeast, whereas Reference TarasovTarasov and others (1999) state that the LGM climate was warmer in this region.
Simulations of Fennoscandian and related ice sheets have been conducted by Reference Fastook and HolmlundFastook and Holmlund (1994), Reference PeltierPeltier (1994), Reference Huybrechts and T’siobbelHuybrechts and T’siobbel (1995), Reference Payne and BaldwinPayne and Baldwin (1999) and Reference Siegert, Dowdeswell and MellesSiegert and others (1999). Simulations of Northern Hemisphere ice-sheet retreat utilizing atmospheric general circulation model (AGCM) results have been carried out by Reference Charbit, Ritz and RamsteinCharbit and others (2002).
Reference Siegert, Dowdeswell and MellesSiegert and others (1999) have utilized the QUEEN project data for ice-sheet extent in their simulations. In their simulation, a zero level of precipitation is set across the Taymyr Peninsula, and resulting reconstructions are then compared with geological evidence. In this way, the most likely glacial scenario is selected.
According to Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others (2001), high flow rates are needed in the Baltic basin to achieve the inferred LGM extent. Several major ice streams are inferred in the area (Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others, 2001; Reference Stokes and ClarkStokes and Reference Clark and MixClark, 2001) (Fig. 1). The Norwegian Channel Ice Stream (NCIS) has recently been identified from geological data (e.g. glacially fed fans) (Reference Sejrup, Landvik, Larsen, Janocks, Eiriksson and KingSejrup and others, 1998, Reference Sejrup, Larsen, Landvik, King, Haflidason and Nesje2000), the Baltic Ice Stream (BIS) and its sub-lobes have been widely discussed (e.g. Reference Kleman, Hättestrand, Borgstrom and StroevenKleman and others, 1997) and the Finnish, the Karelian and the Novgorod ice streams are discussed by, for example, Reference Houmark-Nielsen and HumlumHoumark-Nielsen and Humlum (1994), Reference LarsenLarsen and others (1999) and Reference Lunkka, Saarnisto, Gey, Demidov and KiselovaLunkka and others (2001).
2. Modelling Approach
The dynamics and thermodynamics of the Fennoscandian and adjacent glaciation is investigated with the ice-sheet model SICOPOLIS (SImulation COde for POLythermal Ice Sheets; Reference GreveGreve, 1997a, Reference Greveb). The model describes the material ice as an incompressible, heat-conducting, power-law fluid with thermomechanical coupling due to the strong temperature dependence of the flow law,
where D = sym grad u (velocity v) is the strain-rate tensor, tD the Cauchy stress deviator, 
the effective shear stress, n the power-law exponent and A{T'; the flow-rate factor, which depends exponentially on the temperature T relative to the pressure-melting point and linearly on the water content ω (see Reference Greve, Weis and HutterGreve and others, 1998). The flow-enhancement factor E is equal to unity for pure ice and can deviate from unity due to the softening or stiffening effect of impurities in the ice.
Cold ice with a temperature below the pressure-melting point is distinguished from temperate ice with a temperature at the pressure-melting point, the latter being considered as a binary mixture of ice and small amounts of water. The interface that separates cold and temperate ice is monitored using Stefan-type energy-flux and mass-flux matching conditions. Basal sliding is assumed to be zero for a cold base, and for a temperate base it is described by the Weertman-type sliding law in the form of Reference Greve, Weis and HutterGreve and others (1998),
where usl is the basal sliding velocity, C sl the sliding coefficient, H the ice thickness and h the surface elevation.
The model computes the three-dimensional temporal evolution of the ice extent, thickness, velocity, temperature, water content and age in response to external forcing. The latter must be specified by (i) the mean annual air temperature at the ice surface, (ii) the surface mass balance, which is ice accumulation (here assumed to be snowfall) minus ablation (here assumed to be melting), (iii) the global sea level and (iv) the geothermal heat flux entering the ice mass from below. Isostatic depression and rebound of the lithosphere due to changing ice loads is described by a local-lithosphere-relaxing-asthenosphere (LLRA) model with the isostatic time lag τiso (Reference Le Meur and HuybrechtsLe Meur and Huybrechts, 1996; Reference Greve, Straughan, Greve, Ehrentraut and WangGreve, 2001). A more detailed description of the model is given by Reference GreveGreve (1997a, Reference Greveb).
Since SICOPOLIS is based on the shallow-ice approximation (SIA), it neglects normal stress deviators and shear stresses in vertical planes. The large-scale behaviour of ice sheets is simulated well with the SIA; however, it is not valid locally in the vicinity of ice domes and close to the margin. Further, fast-flowing ice streams with reduced basal drag are not adequately described by the SIA. Nevertheless, SICOPOLIS can produce fast-flow features which arise from sliding over a temperate ice base (Equation (2)). These fast-flow features can resemble real ice streams, but as no specific ice-stream dynamics is included, they are a coarse approximation of the real dynamics.
Ice shelves are not treated explicitly. Instead, the ice sheet is allowed to glaciate the continental shelf below sea level, provided the seabed elevation is larger than a threshold value zthresh. This parameter was introduced in order to prevent the ice sheet from glaciating the deep sea, which is of course not realistic.
The areal coverage in our simulations is from the North Atlantic to the Taymyr Peninsula, and covers an area of 3680 6 3800 km2. This area was chosen such that it covers the reconstructed ice-sheet extent in the east. Any artificial ice-sheet growth further to the east is finally cut off by the margin of the domain. A clear shortcoming of this set-up is that the British Isles are not yet included; this will be done for future simulations. The stereographic projection with a reference latitude of 71 ° is applied. The grid spacing for the simulations in horizontal direction is 40 km, leading to a total of 93 6 96 gridpoints. The vertical resolution is 21 grid-points in the cold-ice region, 11 gridpoints in the temperate-ice region (if existing) and 11 gridpoints in the lithosphere (Reference GreveGreve, 1997b).
Simulation time is from 250 kyr BP to the present, of which the first approximately 50 kyr are required for spin-up (that is, for losing the memory of the arbitrarily chosen ice-free initial conditions). Even as our main focus is on the reconstruction of the LGM around 20 kyr BP, the dynamics of the whole Weichselian glaciation is of interest. Also, the ice-load history is important for the bedrock elevation due to the isostatic time lag. Values of physical parameters are listed in Table 1.
Table 1. Physical parameters of the ice-sheet model
3. Climatic Forcing
In order to provide modern surface temperature and precipitation values, measured and modelled datasets by Reference Legates and WillmottLegates and Willmott (1990) are used. These data have been gathered from the Global Historical Climatology Network (Reference VoseVose and others, 1992) and station records for the years 1950–96. The source resolution is 0.5° x 0.5°, and the data are interpolated for each season with inverse-distance weighing to the 40 km grid of SICOPOLIS.
Our first attempt to parameterize past climate conditions was similar to the method applied by Reference Greve, Wyrwoll and EisenhauerGreve and others (1999). Past temperatures were calculated from present values with a conversion formula that takes into account that climate is generally warmer or colder, and that the surface elevation has changed. However, this approach leads inevitably to a very extensive glaciation because it does not account for different spatial patterns of LGM and present climate conditions, and is therefore not suitable for a detailed study of Eurasian glaciation.
In the current approach, values for temperature and precipitation from the UKMO21 LGM simulations (Reference Hewitt and MitchellHewitt and Mitchell, 1997) are utilized. The UKMO21 dataset was chosen from comparisons of ice-sheet mass budgets in the PMIP simulations (Reference PollardPollard and others, 2000). The coarse atmospheric model output was again interpolated with inverse-distance weighting to the model grid.
The seasonal temperature and precipitation distributions are then interpolated linearly between the present (Reference Legates and WillmottLegates and Willmott, 1990) and LGM (UKMO21) values with a glaciation index g(t). This index scales the Greenland Icecore Project (GRIP) δ18O record (Reference DansgaardDansgaard and others, 1993) to represent LGM (g = 1) and present (g = 0) conditions (Fig. 2). The interpolation is carried out by
Fig. 2. Glacial index g(t) as a function of time. This index scales the GRIP δ18O record (Reference DansgaardDansgaard and others, 1993) to represent glacial (g =1) and present (g =0) conditions.
where x and y are Cartesian coordinates which span the stereographic plane, t is the time, f is the surface temperature or precipitation, and fpresent(x; y) and fLGM(x; y) are the spatially resolved data for the two time-slices, respectively. Conversion from seasonal precipitation P to seasonal snowfall (solid precipitation) S is done with the empirical relation by Reference MarsiatMarsiat (1994),
where Ts is the seasonal surface temperature. Mean annual air temperature is the mean of seasonal temperatures.
Climatic forcing is completed by the degree-day method for surface melting in the form presented by Reference ReehReeh (1991) (parameters: degree-day factor for snowmelt /3snow = 3 mm we. d–1 °C–1; degree-day factor for ice melt Ace = 12 mm we. d–1 °C–1; saturation factor for the formation of superimposed ice Pmax = 0.6; standard deviation of short-term, statistical air-temperature fluctuations σstat = 5°C), sea-level changes derived from the spectral-mapping project (SPECMAP) record (Reference Imbrie, Berger, Imbrie, Hays, Kukla and SaltzmanImbrie and others, 1984) and the global-mean value qgeo = 55 mW m–2 for the geothermal heat flux.
4. Results
4.1. Ice-sheet extent and thickness
At the LGM around 20 kyr BP, the simulated ice sheet (Fig. 3) shows two different parts, namely, the Fennoscandian and the Barents–Kara ice sheets. The Fennoscandian ice sheet lies within its geologically reconstructed boundaries, but the eastern part is heavily over-glaciated in the Kara Sea and Taymyr areas.The total areal coverage of the Eurasian glaciation is 7.46106 km2, and the sea-level equivalent at that date is 39.5 m. If we compare this to the estimates of ice-equivalent sea levels for the LGM ice sheet by Reference Clark and MixClark and Mix (2002), we notice that the maximum Climate Long-Range Investigation Mapping and Prediction (CLIMAP) value for the Fennoscandian and Barents area ice sheets is 34.0 m. This is over 5 m lower than our result. However, looking at our resulting ice-sheet topography at the LGM, we notice that the over-glaciated eastern section makes a major contribution to the ice-sheet volume. The simulated maximum ice-sheet thickness in Fennoscandia at the LGM, approximately 3km, is reached over the Gulf of Bothnia.
Fig. 3. The simulated ice-sheet extent (thick solid line) and thickness (solid contours, spacing 1km) at the LGM (20 kyr BP). The thick dashed line indicates the LGM ice-sheet extent by QUEEN (Reference SvendsenSvendsen and others, 1999; Reference Thiede, Bauch, Hjort and MangerudThiede and others, 2001).
The PMIP results including UKMO21utilize the Ice-4G ice-sheet data (Reference PeltierPeltier, 1994), which have the best available spatial resolution. The PMIP runs are thus based on a heavily glaciated Eurasia which is not supported by the QUEEN reconstruction. It is likely that this leads to biased LGM temperatures which, in turn, are responsible for the over-glaciated Eurasian Arctic of our simulation.
An inverse approach has been used to determine the climate parameters (temperature, precipitation) in order to fit the modelled ice-sheet extent to the QUEEN observations. Preliminary results indicate that if the climate forcing defined by Equations (3) and (4) is modified such that winter precipitations are limited to 150 mm a– 1 and, simultaneously, summer temperatures are kept above 3°C during the glacial, the eastern Eurasian part will not be glaciated. This approach will be investigated in more detail in a future study.
4.2. Ice streams
By definition, an ice stream is a region of grounded ice which flows much faster than the surrounding ice on either side (Reference PatersonPaterson, 1994). The main references for palaeo-ice streams are Reference Denton, Hughes, Denton and HughesDenton and Hughes (1981) and Reference Stokes and ClarkStokes and Clark (2001). In these studies, the hypothesized palaeo-ice-stream locations are identified (see also Fig. 1). The main ice streams were marine, as nowadays in Antarctica and Greenland. However, terrestrial ice streams were active at glacial times as well (e.g. the inferred BIS (Reference Punkari and AberPunkari, 1993, Reference Punkari1997)).
The modelled fast-flow features in our simulations (Fig. 4b) can be seen as coarse approximations of real ice streams. However, it must be kept in mind that, first, the 40 km model grid is too coarse to resolve the small width of real ice streams, and second, as mentioned above, the model uses the shallow-ice approximation which does not include the particular ice-stream physics. In our model, fast flow arises from basal sliding over a temperate base and is therefore favoured by steep surface gradients and topographic depressions of the bedrock where the ice is thick.
Fig. 4. Simulation results for the Norwegian Channel and Baltic Basin ice-stream activity at 21,18,17,16,15 and 14.5 kyr BP. A and B1–B4 are simulated positions of fast-flow features. (a) Ice-sheet thickness (contour spacing 1km). (b) Flow velocity. For white areas 
for light shading 200 m a–1 ≤
dark shading 300 m a–1≤
for short arrows 400ma–1 ≤
for mid-length arrows 700 ma–1 ≤
for long arrows
(c) Basal temperature in °C. Dark (light) shading indicates a temperate base with temperate (cold) ice above. In all panels, the outermost contour shows the ice margin.
We obtain fast-flow regions at St Anna Trough, Franz Josef Trough, Bear Island Trough, the Trondheim area and the Norwegian Channel. These are located in topographic troughs, where the base is temperate (Fig. 4c), and they coincide with locations of marine-based ice streams indicated by geologic evidence (see section 1).
Terrestrial ice streams were inferred for the Baltic area in several stages of glaciation at varying locations (see section 1). In our simulation, ice flow is generally slow in the Baltic area during the advance phase of the ice sheet. However, during the retreat phase after the LGM a temperate ice base forms in that region, and ice flow is significantly accelerated.
The threshold seabed elevation on the continental shelf which is allowed to glaciate (zthresh) is a sensitive determinant of modelled maritime ice-stream locations. Therefore, it is important to choose an appropriate value for this parameter, which was achieved by matching simulated and reconstructed marine ice-sheet margins at the LGM (Fig. 3). Otherwise, maritime ice-stream locations are quite robust to changed input parameters. By contrast, terrestrial ice-stream locations differ significantly with different parameter sets.
5. Discussion
5.1. Norwegian Channel area
The Norwegian Channel ice stream (NCIS) is regarded as important for the Fennoscandian ice-sheet dynamics (Reference Stokes and ClarkStokes and Clark, 2001). The activity of the NCIS has been discussed by Sejrup and others (1998, 2000), who conclude that the NCIS was active during the periods 50–40, 28–22 and 18– 15 kyr BP. Between these periods, the stability and duration of the stream activity is not known. These results are based on lithology, biostratigraphy and geochronological investigations.
The results of our simulation for six time-slices are shown in Figure 4. Evidently, fast flow prevails in the NCIS area (A; Fig. 4). It occurs first at 28 kyr BP, with high surface velocities of 1400 m a–1 in the Norwegian Channel, the velocity near the base being approximately 1000 ma–1. The NCIS is situated near the margin of the ice sheet, with a large surface gradient, and there is a temperate base with temperate ice above. The period 28–18 kyr BP is rather stable in our simulation, with fast ice flow from the northeast to the southwest. The area with a temperate base overlain by a temperate-ice layer is rather constant (around 20 000 km2), with minor areal changes that do not affect the fast flow.
At 18 kyr BP, the simulated ice sheet starts thinning in the western section. There is a large area around southern Sweden with temperate-base conditions but cold ice above. There the velocities are approximately 400–700ma–1.
At 17 kyr BP the simulated NCIS fast-flow area has widened, but the peak velocity has diminished to 900 ma–1. At this time, there are some shifts on the main fast-flow direction to the south. After this, the temperate-base area diminishes considerably. The last time when there is fast flow with an ice sheet terminating in water and velocities of about 500 m a–1 is at 15 kyr BP. Thereafter, the flow is slowly fading, with velocities of 300 ma–1.
In our simulation, the fast flow in the NCIS area is due to temperate basal conditions and resulting basal sliding. From 28 to 17 kyr BP there is a temperate base with temperate ice above in the Norwegian Channel. As the ice sheet thins, the temperate ice above the base first becomes cold as cold ice flows in from the mountains. As the ice sheet thins further, the temperate base disappears and fast flow ceases.
5.2. Baltic region
The simulated BIS (Fig.4) has an extensive spatial and temporal dimension. The fast-flow regions coarsely agree with the ice streams described by Reference Boulton, Dongelmans, Punkari and BroadgateBoulton and others (2001). In our results, lobate structures form when ice is flowing to topographic depressions.
The main simulated advance towards the southern-Sweden/Denmark region starts around 24 kyr BP. At this time, there is a temperate base with cold ice above at the ice-sheet margin in southern Sweden, and the ice sheet ends in the sea. Thevelocities are approximately 200 m a–1. There is stable thickening of the ice sheet until 19 kyr BP, with alternating temperate and cold ice above a temperate base near the margin. At 19 kyr BP the ice in the marginal areas is warm and mobile, with typical velocities around 600 ma–1. From then on, the ice sheet quickly retreats. Nevertheless, and in contrast to the NCIS, the simulated BIS is still active after 15 kyr BP.Velocities of up to 700m a–1 are maintained, and occasionally peak values of 41000 m a–1 are reached.
The simulated fast-flow features connected to separate sub-lobes develop in the following way:
The Danish sub-lobe (B1) advances steadily at 24 kyr BP from areas west of the Baltic Basin with low velocities. At 18 kyr BP it starts to retreat. The ice-flow velocity at this time is approximately 600 m a–1, and there is a wide marginal area (B1–B3) with temperate base and temperate ice above.
The Oder sub-lobe (B2) advances at 23 kyr BP. The flow starts to accelerate at 18 kyr BP, starting from approximately 600ma–1, and attains a peak of 1200ma–1 at
17 kyr BP.This is due to the temperate marginal conditions.
The Wisla sub-lobe (B3) advances at 23 kyr BP. Temperate basin appears in the area at 18 kyr BP, increasing the ice flow. The velocity is 900 m a–1 at 16 kyr BP, 1000 m a–1 at 15 kyr BP and 1400 ma–1 at 14.5 kyr BP.
The Kurshkiy sub-lobe (B4) advances to a maximum position at 21 kyr BP, and has the highest-velocity phase at 14.5 kyr BP, with values of 1400 m a–1.
During the ice-margin retreat, the former advance-phase lobate structure is coarsely preserved. When the ice sheet decreases at 16 kyr BP and the southwestern part ends in the sea, a large ice front is mobilized. A periglacial lake is formed in the Baltic basinat13 kyr BP.This givesrise to a few fast-flow areas (1800 m a–1) at the ice margin.
6. Conclusion
The following conclusions can be drawn from our modelling study:
The extent of the Fennoscandian LGM ice sheet is reasonably reconstructed by our simulation.
In the Kara Sea and Taymyr areas, the glaciation obtained is too extensive. Evidently, in these regions the applied climate forcing is not valid. Therefore, in future work we will seek to use general circulation model (GCM) results for the glacial climate based on the QUEEN instead of the Ice-4G ice-sheet extent. Further, in order to minimize the influence of systematic GCM errors, the glacial climate input will be based on anomalies from GCM results (LGM minus present) instead of absolute GCM output.
Despite the model limitations of the 40 km grid and the simplified ice-sheet dynamics, observed main ice-stream locations can be identified as fast-flow features in the simulation results, and they compare well with the geological evidence.
The modelled fast-flow features arise from basal sliding over a temperate base. This is favoured by, but not necessarily dependent on, topographic depressions of the bedrock.
The lobate structure in our results is a result of the advance phase. The lobate structure is preserved in the following decay phase.
The model domain should be increased for future simulations in order to include the British Isles and to extend further to the east.
Acknowledgements
The authors express their gratitude to M. Saarnisto (Geological Survey of Finland) and J.-P. Lunkka (University of Oulu) for valuable suggestions. The reviews by S. J. Marshall and an anonymous reviewer, as well as the very detailed comments by the scientific editor C. S. Hvidberg, helped to improve the paper substantially. The resources of the Centre for Scientific Computing in Finland made this study possible. Funds from the Academy of Finland (project 47047) are also acknowledged.
