Comparing a thermo-mechanical Weichselian Ice Sheet reconstruction to reconstructions based on the sea level equation : aspects of ice configurations and glacial isostatic adjustment

In this study we compare a recent reconstruction of the Weichselian Ice Sheet as simulated by the University of Maine ice sheet model (UMISM) to two reconstructions commonly used in glacial isostatic adjustment (GIA) modelling: ICE-5G and ANU (Australian National University, also known as RSES). The UMISM reconstruction is carried out on a regional scale based on thermo-mechanical modelling, whereas ANU and ICE-5G are global models based on the sea level equation. The three models of the Weichselian Ice Sheet are compared directly in terms of ice volume, extent and thickness, as well as in terms of predicted glacial isostatic adjustment in Fennoscandia. The three reconstructions display significant differences. Whereas UMISM and ANU includes phases of pronounced advance and retreat prior to the last glacial maximum (LGM), the thickness and areal extent of the ICE-5G ice sheet is more or less constant up until the LGM. During the postLGM deglaciation phase ANU and ICE-5G melt relatively uniformly over the entire ice sheet in contrast to UMISM, which melts preferentially from the edges, thus reflecting the fundamental difference in the reconstruction scheme. We find that all three reconstructions fit the present-day uplift rates over Fennoscandia equally well, albeit with different optimal earth model parameters. Given identical earth models, ICE-5G predicts the fastest present-day uplift rates, and ANU the slowest. Moreover, only for ANU can a unique best-fit model be determined. For UMISM and ICE-5G there is a range of earth models that can reproduce the present-day uplift rates equally well. This is understood from the higher present-day uplift rates predicted by ICE-5G and UMISM, which result in bifurcations in the best-fit upperand lowermantle viscosities. We study the areal distributions of present-day residual surface velocities in Fennoscandia and show that all three reconstructions generally over-predict velocities in southwestern Fennoscandia and that there are large differences in the fit to the observational data in Finland and northernmost Sweden and Norway. These difference may provide input to further enhancements of the ice sheet reconstructions.


Introduction
Fennoscandia has been a key area in the development of theories and models of glacial isostatic adjustment (GIA) due to the unique temporal and spatial coverage of observational data (e.g.Ekman, 2009), and the region remains an important study area; see summaries of recent work in e.g.Plag et al. (1998), Whitehouse (2009) and Steffen and Wu (2011).To model the GIA process two components are needed: an earth model and an ice sheet reconstruction, collectively referred to as a GIA model.Today several reconstructions of the Weichselian Ice Sheet are available, both regional models and as part of global models (e.g.Peltier, 2004;Lambeck et al., 2010;Näslund, 2010).The ice sheet properties that are most important for GIA studies are the areal extent and the distribution of ice thickness.These two are however very different in terms of how difficult they are to constrain.
Published by Copernicus Publications on behalf of the European Geosciences Union.

P. Schmidt et al.: Weichselian Ice Sheet reconstructions
Whereas the areal extent often can be reasonably well determined from geological markers such as moraines, usually very limited data are available on the ice thickness and generally only from mountainous regions.The thickness therefore often has to be determined by indirect methods.Important to notice is that most ice models require an earth model for the reconstruction of ice thickness; a coupling therefore often exists between the ice sheet and the optimal earth model parameters in a GIA model.
Ice sheet reconstructions can broadly be categorised into two groups.The first, classical approach base the reconstruction primarily on geological markers of the extent of the ice sheet at different times (e.g.Denton and Hughes, 1981;Lambeck, 1993;Lambeck et al., 1998;Tushingham and Peltier, 1991;Peltier, 1994).The ice thickness is then adjusted such that the solution to the sea level equation fits available data (mainly relative sea level (RSL) and tide-gauge data, but more recently also GPS data).For global models, estimates of the global mean sea level are used to constrain the total volume of ice at different times.Early models of this type did not extend further back in time than to the last maximum extent of the ice sheet, due to sparsity of older RSL data and geological markers (e.g.Peltier and Andrews, 1976).
The second group reconstructs ice sheets from physical (thermodynamical) principles, often using palaeo-climate data to govern the evolution of the ice sheet (e.g.Forsström and Greve, 2004;Zweck and Huybrechts, 2005;van den Berg et al., 2008).This provides an ice sheet which, within the limitations of the model implementation, behaves as a real ice sheet in terms of basal sliding, ice streams, ice thickness distribution and growth and decay properties.Geological markers and RSL data may be used to constrain the reconstruction.As the interdependence between ice model and earth model varies in the two types of reconstructions, they may provide complementing information on the properties of the earth.
In this study we compare a recent thermo-mechanical reconstruction of the Weichselian Ice Sheet, the University of Maine ice sheet model (UMISM) (Näslund, 2010), with two models constrained by geological markers and RSL observations: ANU (Lambeck et al., 2010, Australian National University, also known as RSES) and ICE-5G (Peltier, 2004).This is not the first study where different ice sheets have been used and compared.Schotman and Vermeersen (2005) used an earlier version of ANU (Lambeck et al., 1998), a modified version of ICE-3G (Tushingham and Peltier, 1991) and a thermo-mechanical reconstruction to investigate the effect on geoid heights to the ice load history in models with shallow low-viscosity earth layers.Wu et al. (2010) used FBK8 (Lambeck et al., 1998) and ICE-4G (Peltier, 1994) to investigate the sensitivity of the present-day rebound to the ice sheet thickness.Steffen et al. (2010) also used the ANU version by Lambeck et al. (1998, referred to as RSES) as well as ICE-3G, ICE-4G and ICE-5G and compared to satellitebased gravity data (GRACE).That study focused on the use of GRACE data in GIA modelling but concluded that both RSES and ICE-5G are adequate for studying the GIA process in Fennoscandia.Zhao et al. (2012) used ANU and ICE-5G and inverted for the optimal earth model parameters using the Lidberg et al. (2010) processing of the Bifrost GPS data, showing that both reconstructions can equally well fit observation of the present-day rebound in Fennoscandia.van der Wal et al. (2013) used ICE-5G as well as a reconstruction described therein, and compared to the Lidberg et al. (2007) processing of the Bifrost GPS data as well as RSL data.The study focused on the effect of the mantle rheology on GIA predictions and concluded that a wet rheology consistent with laboratory data could be found for which both ICE-5G and their own reconstruction could explain RSL data, while uplift rates predicted by their ice sheet that are close to observed present-day uplift demand a dry rheology.
Common to previous studies where multiple ice sheet have been used is a focus on the observational data and/or the earth model parameters.Further, most of the studies only presented the results of a limited number of earth models, thereby limiting the understanding of the differences between the reconstructions.Here we turn the primary attention to the ice sheet reconstructions and their GIA predictions, with the intention to understand how differences in the reconstructions lead to differences in the GIA predictions.The second objective of this study is to investigate how UMISM compares to ANU and ICE-5G, given that the former reconstruction has not been optimised to observations of the GIA process, while the latter two have been.Our comparison stretches back to 69 kyr before present (BP) up until today, which is significantly longer than any of the previous studies.
This paper is organised in three parts, each covered by two sections.In the first part we will focus on the ice sheet reconstructions themselves, with a description of each reconstruction in Sect. 2 followed by an inter-comparison in Sect.3. In the second part we focus on GIA by first describing our GIA model implementation in Sect. 4 and then the GIA predictions of each reconstruction in Sect. 5.In the third and last part of this work we summarise and discuss the results of the first two parts in Sect.6 and finish off with our conclusions in Sect.7.

