Introduction

Ice-sheet models used to solve an array of problems relating to past, present and future ice sheets are continuously growing more complex. As predicted by Waddington (1987), present-day thermomechanical ice-sheet models are no longer primarily limited by computer capacity, but rather by poorly known data required for boundary conditions. At the ice-sheet–bed interface, the most important, and also least known, boundary condition is the geothermal heat flow (or heat-flow density (HFD)). The heat flux at the bed of ice sheets affects their stability because it determines whether the base of the ice sheet will be cold- or wet-based and, if wet-based, also how much of the water generated is available for accelerated glacier sliding.

Ice-sheet modeling studies have shown that basal ice temperatures are also sensitive to relatively small changes in geothermal heat flow (e.g. Greve and Hutter, 1995; Siegert and Dowdeswell, 1996; Näslund and others, 2000; Tarasov and Peltier, 2003). There is strong reason to believe that using realistic non-uniform geothermal heat-flow values in ice-sheet modeling would considerably improve model calculations of basal ice temperatures and basal melt rates, in line with results of sensitivity tests from the Greenland ice sheet (Greve and Hutter, 1995). However, no ice-sheet modeling study up to date has used a detailed, spatially varying geothermal heat-flow dataset as basal thermal boundary condition.

Geophysical and borehole measurements have shown that there may be substantial spatial variation in geothermal heat flow at the surface of the lithosphere, on both regional and local scales. For example, this is the case for the Precambrian shields underlying the former Scandinavian and Laurentide ice sheets (e.g. Lindén and others, 1983; Pollack and others, 1993; Artemieva and Mooney, 2001). Geothermal heat flow under the present Antarctic and Greenland ice sheets probably also varies considerably. In the vicinity of the North Greenland Icecore Project (North-GRIP) drilling site, for example, Dahl-Jensen and others (2003) report spatial variations of geothermal heat flow from ∼50 to 200 mWm^–2.

In Scandinavia, the geothermal heat flow has been known to vary by a factor of at least two (Eliasson and others, 1992) which, in the context of thermodynamic ice-sheet modeling, is a considerable variation. However, the traditional way of specifying the basal thermal boundary condition in ice-sheet modeling has been to apply one uniform value of geothermal heat flow to the entire model domain. The magnitude of this geothermal heat-flow value has often been set according to the general geological settings for the area to be studied. Thermodynamic modeling of Fennoscandian ice sheets has for instance used a single geothermal heat-flow value typical for the Precambrian shield, the most common reported value being 42 mWm^–2 (Boulton and Payne, 1993; Huybrechts and T’siobbel, 1995; Payne and Baldwin, 1999). Others have used higher values, around 50–55 mWm^–2 (e.g. Hindmarsh and others, 1989; Forsström and others, 2003). For a 3 km thick ice sheet in steady state, a 20% error in geothermal flow generates a 6 K error in basal temperature (Waddington, 1987). In several published ice-sheet modeling studies of Fennoscandia, the value used on geothermal flow is not mentioned at all. Recently, some ice-sheet modeling studies have been made for other parts of the world where spatially varying geothermal heat-flow datasets have been used (e.g. Tarasov and Peltier, 2003). However, these are still geothermal datasets with very low resolution, based on a restricted number of point measurements. They do not depict possibly important local and regional variations in heat flow.

For large parts of the regions covered by Fennoscandian ice sheets, data are available for specifying a detailed, spatially variable geothermal heat-flow field. We have done this for the core area (Sweden and Finland) of the Fennoscandian Last Glacial Maximum (LGM) ice sheet, and have applied these data as basal thermal boundary condition for a numerical ice-sheet model to investigate the effects of the variable high-resolution heat flow on ice-sheet dynamics. The overall aims of the study were to: (1) establish what uniform geothermal value should be used for traditional ice-sheet modeling over Fennoscandia, and (2) present a high-resolution non-uniform (i.e. spatially variable) dataset describing geothermal heat-flow variability for future modeling.

