Introduction
Ice from the Eemian period (130–115 ka bp) was found in the central Greenland ice cores (GRIP (Greenland Icecore Project) and GISP2 (Greenland Ice Sheet Project 2)) as well as in the NorthGRIP (North Greenland Icecore Project) ice core. The early part of the Eemian layer was, however, gone in the latter due to a high basal melt rate, and in the central Greenland cores the stratigraphy was broken in the bottom 10% of the cores due to flow over an uneven bed. An undisturbed record of the full Eemian period has therefore not yet been obtained from Greenland. A new ice-core drill site, NEEM (North Eemian), has been suggested at 77.449 ˚ N, 51.056 ˚ W in northwest Greenland, and drilling is planned to begin in 2008. The NEEM drilling project is an international effort with 14 participating nations, and its main purpose is to retrieve a continuous record of the whole Eemian interglacial.
The new drill site is located 365 km downstream from NorthGRIP on the ice ridge that runs north-northwest from GRIP via NorthGRIP towards Camp Century (see Fig. 1). The altitude at NEEM is 2447 m and radar investigations indicate an ice thickness of 2561 m. The accumulation rate and surface velocity are not well known but they are expected to be somewhat larger than at NorthGRIP since NEEM is located further out on the flank at lower altitude and with steeper surface slope. From radio-echo sounding (RES) images it is seen that the very smooth bed found around NorthGRIP does not extend all the way to the NEEM drill site. However, the bedrock undulations at NEEM are on a much smaller scale than at GRIP, so there is no immediate reason to suspect folding of the deep layers. The location of the drill site for the new ice core was selected from studies of the internal structure of the ice, as seen on RES images. The Eemian layer is located too deep in the ice to show up in the existing RES images, but the shape of the younger internal layers indicates that Eemian ice is located relatively high above the bed at NEEM. A modelling effort is needed in order to predict the depth and thickness of the Eemian layer.
Fig. 1. Map of Greenland with indicated drill sites along the ice ridge in northwest Greenland. Starting from the south and moving along the flowline the drill sites are: GRIP, NorthGRIP, NEEM and Camp Century.
Radar Data
Large parts of the Greenland ice sheet have been investigated with airborne radio-echo sounders by the Center for Remote Sensing of Ice Sheets, University of Kansas, USA (Reference Chuah, Gogineni, Allen, Wohletz and LawrenceChuah and others, 1996). The RES images show the ice surface, the ice– bedrock interface and internal layers in the ice. Individual internal layers can be followed over hundreds of kilometres and are generally accepted as isochrones. The shape of the isochrones reveals information about the ice dynamics, especially the basal melt rate. Undulations which increase with depth are an indication of spatially changing basal melt rates. In an area with a high basal melt rate, the isochrones will be pulled down faster than in areas with low or no basal melting.
In order to constrain a Monte Carlo analysis of the inverse problem presented below, we need a dated set of observed isochrones. Twelve layers have been traced from NorthGRIP to NEEM in the RES images. The layers are dated from their depths in the NorthGRIP ice core, and their ages fall between 3.6 and 79.8 ka. The two deepest isochrones are very faint in the RES images, and in some areas they are impossible to trace. However, the visible parts of these two layers are included, as it is crucial to have deep constraints on the Monte Carlo analysis.
Modelling
Forward model
A 224 km long section along the north-northwest-trending ice ridge is considered, and a model approach similar to that used by Reference Buchardt and Dahl-JensenBuchardt and Dahl-Jensen (2007) is taken. The ice flow along the ice ridge is therefore simulated using a two-dimensional Dansgaard–Johnsen model (Reference Dansgaard and JohnsenDansgaard and Johnsen, 1969) that accounts for basal melting and sliding. The model requires surface velocity u sur, ice thickness H, accumulation rate a, basal melt rate w b, kink height h and the ratio of basal sliding velocity to surface velocity F B as input. The horizontal velocity u and the vertical velocity w are calculated as
and
where z is the ice equivalent height above bedrock and
All parameters are allowed to vary horizontally, whereas temporal changes are only considered for H and a.
The spatial variation in H is known from radar surveys (Reference Chuah, Gogineni, Allen, Wohletz and LawrenceChuah and others, 1996), and the temporal changes are calculated using the SICOPOLIS (SImulation COde for POLy-thermal Ice Sheets) ice-sheet model for Greenland (Reference GreveGreve, 2005). The accumulation rates in the area are not well known. Accumulation rate values along the ice ridge are inferred from Reference Ohmura and ReehOhmura and Reeh (1991) but tuned to match the known values at the drill sites. The accumulation maps from Reference Ohmura and ReehOhmura and Reeh (1991) are based on interpolations between measurements with a coarser resolution than optimal for use in this study. We therefore wish to allow for the possibility that the overall pattern with low values in the centre of the ice sheet and higher values closer to the coast is more pronounced in the area of study than suggested by the current dataset. This is done by adding to the data a contribution that grows linearly with the distance from NorthGRIP (where the accumulation rate is well determined). The speed γ with which the contribution grows with distance is left to be determined from the Monte Carlo analysis.
The accumulation history is calculated from the dated δ18O record from NorthGRIP using a model similar to that developed by Reference Johnsen, Dahl-Jensen, Dansgaard and GundestrupJohnsen and others (1995). We therefore have high accumulation in warm periods and low in cold periods. It is assumed that the ratio of the accumulation rate at any point along the line to that at NorthGRIP is constant in time. The oldest ice found in the NorthGRIP ice core has been dated to 123 ka. When starting the forward model before this time, the glacial index from Reference GreveGreve (2005) is used to scale the accumulation further back in time.
The surface velocities in the area around NEEM are not known. However, from empirical studies of the area around NorthGRIP where a strain net was established, a linear relationship between surface slope and along-ridge surface velocities has been found. Therefore, the value of u sur is in this
work calculated as a linear function of the observed surface slope, i.e.
where C is a constant and dS/dx is the surface slope. The value of C is chosen so that the calculated surface velocity at NorthGRIP matches the observed value of 1.3 ma − 1 (Reference Hvidberg, Keller and GundestrupHvidberg and others, 2002). At all times, the value of u sur calculated as above is scaled with the factor a(t)/a 0 to compensate for the changes in accumulation rate.
The value of w b is allowed to change for every 8 km, giving rise to 28 different unknown values along the section. The kink height h and the fraction of basal sliding F B are tied linearly to the melt rate to restrict the number of unknowns (Reference Buchardt and Dahl-JensenBuchardt and Dahl-Jensen, 2007).
To investigate the importance of including the temporal changes of H, the model is also run with dH/dt = 0 for comparison.
Solving the inverse problem
Thirty-four parameters from the forward model are unknown and must be determined from a Monte Carlo analysis. The forward model is run several hundred thousand times with various combinations of values for the model parameters. After each run, the fit between observed isochrones and isochrones calculated by the forward model is used to decide whether the used model parameters should be accepted or rejected. The values of the model parameters are changed between each run by using random numbers (a random walk in the model space). By examining the statistical properties of the accepted values for the model parameters, best estimates for these 34 unknown parameters can be inferred. A more thorough description of the Metropolis algorithm used to solve this problem is given by Reference Buchardt and Dahl-JensenBuchardt and Dahl-Jensen (2007).
Results And Discussion
Time-dependent ice thickness
The fit between the observed and the modelled isochrones is satisfactory. The modelled isochrones recreate the overall shape of those observed (Fig. 2). About 50 km upstream from NEEM, the observed isochrones display high-amplitude undulations. These variations happen over too short a distance to be resolved by the model, which has melt rate intervals 8 km long.
Fig. 2. Observed (solid lines) and modelled (dotted lines) isochrones. The NEEM drill site is indicated by the vertical line.
Most model parameters are well defined by the Monte Carlo analysis, i.e. the histograms of the accepted values for the parameters resemble Gaussian distributions. Only the melt rate values in the intervals where there are no data for the two deepest isochrones are not well defined. In these intervals, the distributions for the accepted values for the melt rate are broader and do not show a clear single maximum (Fig. 3). This illustrates the importance of including deep isochrones if information on the melt rate is required.
Fig. 3. Histogram of accepted values for the melt rate (a) in the interval at NEEM and (b) 75 km upstream from NEEM, where there are no data for the lowest two isochrones. In (a), the histogram shows a single maximum, so the melt rate in this interval is well defined by the solution. However, it can be seen that the solution in (b) does not provide a well-defined value for the melt rate.
The melt rate at NEEM is found to be 1.0±0.6mma − 1. However, about 50 km upstream from the drill site where the observed layers show undulations of a very high amplitude, values of almost 11 mma − 1 are found. Using the value of γ found from the Monte Carlo analysis, we find the accumulation rate at NEEM to be 0.26ma − 1.
Running the forward model with the parameter estimates found from this analysis results in a modelled Eemian layer of 60±10m thickness located 100±40m above bedrock at NEEM (Fig. 4). The age of the ice at the base is estimated to be around 200 ka. Eemian ice has been transported approximately 50 km along the ice ridge since deposition. It is therefore not likely that the Eemian layer has been significantly affected by the higher melt rates in the previously mentioned area with large undulations of the isochrones.
Fig. 4. Observed isochrones (solid lines) and modelled top and base of the Eemian layer (dashed lines) calculated from the model with time-dependent ice thickness. The dotted lines show the paths that the top and the base of the Eemian layer at NEEM have taken through the ice. The NEEM drill site is indicated by the vertical line.
Constant ice thickness
The general observations regarding fit and well-determined parameters for the model with constant ice thickness are the same as for the analysis of the model with time-dependent ice thickness. The thickness of the Eemian layer is found to be 75±10m and the base of the layer is found to be located 150±40m above bedrock. The accumulation rate in this case is found to be 0.27ma − 1 and the melt rate is 1.0±0.5mma − 1.
In the upper part of the ice sheet, solving both inverse problems results in modelled isochrones that are too shallow close to NEEM and too deep upstream from the drill site (Fig. 2). A possible cause for this is that the accumulation pattern used in the model may not be a good estimate of the true accumulation pattern. Furthermore, changes over time in the accumulation pattern are not accounted for in the model but may have significantly influenced the shape of the observed isochrones.
Conclusion
Both models indicate that a full record of the Eemian period can be obtained at NEEM, and both models predict the layer to be located well above bedrock. The model with time-dependent ice thickness predicts an Eemian layer of 60±10 m; the model with constant ice thickness over time predicts a thickness of 75±10 m. Both solutions indicate basal ice of an age of at least 200 ka. The melt rate at NEEM is estimated to be 1 mm with an uncertainty of around 0.5 mm. This is good news for the drilling project, since a small melt rate keeps the layer thickness at the base larger than if there is no melting.