Ice sheet reconstructions
Observations of the GIA process, as well as geological markers, are usually dated using the C-14 method.Most of the early reconstructions are therefore given in C-14 years rather than calender years.As these two timescales differ by up to 3.5 kyr around the time of the last glacial maximum (Bard et al., 1990), it is important to note that all three reconstructions considered here, as well as all observational data used herein, are dated in calender years.An event commonly referred to when discussing ice sheet reconstructions is the last glacial maximum, LGM.This often refers to the last maximum advance of the ice sheet.However, the conditions governing the advance and retreat of an ice sheets will vary from place to place.Therefore the maximum advance will not be a synchronous event in all parts of a large ice sheet.In the case of the Weichselian Ice Sheet the time span enclosing the last maximum advance may be as long as 10 kyr (e.g.Boulton et al., 2001;Forsström and Greve, 2004).
In this study we have chosen to define the LGM as the point in time of the last occurrence of the maximum volume of the ice sheet, remembering that this may not coincide with the maximum areal extent, thickness or advance of the ice front.

UMISM
The UMISM ice sheet reconstruction (Näslund, 2010) uses the October 2004 version of the thermo-mechanical University of Maine ice sheet model (Fastook and Chapman, 1989;Fastook, 1994;Fastook and Prentice, 1994).UMISM was part of the "European ice sheet modelling initiative model inter-comparison experiment" and yielded output in agreement with other thermodynamic ice sheet models (Huybrechts et al., 1996;Payne et al., 2000).The reconstruction has previously been used in GIA modelling for assessment of shoreline migration (Whitehouse, 2006) as well as fault stability (Lund et al., 2009).
In the present simulation UMISM was used for a regionalscale reconstruction of the Weichselian Ice Sheet on an equidistant-resolution (50 km × 50 km) grid every 100 yr since 120 kyr BP.The ice sheet constitutes three main subsystems: mass balance, ice movement and ice temperature for which the model solves the conservation of mass, momentum and energy equation respectively.The UMISM model used for the present reconstruction uses the shallow-ice approximation for solving stresses and ice velocities.The model includes a sub-glacial hydrology model (Johnson, 1994) that transports melt water under the ice sheet according to prevailing pressure potentials, governed by ice sheet thickness and basal topography.The response of the solid earth is modelled in a simplified way using a hydrostatically supported elastic plate model, adequate for the purpose of placing the ice sheet surface at an appropriate altitude, and hence obtaining an appropriate air temperature for the ice surface mass balance calculation.The isostatic response of the earth is modelled in a different way in UMISM than in the GIA model we use in this study, which leads to a slight inconsistency in modelling approach.Investigating the degree to which this affects the modelled ice sheet configuration is beyond the scope of this study.
The simulation is run using a palaeo-temperature record, from which the spatial pattern of air temperature is obtained, and precipitation is calculated through a mass balance relationship, developed from the Antarctic Ice Sheet (Fastook and Prentice, 1994).This precipitation is further dependent on distance from the pole, saturation vapour pressure (function of altitude and lapse rate) and surface slope.As a proxy for the air temperature record, data for the last 120 kyr from the Greenland Ice Core Project (Dansgaard et al., 1993) have been used.The ice sheet reconstruction was calibrated against dated ice-marginal positions for Weichselian stadials (e.g.Lokrantz and Sohlenius, 2006) by making slight systematic adjustments to the temperature curve (Näslund, 2010).This ice sheet calibration process did not focus on the northern ice sheet margins, resulting in a major uncertainty in ice margin position around e.g. the Barents Sea.For the basal boundary conditions the ETOPO2 digital elevation model was used as well as the geothermal heat flux model by Näslund et al. (2005).A sea level curve estimated from a previous UMISM reconstruction of all Northern Hemisphere ice sheets, representing a 100 m lowering of global mean sea level at LGM, was used, together with the prevailing bathymetry, including the location of the continental shelf, for determining the changes in position of the ice sheet grounding line for marine parts of the ice sheet.