Geothermal Heat Flow

The geothermal conditions of the Baltic shield and the calculation of geothermal heat flow are developed in a separate paper by Näslund and others (unpublished information). Therefore, only a summary of the calculation methods is provided here.

The geothermal heat flow, or surface HFD, in cratonic areas consists of two components: (1) heat produced within the mantle and core of the Earth and (2) heat produced within the crust itself. The contribution from the Earth’s interior (so-called Moho or reduced heat flow) arises from the cooling of the Earth and formation of a solid core, and from radiogenic heat production (e.g. Turcotte, 1980; Pollack and others, 1993). The crustal component consists of a radiogenic heat production where heat is produced by the natural radioactive decay of primarily ²³⁸U, ²³²Th and ⁴⁰K (e.g. Furlong and Chapman, 1987; Nielsen, 1987). The Moho heat flow has a smooth spatial variation, possibly depending on mantle convection cell distribution (e.g. Beardsmore and Cull, 2001), whereas the spatial variation in radioactively decaying nuclides in the lithosphere generates a heat flow with large spatial variations. The surface HFD can be estimated by a heat-flow–heat-production (Q–A) relationship of the form

(Birch and others, 1968; Lachenbruch, 1968; Beardsmore and Cull, 2001), where Q is the surface heat-flow density (or geothermal heat flow), q ₀ is the Moho heat flow, D represents the vertical distribution of heat-producing elements in the lithosphere and A ₀ is the radiogenic heat production from the surface near rocks.

Cratons constitute the oldest part of continental crust and usually exhibit low heat-flow values (e.g. Čermak and others, 1993). The east European craton is an area with especially low HFD values and is surrounded by normal to high values in western, central and southeastern Europe. The lowest heat-flow values are found in the Ukrainian and Baltic shields (e.g. Čermak and others, 1993). Artemieva and Mooney (2001) present a compilation of 300 observed HFD values for the Baltic shield, sufficient to show the general heat-flow pattern but not regional or local variations. Their compilation shows increasing HFD values from the northeast (35–50mWm^–2) to the south and southwest (60–70mWm^–2). The regional HFD pattern does not correlate with gravity variations (Balling, 1984), magnetic anomalies (e.g. Riddihough, 1972) or crustal thickness (Meissner and others, 1987; Čermak and others, 1993). However, within the Baltic shield, as in other areas with similar geological settings, there is close correlation between HFD and regional geological units, with higher heat flow from acid (commonly granitic) areas and lower heat-flow values from basic areas (e.g. Landström and others, 1980; Malmqvist and others, 1983).

The linear Q–A relationship in Equation (1) has been regarded as an oversimplification for some purposes (e.g. Furlong and Chapman, 1987; Nielsen, 1987). However, several studies have presented information on heat flow and heat production for the Baltic shield with most of the regression lines falling into relatively homogeneous groups, which led Kukkonen (1993) to suggest that a linear Q–A relationship may be valid for the Baltic shield. In any case, the linear Q–A relationship is sufficient for the purposes of the present study.

Our calculation of HFD values for Sweden and Finland is based on detailed datasets from numerous g-emission measurements from bedrock and till. Airborne surveys of γ emission have been carried out by the Geological Survey of Sweden (SGU), sampling data at 70 m intervals along flight-lines with 17km separation. In order to avoid shielding effects from vegetation and lakes, only data from exposed bedrock and till surfaces were used in our calculations. The Finnish dataset is based on radiometric g-emission measurements of 1054 till samples providing a full spatial coverage of the country (Kukkonen, 1989).

The calculation of HFD is performed in several steps. First, the concentration of ²³⁸U, ²³²Th and ⁴⁰K is calculated from the g-emission measurements, using information from detailed reference measurements over calibration plates with well-known isotopic concentrations. Near-surface heat production was then calculated from the concentration of radioactive nuclides by