ICE-5G (VM2)
The global ICE-5G model (Peltier, 2004) is built upon successive refinements of models of the last Pleistocene deglaciation.The initial model in the suite, ICE-1 (Peltier and Andrews, 1976), tabulated ice thicknesses of the Laurentide, Greenland and Fennoscandian ice sheets from 18 kyr BP and onward.Updates include the widely used ICE-3G (Tushingham and Peltier, 1991) and ICE-4G (Peltier, 1994) models.A new model (ICE-6G) is under development (Toscano et al., 2011) but has so far not been made publicly available.
The ICE-x suite consists of global models based on dated observations of ice sheet margins, RSL curves and the global mean sea level curve.As such, the ICE-x models critically depend on the use of the sea level equation to compute RSL estimates.ICE-1 was based on analytical relations between the distance from the ice margin and the ice thickness, assuming dynamical equilibrium of the ice sheet, as well as estimates of the ice history in some central areas considered critical.In later versions the ice thicknesses have been manually adjusted to optimise the fit to the growing body of observational constraints.Of the individual ice sheets, Antarctica is the least well constrained in that no or very little adequate data are available from this continent (Peltier, 1998).Therefore, Antarctica has mainly been used as a buffer to ensure that the fit to the global mean sea level is maintained, as well as the fit to the sparse sample of RSL records from the Southern Ocean.

P. Schmidt et al.: Weichselian Ice Sheet reconstructions
Parallel to the ICE-x development, the VMn mantle viscosity models have been developed using GIA modelling constrained by RSL data, rebound relaxation spectra, earth rotation anomalies and polar wander (Peltier and Jiang, 1996;Peltier, 1998).In the inversion for the VMn models, the ICEx models have been used as predefined loading.The derived viscosity models have then been used in constructing the next-generation ICE-x models.
The latest published version of the ICE-x suite is the ICE-5G (VM2) model (Peltier, 2004).The theoretical framework and methodology of ICE-5G is the same as that employed for its closest predecessors, but the viscosity structure of the earth model has been updated to the more advanced VM2 model (Fig. 1), which was constructed based on the ICE-4G model (Peltier and Jiang, 1996).In the original VM2 model the thickness of the lithosphere was prescribed to 120.6 km; this was however reduced to 90 km in Peltier et al. (2002) to better fit GIA data from the British Isles.As displayed in Fig. 1 the VM2 model resolves the mantle viscosity in several layers; however, Zhao et al. (2012) found that the GPS in Fennoscandia alone can only resolve three layers: the lithosphere, the upper mantle and the lower mantle.Similarly, Paulson et al. (2007) found that gravity data (GRACE) and RSL data from Hudson Bay, Canada, can also only resolve two mantle layers.We therefore note that VM2 has a mean viscosity of about 5 × 10 20 Pa s in the upper mantle and about 1.6 × 10 21 Pa s in the uppermost part of the lower mantle (670-2000 km depth).Peltier and Fairbanks (2006) presented an extension of ICE-5G to 120 kyr BP, based on the assumption of a constant areal coverage up until the LGM and scaling the pre-LGM ice thickness by the SPECMAP δ 18 O record by Martinson et al. (1987).In this study we use the extended version of ICE-5G (VM2), which is sampled on an approximately 0.7 • × 0.7 • grid with a temporal resolution of 500 yr from 17 kyr BP to present, 1 kyr between 32 and 17 kyr BP and 2 kyr at earlier times.For simplicity we will in what follows refer to this model as ICE-5G rather than ICE-5G (VM2).

ANU
The ANU model, also known as RSES, is best considered a collection of models of individual ice sheets, together comprising a global model.As in the case of the ICE-x suite, the ANU model has been developed in a series of papers, starting with Nakada and Lambeck (1987, 1988, 1989).However, in contrast to the ICE-x models, where the entire global model is updated, ANU has evolved from successive reconstructions of individual ice sheets.
The first version used ICE-1 and ICE-2 for the Laurentide and Fennoscandian ice sheets, constrained in the Barents and Kara Sea by the model by Hughes et al. (1981), and a model of the deglaciation of Antarctica were constructed based on the work by Hughes et al. (1981), Drewry (1982) and Wu and Peltier (1983).A regional model of the British Ice Sheet was added by Lambeck (1993Lambeck ( , 1995)), and the Fennoscandian Ice Sheet was modified in Lambeck et al. (1998).
The latest ANU model was presented in Lambeck et al. (2010).This revision present a new reconstruction of the Fennoscandian Ice Sheet with exception for the Barents and Kara Sea region, where the solution from Lambeck (1996) is used.A new reconstruction of the Laurentide Ice Sheet has also been generated (Lambeck et al., 2010) although not yet published.For the period preceding the LGM, the reconstruction is mainly controlled by available data on ice sheet margins and an assumption of ice sheet basal conditions equal to those at the LGM.At times prior to 64 kyr BP the reconstruction by Lambeck et al. (2006) is used.Loading of icedammed lakes and marine limit data have been added to the computation and the density structure, and elastic parameters of the earth are adopted from PREM (Preliminary Reference Earth Model Dziewonski and Anderson, 1981).
In the reconstruction the sea level equation is solved, constrained by geological markers of ice sheet extent and RSL data.As a starting model, ice thicknesses are computed from simple glaciological assumptions leading to analytical expressions for the relation between the thickness, distance from the margin and basal shear stress (Paterson, 1994), with basal stress determined from the reconstruction between 23 and 21 kyr BP in Lambeck et al. (2006).The final solution is obtained through a series of iterations involving fit to different parts of the constraining data set or introduction of new data while optimising either via a spatially and temporally varying scale factor or via the earth model parameters.Therefore, in addition to the ice sheet the reconstruction also produces an estimate of the thickness of the lithosphere (65-100 km) and the viscosity of the upper and lower mantle beneath Fennoscandia (3-4 × 10 20 and 5-20 × 10 21 Pa s respectively, Fig. 1).This can be compared to the ICE-5G reconstruction, where the earth structure is assumed known prior to the reconstruction.The spatial resolution of the model is 0.5 • in longitude and 0.25 • in latitude.In time the model is sampled on varying length intervals (450-5000 yrs), capturing the timing of important changes in the evolution of the ice sheet.

Comparison of the ice sheet reconstructions
In this section we compare the three ice sheet reconstructions directly, first in terms of integrated quantities such as volume, area and mean thickness, after which we look more closely at the details of the reconstructions at the LGM and a few selected post-LGM snapshots.We will compare the ice sheet reconstructions from 69 kyr BP onward since the UMISM model is less well constrained at prior times.

Ice volumes, areal extent and thickness
A comparison of ice volume, areal extent, mean and maximum thickness of the three ice sheet reconstructions is shown in Fig. 2. In terms of volume and areal extent we see that both ANU and UMISM display a period of small ice cover preceding the LGM by some 13-16 kyr, while ICE-5G only displays minor fluctuations in volume and nearly constant areal extent pre-LGM.The LGM in ICE-5G occurs at 26 kyr BP although the decline in volume up until 21 kyr BP is only about 7 % while the areal extent over the same period is more or less constant.In ANU and UMISM the LGM occurs at 21 and 18.2 kyr BP respectively, but we note that the maximum areal extent occurs slightly earlier (21 533 and 18 400 yr BP respectively).In general, UMISM has the smallest areal extent and ICE-5G the largest.During its two periods of extensive ice sheets, UMISM displays the greatest mean thickness of the three ice models.As seen in Fig. 2, the greatest thickness in UMISM occurs some 3.9 kyr after the LGM.
A notable difference between the reconstructions during deglaciation is a 2 kyr long hiatus in UMISM some 1.5 kyr after the LGM.This is followed by very rapid deglaciation until the Younger Dryas, at about 13 kyr BP, when the ice sheet increases slightly both in volume and extent (Fig. 2).Neither ICE-5G nor ANU displays the hiatus or the growth seen in UMISM, although the deglaciation rate is notably impeded in ANU from about 2 kyr prior to the Younger Dryas.In ANU the deglaciation rate increases again after the Younger Dryas, while in ICE-5G the deglaciation continues at an approximately unchanged rate until the end of deglaciation.

Snapshots of the Weichselian Ice Sheet reconstructions at the
LGM can be seen in Fig. 3.We see that the maximum thickness in UMISM is centred over the Gulf of Bothnia, whereas the maximum thickness in ANU and ICE-5G, at their respective LGM, is located slightly further south.Both ANU and ICE-5G display double ice domes over Fennoscandia, with one of the domes approximately co-located with the presentday centre of uplift.In ICE-5G the second dome is located over the northern part of the Gulf of Bothnia, while in ANU it is located just north-west of lake Vänern, close to the border between Norway and Sweden.To the north, ANU and ICE-5G show extensive ice sheets over the Barents and Kara seas, whereas UMISM mainly shows ice coverage in an area east of Svalbard.
From Norway toward the British Isles, both ICE-5G and UMISM have a connecting ice bridge, albeit thinner and wider in the case of ICE-5G.In ANU no such feature is seen in Fig. 3, although a continuous ice sheet from Norway to the British Isles exists in ANU between 29.5 and 27 kyr BP.A notable thinning of ice in ICE-5G is seen across a NE-SW divide running parallel to the southern shoreline of Finland down through south-central Sweden.Such a feature is absent in ANU and UMISM.Southeast of the divide, ICE-5G is significantly thinner than ANU and UMISM.At the LGM the ice edge is located slightly further inland over the Baltic countries in ANU and ICE-5G than in UMISM, and west of Norway, ICE-5G extends further out with greater thicknesses than ANU and UMISM.

Post-LGM ice sheets
The post-LGM ice sheets in the three reconstructions can be seen at three selected times in Fig. 4.Over Fennoscandia the areal extents of the ice sheets agree well, although the timing of the end of glaciation differs slightly, with ice-free conditions at 10 100 BP in UMISM, 9650 BP in ANU and 8 kyr BP in ICE-5G.The non-zero thickness of ICE-5G seen in Fig. 2 from the end of glaciation and until present day is associated with a non-vanishing small ice cover in the northernmost part of Novaya Zemlya, Russia.Similar to the ice sheets at or just prior to the LGM, we see that ICE-5G and UMISM have the greatest thickness over the Gulf of Bothnia, whereas ANU displays two local maxima further south.We also note that in the last stages of ice retreat, the ice has migrated into the mountains in UMISM, while in ICE-5G it is centred on the Bay of Bothnia and in ANU over inland northern Sweden.

Glacial isostatic adjustment modelling
To model the GIA process we need two components: an ice model and an earth model.Collectively we refer to these as a GIA model.Combined with observational data, GIA modelling makes it possible to infer rheological properties of the earth, such as the thickness of the lithosphere and the viscosity structure of the mantle.However, since many ice models require an earth model for the reconstruction of ice thickness, the two are often not independent.

Earth model implementation
Our GIA model is implemented in the commercially available finite element (FE) code Abaqus, following the recipes of Wu (2004) and Schmidt et al. (2012).We use the incompressible flat-earth approximation and neglect selfgravitation.Our implementation was benchmarked in Spada et al. (2011), where it was found that the vertical displacements and displacement rates agreed well with those of spherical models including self-gravitation, thus confirming the findings by the previous study by Schotman et al. (2008).However, Schotman et al. (2008) also found that the horizontal displacement rates predicted by the flat-earth FE model were generally larger than those predicted by an incompressible self-gravitating spectral model, unless material compressibility was included in the FE model.We therefore include material compressibility in our model.However, as we have not yet properly benchmarked our implementation with respect to the predicted horizontal displacements and displacement rates, we can not use these and will therefore The central part of our FE model covers the formerly glaciated Fennoscandian and Barents Sea regions, with a resolution of 50 ×50 km horizontally, identical to that of UMISM.The sub-surface of the model is expanded to a halfsphere of radius about 10 times the central region using a coarser mesh, and at the outer edges we apply semi-infinite elements as outer boundary conditions.Material boundaries in terms of density and elastic parameters are included at 15 and 50 km depth, at the base of the lithosphere and at 410 and 670 km, using PREM (Dziewonski and Anderson, 1981) volume averages as summarised in Table 1.Several test models have been run in order to check that our choice of outer boundary conditions as well as our averaging of PREM do not affect the predictions of the model in the region of interest.In what follows we will only vary the thickness of the lithosphere in the range 120-160 km and the viscosities of the upper and lower mantle in ranges 1-50 × 10 20 Pa s and 10-500 × 10 20 Pa s respectively.