where A ₀ is the radiogenic heat production, ρ is the bedrock density, C _U and C _Th are the ²³⁸U and ²³²Th concentrations in ppm, and C _K is the ⁴⁰K concentration in %. An average rock density of 2.65 kgm^–3 was used (Lindén and others, 1983).

The A ₀ values calculated by this method only describe the heat production near the surface and not the geothermal conditions at depth. The most important property to consider here is the Moho heat-flow contribution, which must be added in accordance with Equation (1). We added the distributed Moho heat flow by Artemieva and Mooney (2001). Furthermore, in Equation (1), the thickness of the upper part of the crust (D) was set to 13 km. This is the mean value of previously used D values for Finland and Sweden (11 and 15 km respectively) (Lindén and others, 1983; Kukkonen, 1989). The results are shown in Figure 1.

Fig. 1. Geothermal heat flow (or HFD) distribution for Sweden and Finland based on (a) new calculated values partly based on measured g-emission and a heat-flow–heat-production relationship (Q _mean = 44.9 mWm^–2), and (b) previously published observed values obtained from borehole data (H.N. Pollack and others, data downloaded from US National Geophysical Data Center, ftp://ftp.ngdc.noaa.gov/Solid_Earth/Global_Heatflow/) (Q _mean = 45.9 mWm^–2). The calculated high-resolution dataset (a) has a larger range, with minimum and maximum values of 30 and 83 mWm^–2, while corresponding values for the observed low-resolution dataset (b) are 36 and 60 mWm^–2.

The difference in g-emission sampling method may have an impact on calculated HFD values, with airborne measurements generally underestimating true ground-surface g-emission values (e.g. Lindén and others, 1983). A correction of the Swedish data against Finnish data (based on laboratory measurements of field samples) was therefore applied by matching the datasets along the Swedish–Finnish border. Finally the HFD dataset was resampled to a grid with 5 km resolution using cubic convolution.

To provide HFD coverage for the entire model domain, we have added surrounding data from the much coarser observed global HFD dataset by H.N. Pollack and others (data downloaded from US National Geophysical Data Center, ftp://ftp.ngdc.noaa.gov/Solid_Earth/Global_Heatflow/) to our calculated high-resolution dataset for Sweden and Finland (Figs 1 and 2). Most of Pollack and others’ data originate from drillhole observations.

Corrections of Geothermal Heat

It should be noted that the HFD dataset over Scandinavia does not yet include the Barents Sea region. Furthermore, in the new HFD dataset the Baltic Sea between Sweden and Finland is not covered by actual measurements or observations but relies on the Pollack interpolation between surrounding land areas. Furthermore, large areas of the Baltic Sea bed comprise sedimentary rocks overlying basement. The sedimentary units do not contribute significant radiogenic heat production themselves (Slichter, 1941). The effect of the overlying sediment strata on geothermal heat flow has not yet been accounted for in the dataset. Other factors may also affect the amount and distribution of HFD, such as changes in climate (including the effect of the presence of ice sheets), groundwater flow, sedimentation rates in basins (e.g. the Baltic depression), basement relief and topography, and erosion on the lithospheric surface. In a coming version of the dataset, some of these corrections will be made. Until then, the presented HFD dataset should be considered as preliminary. However, the quality of the presented dataset is high enough for evaluation of the effects of such a dataset in ice-sheet modeling, and also to extract first-order statistics on HFD over Fennoscandia.

The Ice-Sheet Model

The ice-sheet model used to simulate the Weichselian glaciation is the University of Maine Ice Sheet Model, which is a time-dependent, non-steady-state thermomechanical ice-sheet model (Fastook and Chapman, 1989; Fastook, 1994; Fastook and Holmlund, 1994; Fastook and Prentice, 1994; Johnson and Fastook, 2002; Näslund and others, 2003). The model uses the finite-element method to solve the mass-continuity equation with a vertically integrated momentum equation. Isostatic changes of the bed and eustatic sea-level changes are taken into account by the model. Important input to the ice-sheet model is a digital elevation model (DEM) of landscape topography and a parameterized mass balance at each grid node (Fastook and Holmlund, 1994). A finite-element temperature solver, with 20 layers through the ice thickness, is coupled to the mass-continuity solver to generate internal ice temperatures. Derived ice temperatures, together with density variations with depth, control ice hardness and ice flow. The thermodynamic calculation accounts for vertical diffusion, vertical advection, and heating caused by internal shear. Heating by internal shear is important for ice deformation near the bed where most of the shearing occurs. Horizontal diffusion and advection of heat within the ice sheet is not considered. Especially in cases of fast ice flow, horizontal advection is of importance and should probably not be neglected as it is here. Another thing to consider in our study is that the model does not incorporate heat conduction in the upper layer of bedrock. These weaknesses in the present approach will be remedied in future work.

The numerical model has been augmented with a nesting capability in which high-resolution runs of small sub-domains can be modeled using the coarser full ice-sheet run as boundary condition for the selected sub-domain.

Experimental Procedure

The ice-sheet model was first run to simulate the entire Weichselian ice sheet through the last glacial cycle, using the Greenland Icecore Project (GRIP) temperature record as input climate signal and 50 × 50 km resolution topography extracted from the ETOPO2 DEM. Ice configurations were calibrated against dated LGM and Younger Dryas ice-marginal positions known from geological observations. Model simulations were made using two different HFD scenarios employing (1) the commonly used estimated HFD of 42 mWm^–2 and (2) the full new distributed HFD dataset. This experiment provides information on the resulting differences in basal temperature distributions and melt rates between these two cases. Furthermore, we employ nested high-resolution modeling to study how the distributed high-resolution HFD data affect ice-sheet basal conditions on a local scale. The area chosen for this study is a ∼115×100km² area located in southern Sweden, which includes one of two suggested sites for a future Swedish nuclear-waste depository. Full-resolution data from the ETOPO2 DEM (∼2×3.5km grid resolution) were used as topographical input to this model simulation.

Results

The calculated distributed geothermal dataset for Sweden and Finland is shown in Figure 1a, while Figure 1b shows interpolated observed HFD data from boreholes from the global dataset of H.N. Pollack and others. The general trends in the two datasets are the same. In addition, the mean values of the datasets are also in effect the same (45 and 46 mWm^–2). This shows that the methodology used in calculating the new distributed data for Sweden and Finland is appropriate, and also verifies the quality of the new data. One important difference seen in comparing the datasets in Figure 1 is that the new calculated high-resolution dataset has a larger range, the maximum and minimum values being 30 and 83 mWm^–2 respectively. The corresponding values for the observed low-resolution dataset are 36 and 60 mWm^–2. This is a result of the high-resolution dataset including distinct regional to local HFD anomalies not present in the coarser global dataset.

The new complete HFD dataset for Scandinavia, including the low-resolution global data surrounding Sweden and Finland, is presented in Figure 2. As expected, large local variations in heat flow occur, as readily seen in the data for Sweden and Finland. Within the area covered by the LGM ice-sheet margin, the HFD varies between 30.2 and 83.2 mW m^–2, with a mean of 48.7 mWm^–2. This shows that the HFD within the LGM margin varies by a factor of 2.8, which is considerably more than the previously known, but typically neglected, variation of 2. For the entire ice sheet, the total basal meltwater production for the whole glacial cycle is 1.4 times higher when using the full HFD dataset shown in Figure 2 as geothermal boundary condition, than using a traditional uniform value of 42 mWm^–2 (265Mkm³ of water instead of 190Mkm³).

Fig. 2. Heat-flow density distribution for the area covered by the LGM ice sheet. The data were used as basal thermal boundary condition for ice-sheet modeling focusing on basal temperature and hydrology. High-resolution data were calculated for Sweden and Finland (Fig. 1b), while the much coarser data for surrounding regions are based on observed HFD (H.N. Pollack and others, ftp://ftp.ngdc.noaa.gov/Solid_Earth/Global_Heatflow/). Within the LGM ice-margin position the geothermal heat flow ranges between 30 and 83 mWm^–2, with an average value of 49 mWm^–2. This is the most appropriate value to use in modeling the Fennoscandian ice sheet if adopting the traditional simplified treatment of geothermal heat using one uniform value. The square denotes the area selected for high-resolution modeling (Fig. 3).

For the high-resolution study area in southern Sweden, the uniform and distributed HFD cases are shown in Figure 3b and c. In the distributed case (Fig. 3c), HFD values are low in the northwest of the area (<45mWm^–2), and generally higher towards the southeast, closer to the coast. An exception, with low HFD values also at the coast, occurs in the south-central part of the study area. A strong anomaly of exceptionally high HFD is seen in the central part of the study area (Fig. 3c).

Fig. 3. (a) Topography and bathymetry of a ∼115×100km² subdomain in southern Sweden. (b) Uniform HFD basal boundary condition (average of values in (c)). (c) New spatially variable high-resolution HFD basal boundary condition. (d, e) The period of cold-based conditions expressed as fraction of total ice-covered time using the HFD basal boundary conditions in (b) and (c), respectively.

Examples of results of the high-resolution model runs for a uniform and a distributed HFD case are shown in Figures 3d and e and 4.

Discussion

The ice-sheet model was calibrated against geological data so as to fit in time and space with both the LGM position and the Younger Dryas zone. Using the new high-resolution dataset on geothermal heat flow (with a mean value of 49 mWm^–2 instead of a uniform value of 42 mWm^–2) in modeling the Fennoscandian ice sheets has several effects. First, the area of temperate bed conditions increases with the new dataset. This difference is also true through time. Second, the total amount of water produced by basal melting increases significantly. Our model runs yielded a 1.4 times higher total water production from the ice sheet for the full high-resolution HFD case than for the uniform 42 mWm^–2 case, showing the importance of using a proper value on HFD in modeling simulations aimed at studying basal hydrological conditions.

Figure 3d and e show the period of cold-based ice cover expressed as a percentage of total ice-covered time for the high-resolution modeling. In both panels, high values in the northwest of the study area indicate prolonged cold-based conditions during the glacial cycle, while low values in the southeast indicate prolonged wet-based ice coverage.

In Figure 3d, the more cold-based conditions in the northwest are to a large degree caused by the higher elevations of the bed in this area (Fig. 3a), while the warmer bed conditions in the southeast are caused by low bed elevations. One effect of including the new realistic HFD dataset is shown in Figure 3e. The temperature gradient at the ice-sheet bed has increased significantly over the area, in the northwest resulting in longer periods of cold bed conditions and larger areas covered by cold-based ice. This is because both the high elevation and the contribution by the low HFD in this region (Fig. 3c) act in the same direction, to produce lower basal ice temperatures. At the same time, the areas in the southeast experience the opposite situation: larger areas have become warm-based during the glacial cycle. Especially large changes in basal temperature are observed over the area of the warm HFD anomaly in the central part of the area (Fig. 3c and e). One should note that our simulation does not include temporal variations in the HFD induced by the advancing and retreating ice sheet. Such effects will be investigated in forthcoming studies.

In this specific study area, the use of realistic HFD boundary conditions resulted in a clearly accentuated basal temperature pattern. Furthermore, looking at the time-transient behavior of the ice, areas with high HFD also experienced more rapid heating of cold beds during some phases of the glaciation, while the opposite was the case for areas with low HFD values.

Basal melt rates were also clearly affected by the use of the new HFD dataset in the high-resolution model run. Areas of high HFD yielded regionally higher local meltwater rates (two or more times higher than when using the uniform HFD value), while areas of low HFD sometimes reduced basal melt rates by a factor of two. As shown in Figure 1, the spatial variations in HFD are greater for other regions in Scandinavia, such as northern Sweden, than in the high-resolution study area examined here. Similar ice-sheet modeling studies in such areas would most likely yield results with larger variations in basal conditions and hydrology when using realistic high-resolution HFD data.

The results of the study thus clearly show that numerical ice-sheet modeling with a realistic spatially variable heat-flow dataset is essential for proper identification and description of thermal and hydraulic conditions beneath ice sheets. Furthermore, in cases where the amount of basal sliding in the model is coupled to the amount of water at the bed (e.g. Johnson and Fastook, 2002), the use of realistic, distributed HFD data becomes increasingly important for ice flow, for instance if ice streams are to be studied.

One example of the modeled basal melt-rate distribution for the area in southern Sweden is shown in Figure 4a, with associated surface velocities presented in Figure 4b. The high-resolution model run, using the distributed HFD boundary condition, provides detailed information on local flow events. One example of such an event is an ice-stream-like feature in the trough between Öland and the Swedish mainland during the last deglaciation (Fig. 4b). This feature is associated with the previously known generally higher ice-flow values of the Baltic depression. However, this is a short-lived (∼14 500–14 100 BP) but very intense feature, with surface velocities exceeding 500ma^–1. A comparison between this model result and a geomorphologic interpretation of landform assemblies shows that the fast-flow feature from the model corresponds both geographically and directionally to a specific flow fan identified by Kleman and others (1997, fig. 7). This ice-stream-like feature was clearly more active in the model run using realistic HFD than in the uniform HFD case, largely because the presence of the high-HFD anomaly in the central part of the study area produces more basal meltwater. The new distributed dataset thus strongly enhances a flow feature that can be observed in the geomorphologic record, suggesting both a cause and timing for the event producing the landforms. This strongly suggests that ice-sheet modeling aimed at studying the presence and distribution of paleo ice streams (e.g. in former Weichselian and Laurentide ice sheets) should employ realistic high-resolution HFD data. To this end, existing observed global HFD datasets are too coarse to capture highly significant regional variations in HFD (Fig. 1). For Scandinavia, and possibly also for North America, the HFD datasets must have a resolution of 10–15km or higher to capture such geothermal heat anomalies.

Fig. 4. (a) Basal melt rate for the sub-domain in Figure 3a during a 300–400 year long local ‘ice-streaming’ event at ∼14 500–14 100 _BP, during the last deglaciation; (b) associated ice surface velocities. The model run uses the spatially variable HFD basal boundary condition in Figure 3c.

Conclusions

The construction of a spatially distributed HFD dataset for Scandinavia shows large local variations in geothermal heat flow. For the area covered by the LGM Weichselian ice sheet, HFD values range between 30 and 83 mWm^–2, i.e. the geothermal heat flow within the LGM margin varies by a factor of as much as 2.8. The average HFD value within the LGM margin is 49 mWm^–2, which is 17% higher than the typical value used in several previous ice-sheet modeling studies. The introduction of a realistic high-resolution dataset on geothermal heat flow significantly affects the absolute values and distribution of basal temperature and basal melt rates beneath the modeled ice sheet. For example, using the new distributed dataset on geothermal heat flow, instead of a traditional uniform value of 42 mWm^–2, yields a 1.4 times larger total basal meltwater production from the Fennoscandian ice sheet for the last glacial cycle. In sheet-modeling studies where basal ice temperatures are important (e.g. studies focusing on basal meltwater production or erosional capacity), a correct and spatially varying dataset on geothermal heat flow should be used. Most likely, regional to local variations in geothermal heat flow also need to be considered for proper identification and treatment of thermal and hydraulic bed conditions when studying Laurentide, Greenland and Antarctic ice sheets.