Ice sheet implementation
The ice sheet over Fennoscandia, the British Isles and the Barents and Kara seas is implemented as a pressure source in our GIA models.We transfer the spatial sampling of ICE-5G and ANU to that of the earth model by bilinear interpolation.Thus the earth models loaded by different ice models will use the same mesh, and for earth models with identical layering and material parameters the differences in predicted displacements will only depend on differences in the load history.
Each ice model is linearly ramped up from zero thickness at 69 kyr BP to the thicknesses in the reconstructions at the closest available snapshots in time after 69 kyr BP.During the simulations we assume linear growth and decay of ice thicknesses between adjacent times, while the areal extent is assumed constant at the value of the closest preceding snapshot.This assumption leads to a smooth evolution of the ice sheet.To test whether or not excluding the ice history prior to 69 kyr BP affects the postglacial uplift, we have run two suites of test models.In the first suite we run a series of models that only include the ice sheet reconstructions between present day and 21 kyr BP, 25 kyr BP, 30 kyr BP, . .., 60 kyr BP.In each run we assume initially ice-free conditions and linearly ramp up the ice thicknesses from 0 thickness to that given by each reconstruction over a period of 1 kyr.In the second suite we assume an earth model in (close to) isostatic equilibrium at the LGM of respective reconstruction.Note that we here use 21 kyr BP as the timing of the LGM in ICE-5G.To achieve (close to) isostatic conditions, the earth models are loaded with the LGM ice sheets already at 68 kyr BP (linearly increased from 0 at 69 kyr BP), after which the ice sheets are held constant until the onset of deglaciation at the LGM.Both test suites are combined with earth models with a 120 km lithospheric thickness and uniform mantel viscosities in the range 10-60 × 10 20 Pa s.

Observational data
We compare the predicted displacement rates from our GIA models to present-day Bifrost GPS data (Fig. 5) processed by Lidberg et al. (2010) in the ITRF2005 reference system (Altamimi et al., 2007).The objective of comparing the predicted uplift rates to observational data rather than performing a direct inter-comparison between the predictions of the different reconstructions is twofold.Firstly this allows quantifying the uplift rates of each choice of earth model parameters by a single number, hence making it possible to easily present and analyse the response of a large range of earth models.Secondly this allows for a reasonable criteria to select a limited subset of earth models for which a more detailed comparison of uplift histories and residual velocities may be of interest.
The uncertainties of the vertical component in the Lidberg et al. (2010) data are in the range 0.15-0.83mm yr −1 (mean 0.34 mm yr −1 ).A prior processing of the Bifrost data in the ITRF2000 reference frame (Altamimi et al., 2002) was presented by Lidberg et al. (2007).The two realisations differ by up to approximately 1 mm yr −1 in the vertical component, with generally greater velocities in the more recent processing.As the major difference between the two solutions is the choice of reference frame (ITRF2000 vs. ITRF2005), this indicates that the uncertainties in reality are of the order of 1 mm yr −1 due to the reference frame realisation (Lidberg et al., 2010).We also note that previous studies (e.g.Wu et al., 2010;Zhao et al., 2012) have used significantly larger uncertainties than the ones presented in Lidberg et al. (2010).Here we will simply use the uncertainties given in Lidberg et al. (2010) as our aim in using the Bifrost data is not to define optimal ranges for the earth model parameters but as a familiar mean of characterising the predicted uplift rates by a single number, allowing for a simple visualisation and interpretation in the earth model parameter space.Based on the fit in parameter space we can then single out a limited subset of earth models for which more detailed comparison of predicted uplifts and uplift rates are of interest.
The fit of the model predictions to the GPS data is computed using the normalised chi-squared value where v mod i and v obs i are the vertical velocities predicted by the model and observed by GPS respectively, σ i are the uncertainties of the observed velocities, N is the number of data points and M the number of free parameters in the GIA model (in this study the lithospheric thickness and the mantle viscosity).

Impact of the pre-LGM ice history on present day uplift rates
To evaluate the impact of the pre-LGM ice history on the present-day uplift rates, we compute the χ 2 v fit (Eq. 1) between uplift rates computed using perturbed ice reconstructions (described in section 4.2) and uplift rates computed using the full reconstructions back to 69 kyr BP.We use uplift velocities predicted at the location of the GPS stations in the Bifrost network (see Fig. 5) assuming σ = 0.34, equal to the mean vertical uncertainty in the Bifrost data (Lidberg et al., 2010).A perturbation to an ice history resulting in a χ 2 v of 1 or less can then be interpreted as not resolvable by the current Bifrost data.For simplicity we only use earth models with uniform mantle viscosity.The results are summarised in Fig. 6.
For a mantle viscosity of 10 21 Pa s we find that the ice history prior to 25 kyr BP in all the reconstructions does not affect the current uplift rates enough for a χ 2 v of more than 1.If instead a mantle viscosity of 6 × 10 21 Pa s is used, changes in the ice history more recent than 40 kyr BP for UMISM and 50 kyr BP for ANU or ICE-5G may result in χ 2 v greater than 1.From Fig. 6b we find that the assumption of isostatic equilibrium at the LGM is not valid for ANU or UMISM if the mantle viscosity is 10 21 Pa s or greater.For ICE-5G the assumption can be accepted for mantle viscosities up to approximately 3 × 10 21 Pa s.

Predicted present-day uplift rates
Figure 7 displays the χ 2 v fit of the uplift rates predicted by each reconstruction to GPS data as a function of viscosity of the upper and lower mantle for lithospheric thickness in the range 120-160 km.We find that only the misfit plots of ANU display a simple region of well-fit earth model parameters.For UMISM and ICE-5G the misfit plots instead display a doughnut-shaped region of good fit, with poorer fit at the centre and outside this region.This is a bifurcation in the optimal viscosities in the sense that, given a fixed viscosity in the lower mantle, the fit to the GPS data may be optimised for two different upper-mantle viscosities, and vice versa.
The best fit among the tested earth models for each reconstruction are presented in Table 2 and marked by yellow circles in Fig. 7.We note that for ANU the misfits decreases with decreasing lithospheric thickness.For UMISM and ICE-5G the minimum misfit changes by very little over the tested range of lithospheric thickness, but we observe that the viscosity bifurcation becomes more pronounced with decreasing lithospheric thickness.
Analysing the predicted velocities we find that earth models with viscosities of about 10 21 Pa s in the upper mantle and about 7-10 × 10 21 Pa s in the lower mantle predict the greatest present-day uplift rates.In addition we see that the uplift rates generally increase with decreasing lithospheric thickness.We further find that for identical earth models, ICE-5G predicts the highest present-day uplift rates over the previously glaciated regions while ANU predicts the lowest.

Residual velocities
Figure 8 shows the residual velocities after subtraction of the uplift rates predicted by the best-fit earth models for each ice sheet from the observations.We find that ANU underpredicts the uplift rates at most stations with the exception of stations in central Sweden and Northern Denmark.UMISM over-predicts the velocities in an approximately east-westoriented band stretching from southernmost Finland and the Baltic states through southern Sweden, southernmost Norway and Denmark, with the greatest over-prediction occurring over eastern Denmark-western Sweden.In the northern parts of Norway, Sweden and Finland UMISM underpredicts the uplift rate.ICE-5G is found to over-predict the velocity over more or less all of Finland, southern Sweden and Denmark, while under-predicting the velocities in central to northern Sweden and more or less all of Norway.Common to all reconstructions is an under-prediction of the velocities at the stations south of Denmark and in the British Isles.
Given the viscosity bifurcation seen in misfit plots of UMISM and ICE-5G the residual velocities displayed in Fig. 8 may not be representative for the full parameter range of well-fit earth models.We therefore select a few more well-fitting earth models to study for these two reconstructions (see Fig. 7). Figure 9 displays the residual velocities of UMISM and ICE-5G for three different, but to GPS data well fit, earth models.Although the details differ at individual stations, the general trends observed for the best-fit earth models in Fig. 8 are also seen for each of the earth model shown in Fig. 9.In fact, extending the selection of models to all well fit models (not shown here) confirms that the observed trends are robust features of the reconstructions as implemented in our GIA model.v fit to the vertical displacement rates of the Bifrost data for GIA models using the ANU (upper row), ICE-5G (middle row) or UMISM (lower row) ice sheet reconstructions, as a function of upper-mantle and lower-mantle viscosities (η um and η lm ) and thickness of the lithosphere of 120 km (left column), 140 km (middle column) and 160 km (right column).Tiny "x"s mark the locations of tested parameter combinations, yellow circles mark the locations of best-fit models and yellow triangles mark the location of the UMISM and ICE-5G models compared in Fig. 9.

Postglacial uplift curves
Sweden (see Fig. 5).For the best-fit models we find that UMISM and ICE-5G predict comparable uplift curves at both Tromsö and Ångermanälven, while at Blekinge the uplift curve predicted by UMISM is almost twice as steep as the uplift curves predicted by ICE-5G and ANU.The uplift curves of the best-fit ANU model are further found to be significantly steeper at Ångermanälven but less step at Tromsö than those of UMISM and ICE-5G.
Comparing the uplift curves for similar earth models loaded by UMISM and ICE-5G (Fig. 10d-e) we find that at Ångermanälven similar earth models give rise to similar uplift curves generally differing by less than 10 m.At Tromsö and at Blekinge on the other hand the uplift curves generally differ by more than 15 m.ICE-5G further predicts steeper uplift curves at Tromsö than UMISM, whereas at Blekinge, UMISM predicts steeper uplift curves than ICE-5G, given similar earth models.Interestingly, the uplift curve of ICE-5G steepens with increasing upper-mantle viscosity at Blekinge.This is in contrast with the behaviour of the uplift curve at all other sites for both ICE-5G and UMISM as well as for UMISM at Blekinge (note that the apparent steepening of UMISM at Tromsö in Fig. 10d vanishes if the viscosity of the lower mantle is increased to 80 × 10 20 Pa s.).

Summary and discussion
The three ice sheet reconstructions studied here display both differences and similarities.Common to all three models are the use of dated ice-marginal positions in constraining the extent of the ice sheet at different times; however the data sets used only partly overlap.Whereas UMISM is driven by palaeo-climatic data and focuses on physically viable ice sheet dynamics, ANU and ICE-5G focus on matching post-LGM RSL data by adjusting the thickness of the ice sheet.Although both ANU and ICE-5G are based on the solution of the sea level equation, differences, especially the pre-LGM constraints, have resulted in significant differences between the two reconstructions at pre-LGM times.As both of these reconstructions have been optimised against RSL data, this indicates that RSL data alone may not be sufficient to constrain the earlier history of the Weichselian Ice Sheet.
The importance of the pre-LGM ice sheet on the presentday uplift rates increases with increasing mantle viscosity but also depends on the details of the reconstruction (Fig. 6).Whereas the assumption of isostatic equilibrium at the LGM may be reasonable for ICE-5G with its massive pre-LGM ice sheet, it will not be valid for UMISM or ANU, where the ice sheet grows from almost ice-free conditions at about 34 kyr BP up until the LGM (Fig. 2) unless the mantle viscosity is less than 10 21 Pa s.To investigate the response of a mantle with a viscosity of 6 × 10 21 Pa s, the ice history from at least 50 kyr BP needs to be taken into account, for all reconstructions (Fig. 6a).The start of our simulations at 69 kyr BP is therefore sufficient for the range of viscosities studied here.
The style of deglaciation in UMISM differs notably from ANU and ICE-5G, Fig. 2. Specifically, the initially increasing mean thickness during periods of deglaciation indicates that the ice sheet melts preferentially from the edges inwards in UMISM, whereas the correlation between mean thickness, ice volume and extent in ANU and ICE-5G indicates that the ice melts more or less uniformly over the ice sheet in these reconstructions.From about 40 kyr BP up until the LGM, UMISM and ANU display similar evolutions in terms of integrated characteristics such as volume and areal extent (Fig. 2).However, the fit of the uplift rates predicted by UMISM to the Bifrost data is more similar to that of ICE-5G.This is mainly caused by the about 2.8 kyr later LGM in P. Schmidt et al.: Weichselian Ice Sheet reconstructions UMISM than in ANU, resulting in greater present-day uplift velocities in UMISM and therefore predictions closer to those of the more massive ICE-5G.
Despite a fundamentally different approach used in the reconstruction, UMISM can fit present-day uplift rates in Fennoscandia as well as both ANU and ICE-5G.It is therefore not feasible to claim one reconstruction to be more successful than the others in reproducing the GIA observations we have used herein.However, based primarily on physical principles rather than inversion of GIA data, the UMISM reconstruction is bound to evolve as a real ice sheet (within the limits of the modelled physics), which is not the case for the ICE-5G or the ANU reconstructions.In addition, the modelling of the earth response is handled in a simplified way in UMISM.This indicates a slightly lower degree of coupling between the earth and the ice model than offered by ICE-5G and ANU, as the reconstruction in these latter two models largely rests on a more realistic modelling of the earth response.
The viscosity bifurcation seen in the misfit plots of ICE-5G and UMISM (Fig. 7) are not unique to this study.A similar feature can be seen in the misfit plots in Steffen et al. (2010) and Lidberg et al. (2010), as well as possibly hinted at in the misfit plots in Milne et al. (2004).The doughnutshaped region of well-fitting models can be readily understood from a simple model of postglacial uplift.During the rebound process, the vertical displacement, w, in a formerly glaciated region can be described by a function on the form (e.g.Turcotte and Schubert, 2002) where W is a constant proportional to the maximum depression and A is a site-specific constant.The vertical displacement rate is then given by a function of the form This function has a maximum velocity, v * (t), at a viscosity of The equations above show that at any given time, t, there exists a viscosity, η * (t), that will give rise to the greatest uplift rate, v * (t).If however the actual viscosity is either higher or lower than η * (t), the uplift rate will be smaller than v * (t).
The time dependence of η * (t) is a reflection of the fact that a low-viscosity earth will rebound fast, with the rebound velocity decaying fast with time.In a high-viscosity earth the initial rebound velocity is small, but the decay with time is slow.With time, the rebound velocity in a low-viscosity earth will therefore become smaller than the rebound viscosity in a high-viscosity earth, subjected to an identical loading history.In a multilayer model each viscous layer will have its own v * (t) and corresponding η * (t), potentially giving rise to a bifurcation in the optimal viscosity of each layer as seen in the ICE-5G and UMISM misfit plots.In the ANU case the predicted present-day uplift rates are very close to v * (now).
The above analyse is valid for all stations well within the region that was once covered by the Weichselian Ice Sheet.In fact, plotting the misfit at single stations yields the same viscosity bifurcation as seen in Fig. 7 for all stations located well inside the LGM extent of the ice sheets, whereas at stations close to or outside the LGM margins display a more complex misfit field.As the uplift rate is of greater magnitude beneath the former ice cover, the misfit at stations located here will also dominate the total misfit plot.Given identical earth models, ICE-5G in general predicts the fastest present-day uplift rates and ANU the slowest.Analysing the residual velocities for well-fitting earth models, with emphasis on trends independent of the earth model, we find that improvements can be made in all three reconstructions.In particular, we find that ICE-5G tends to overpredict the velocities over Finland whereas ANU tends to under-predict them (Fig. 8 and 9).Although inspection of the respective ice sheet thicknesses shows that the post-LGM ice sheet in ANU is relatively thin over Finland while ICE-5G displays a thick ice coverage over Finland, stretching well into westernmost Russia (Fig. 4), this may not be the sole explanation.Also the post-LGM ice sheet in UMISM is relatively thick over Finland and western most Russia (Fig. 4), yet the predicted uplift velocities in Finland lies in between those predicted by ANU and ICE-5G (Fig. 8).Common to UMISM and ANU is instead a continuation of the ice thickness south of Finland, whereas in ICE-5G a clear divide is seen along the southern border of Finland, with great ice thicknesses to the north and thin ice to the south.Further, the centre of the ice sheet in both UMISM and ANU migrates westward in the last stages of the deglaciation phase, whilst in ICE-5G the centre stays more or less fixed over the Gulf of Bothnia.It is therefore likely that ANU would benefit from greater post-LGM thickness to the east whereas ICE-5G would benefit by slightly reduced thicknesses over inland Finland, greater post-LGM thickness over the southern half of the Baltic Sea and the western shores of the Baltic states, as well as a westward migration of the ice centre in the final stages of deglaciation.UMISM generally over-predicts the velocities over Denmark and southern to central Sweden; similar trends can also be seen in ANU and ICE-5G.This may indicate that the centre of mass in all three reconstructions at the LGM and onward is placed slightly too far south.
The trends we have identified here apply to a wide range of earth model parameters that yields uplift rates in reasonable agreement with observed present-day uplift rates in Fennoscandia and are therefore robust features of the ice sheet reconstructions as implemented in our GIA model.However, it should be noted that we do not include the ocean load from the melting glaciers, and this may affect our conclusions.Further, several studies using different methods have shown that the lithospheric thickness changes from great thicknesses in the old cratonic parts in the east to a thinner lithosphere in the younger western part (e.g.Pérez-Gussinyé and Watts, 2005; Priestley and M c Kenzie, 2006;  Artemieva and Thybo, 2008).It is therefore to be expected that also the GIA process will "see" an approximately eastwest variation of the lithospheric thickness.Whitehouse et al. (2006) showed using the previous version of ANU (Lambeck et al., 1998) that this will result in reduced uplift rates over Finland and central Sweden, while the uplift rates increase in northern Finland, Norway, Denmark and the Baltic states.Although not shown here we can confirm these general trends for all three reconstructions used in this study (for UMISM see Lund et al., 2009), but we note that the effect diminishes with increasing mantle viscosity.Thus a laterally varying lithosphere may at a first glance seem to counteract the over-predictions seen for ICE-5G in Finland and the under-predictions in central Norway; however it would enhance the over-predictions seen in northern Finland and northern Norway as well as in Denmark.Similarly, regions can be found for both UMISM and ANU where a lateral variation in the lithospheric thickness would either improve or degrade the fit of the predicted uplift rates.Lateral variations in the lithospheric may therefore also affect our suggested improvements to the ice sheet reconstructions.
While the fit of the predicted uplift rates to the Bifrost data fails to single out a "best" ice reconstruction, comparison of predicted uplift curves indicates that inclusion of data measuring the past uplift (e.g.RSL data) may resolve this (Fig. 10).Further, comparing the uplift curves of similar well-fit earth models for UMISM and ICE-5G indicates that also the viscosity bifurcation observed in the uplift misfits may be resolved by inclusion of such data.Unfortunately, as our current GIA model does not solve the sea level equation nor compute geoid heights, this is out of the scope of this study.We finally note that the flatter uplift curves at Tromsö and steeper uplift curves at Blekinge by UMISM compared ICE-5G correlate well with differences in the ice sheet from the LGM and on.As seen in Fig. 3 and 4, ICE-5G has a more extensive ice sheet in the Tromsö region and a thinner ice sheet in the Blekinge region than UMISM.

Conclusions
We have compared a thermo-mechanical ice sheet reconstruction, UMISM, to two reconstructions based on the relative sea equation, ANU and ICE-5G, and commonly used in GIA modelling.We find that given appropriate earth models the predictions by UMISM fit observational data equally well as the predictions by ANU and ICE-5G.Although it is not possible to claim one model better than the others in reproducing the present-day uplift rates, UMISM has the benefit of being based on thermo-mechanical modelling and is therefore a more physically viable reconstruction than ANU or ICE-5G, where the ice sheet thickness has been optimised to yield postglacial uplift curves in accordance with observations.We find that while the characteristics of the UMISM reconstruction is more similar to ANU, the overall fit of the predicted present-day uplift rates to the Bifrost data is more similar to ICE-5G.ANU yield a relatively well constrained best-fit model when compared to observed present-day uplift rates.However, both ICE-5G and ANU bifurcations in the optimal upper-and lower-mantle viscosity give rise to a range of well-fitting models.In general, though, with identical earth models, the present-day uplift velocities predicted by UMISM are smaller than those predicted by ICE-5G but greater than ANU predictions.Given the large difference between ICE-5G and ANU ice sheets it is clear that there exists a large freedom in reconstructing the Weichselian Ice Sheet based on an optimisation to observational RSL data.We find that improvements can be made to all three reconstructions to better fit the observed present-day uplift rates.More specifically, the post-LGM ice sheet in ANU would benefit from greater thickness to the east, while the post-LGM ice sheet in ICE-5G would benefit from reduced thickness over Finland, increased thickness over the southern half of the Baltic Sea and the western shores of the Baltic states and a westward migration of the ice centre in the final stages of deglaciation.In addition, in all three reconstructions, the mass centre of the ice sheet appears to be located slightly too far south from the LGM and onward.These suggested improvements may however be affected by our neglect of the ocean load as well as our assumption of a uniform lithospheric thickness.

Figure 1 .
Figure 1.Viscosity profile VM2 (black line) used in the ICE-5G reconstruction and optimal two-layer viscosity range (grey regions) found in the ANU reconstruction.

Figure 2 .
Figure 2. Temporal development of (a) ice volume, (b) areal extent, (c) mean thickness and (d) maximum thickness of the UMISM (red), ICE-5G (green) and ANU (black) ice sheet reconstructions as implemented in this study.The simulations are started by linearly increasing the load from 0 at 69 kyr BP to the snapshot of each reconstruction closest in time as shown by the unfilled circles in the figure.The non-zero thickness in ICE-5G from about 8 kyr BP until today is associated with a non-vanishing small ice cover in the northernmost part of Nova Zemlya, Russia.

Figure 3 .
Figure 3. Ice sheet extent and thickness at 26 kyr BP (a-c, LGM of ICE-5G), 21 kyr BP (d-f, LGM of ANU) and 18 200 BP (g-i, LGM of UMISM) for ICE-5G (a, d and g), ANU (b, e and h) and the UMISM model (c, f and i).Bold font header indicates the LGM of respective ice sheet reconstruction (a, e and i).ICE-5G and ANU ice thicknesses at 18.2 kyr BP have been linearly interpolated from adjacent time frames (16.5 and 20 kyr BP for ANU, 18 and 19 kyr BP for ICE-5G), while the extents have been inherited from the closest preceding snapshot.

Figure 4 .
Figure 4. Ice sheet extent and thickness at 16.5 kyr BP (a-c), 14 kyr BP (d-f) and 10.5 kyr BP (g-i) for ICE-5G (a, d and g), ANU (b, e and h) and the UMISM model (c, f and i).Note that the ANU model at 10.5 kyr BP has been generated by linear interpolation in thickness between the solution at 10 910 and 10 274 yr BP, the extent of the ice sheet is inherited from 10 910 yr BP, while the ice sheet displayed at 14 kyr BP is the snapshot at 13 940 yr BP.

Figure 5 .
Figure 5. Present-day uplift rates in Fennoscandia as measured by GPS in the Bifrost project as processed by Lidberg et al. (2010).Triangles mark the location of the GPS stations; circles abbreviated with a, b and c show locations at which predicted uplift curves are displayed in Sect.5.4

Figure 10 Figure 6 .
Figure 10 displays uplift curves since 10 kyr BP for all GIA models presented in Figs. 8 and 9 at three selected sites from Tromsö in northern Norway through Ångermanälven in Sweden and close to the uplift centre to Blekinge in southern

Figure 8 .Figure 9 .Figure 10 .
Figure 8. Residual vertical Bifrost velocities after subtraction of the best-fit earth models of ANU (left), UMISM (middle) and ICE-5G (right).Positive residuals (green colours) show where the reconstructions under-predicts the present-day uplift, whereas negative residuals (blue colours) show regions of over-predictions.Red triangles indicate the location of the stations in the Bifrost network.

Table 1 .
(Dziewonski and Anderson, 1981) densities used in the earth models as derived from volume averages of PREM(Dziewonski and Anderson, 1981).

Table 2 .
Parameters for the earth models that best fit the Bifrost vertical velocities, from Fig.7, for all three ice sheet reconstructions.