Introduction

The NorthGRIP (North Greenland Ice Core Project) ice core was drilled during 1996–2004 at 75.108 N, 42.328 W, 316km NNW of the GRIP drill site in Central Greenland. The ice in the area flows along a NNW-trending ice ridge from GRIP towards NorthGRIP. The surface velocity at NorthGRIP is 1.3 ma^–1 (Reference Hvidberg, Keller and GundestrupHvidberg and others, 2002), the ice thickness is 3090m and the present mean annual temperature is –31.58C. The aim of the drilling was to retrieve ice from the Eemian interglacial period 130–115 kyr ago. Before drilling was initiated, it was predicted that the Eemian layer would be found at depths of 2750–2850m (Reference Dahl-JensenDahl-Jensen and others, 1997). However, as bedrock was approached it became evident that the ice was melting at the bottom. The basal layers did not thin as fast as was expected, and Eemian ice was not encountered until 80 m above bedrock (North Greenland Ice Core Project members, 2004). When it had been established that there is basal melting at NorthGRIP, it was concluded from radio-echo sounding (RES) images that the ice must be melting at the base in a large area in Northern Greenland. As the geothermal heat flux in the area is unknown, it is not straightforward to calculate the melt rate at NorthGRIP. Furthermore, the shape of the RES layers suggests that the melt rate varies significantly over short distances in the area (Reference Dahl-Jensen, Gundestrup, Gogineni and MillerDahl-Jensen and others, 2003). Using a Monte Carlo method to invert an ice flow model for the flow line from GRIP to NorthGRIP allows us to estimate the melt rates in the area around NorthGRIP.

Modelling the Iceflow

A Dansgaard-Johnsen model (Reference Dansgaard and JohnsenDansgaard and Johnsen, 1969) is used to simulate the ice flow along the flowline from GRIP to NorthGRIP. Both one- and two-dimensional approaches are used. For this purpose, a coordinate system is adopted with a horizontal x axis along the NNW-trending ice ridge in the direction of the flow at NorthGRIP, and a vertical z axis pointing upwards. The origin of this coordinate system is located at GRIP at sea level. This study has been concerned with a 104 km long section of the ridge starting 82 km upstream from NorthGRIP and ending 22 km downstream. Accounting for melting and sliding at the base, and assuming constant ice thickness with time, the horizontal velocity u and the vertical velocity w are given by

(1)

and

(2)

respectively. Here u _sur is the horizontal surface velocity, z the ice equivalent height above bedrock, F _B = u _bed /u _sur the fraction of basal sliding, H the ice thickness in ice equivalent and h is called the kink height. The vertical velocity at the base is w ₀, and

(3)

where a is the annual ice equivalent accumulation. The basal melt rate w _b is given by w _b = –w ₀. The one-dimensional model is obtained by disregarding the horizontal movement, u = 0.

In the two-dimensional model, the basal melt rate is allowed to vary along the flowline, changing value every 4 km. The melt rate is considered constant within each of these 4 km intervals. In order to limit the number of parameters to be determined by the Monte Carlo inversion, the kink height h and the fraction of basal sliding F _B are considered linear functions of the melt rate:

(4)

(5)

Thus h and F _B also vary from one 4 km interval to another. Reference Dahl-JensenDahl-Jensen and others (1997) obtained estimates for the present accumulation rates along the flow line from shallow ice core studies. In the present work, it is assumed that the ratio of the accumulation rate at any point along the line to that at NorthGRIP is constant in time. Thus, the accumulation history at any point along the line can be inferred from the accumulation history at NorthGRIP aðtÞ which is calculated from the accumulation model presented below.

The ice thickness at NorthGRIP is assumed to be constant in time in agreement with model results (Letréguilley and others, 1991; Reference Marshall and CuffeyMarshall and Cuffey, 2000). All the parameters of the flow model except the accumulation rate a are thus assumed to be constant in time. α, ß, h ₀ and the value of w _b within each 4 km interval are unknown and will be estimated using a Monte Carlo inversion of the flow model.

The Accumulation Model

The accumulation history at the NorthGRIP drill site is calculated using a model of the same type as that used by Reference Johnsen, Dahl-Jensen, Dansgaard and GundestrupJohnsen and others (1995) to date the GRIP ice core. The time-dependent ice equivalent accumulation rate aðtÞ is calculated from the measured d¹⁸O values:

(6)

(7)

where a ₀ is the present ice equivalent accumulation rate at

NorthGRIP and d¹⁸O_w = –35.2% and d¹⁸O_c = –42% are typical d¹⁸O values for warm and cold climate at NorthGRIP, respectively. c ₁ and c ₂ denote the relative slopes of a in warm and cold climate, respectively, and are defined as

(8)

The parameters a ₀, c ₁ and c ₂ are found from the Monte Carlo inversion. The GICC05 timescale (Reference RasmussenRasmussen and others, 2005; Reference VintherVinther and others, 2006) is used for the d¹⁸O curve back to 42 kyr b2k (before 2000 AD), and further back in time the ss09sea timescale (Reference JohnsenJohnsen and others, 2001) is used. However, the ss09sea timescale has been shifted to agree with the GICC05 at 42 kyr b2k. The measured d¹⁸O values have been corrected for the changes in the isotopic composition of seawater due to the build up of ice on the continents during the glacial period (Reference WaelbroeckWaelbroeck and others, 2002).

Monte Carlo Inversion

In the one-dimensional model, the horizontal velocity is u = 0 and only the basal melt rate at NorthGRIP is included. Thus the kink height h and the fraction of basal sliding F _B are included directly as model parameters instead of α, h ₀ and ß (see Equations (4) and (5)). This reduces the number of model parameters to be determined by the Monte Carlo inversion to 6: c ₁, c ₂, a ₀, F _B, h and w _b.

In the two-dimensional model, the basal melt rate w _b has 26 unknown values, one for each 4km along the 104 km long flowline. Together with α, h ₀ and ß from Equations (4) and (5) and c ₁, c ₂ and a ₀ of the accumulation model it adds up to a total of 32 unknown model parameters.

An observed data set exists d ^obs consisting of 20 internal layers identified in the RES images (Reference Chuah, Gogineni, Allen and WohletzChuah and others, 1996; Reference Dahl-JensenDahl-Jensen and others, 1997; Reference Gogineni, Chuah, Allen, Jezek and MooreGogineni and others, 1998, Reference Gogineni2001; Reference Fahnestock, Abdalati, Joughin, Brozena and GogineniFahnestock and others, 2001; Reference Kanagaratnam, Gogineni, Gundestrup and LarsenKanagaratnam and others, 2001). The layers are generally accepted to be isochrones. They have been dated from their depths (600– 2700 m) in the NorthGRIP ice core using the same timescale as for the d¹⁸O record. This gives isochrone ages from 3.5– 79.6 kyr. We will now use the ice flow model and the observed data to calculate the unknown model parameters as an inverse problem. Since the problem is highly nonlinear we turn to a Monte Carlo method in order to solve it. The model space is investigated through a random walk. For each step in the random walk, a modelled data set d(m) is created by running the forward flow model with the combination m of unknown model parameters. This is compared to the observed data set by calculating the misfit function S:

(9)

where i = 1–20 as there are 20 isochrones and j runs through the 81 data points followed on each isochrone. s_ij denotes the uncertainty in a data point d ^obs _ij . This uncertainty is given by the vertical resolution of the radar used to measure data. The starting point of the forward model is 79.6 kyr ago, since we do not have older isochrones to compare. The model is run to the present time in steps of 100 years. The likelihood function L is given by

(10)

where k is a normalization constant. Each step of the random walk is accepted or rejected according to the Metropolis criterion

(11)

where m _current is the most recently accepted model and m _test is the model being tested. It can be shown that this random walk samples the posterior probability density in the model space (Reference Mosegaard and TarantolaMosegaard and Tarantola, 1995). The final result is independent of the choice of initial values for the unknown model parameters.

Results

One-dimensional inversion

The random walk in the model space was continued until a reasonable statistic was obtained. In the results presented here, 300 000 models were accepted. The distributions of the accepted values for each model parameter are shown in Figure 1. The mean and standard deviation for each distribution are displayed above the histograms. All distributions are seen to resemble Gaussian distributions, with strong single maxima. This means that the parameters are well defined by the Monte Carlo inversion.

Fig. 1.

Histograms of accepted values for parameters of the one-dimensional model; mean and standard deviation are displayed above each distribution. (a), (b) The relative slopes of the accumulation rate in warm (c ₁) and cold (c ₂) climate, respectively (Equation (8)). (c) The present ice equivalent accumulation rate at NorthGRIP a ₀. (d) The fraction of basal sliding, F _B. (e) The kink height h from the Dansgaard-Johnsen model (Equations (1) and (2)). (f) The basal melt rate at NorthGRIP, w _b.

The result for the melt rate at NorthGRIP is found to be 8.2 ±0.9 mma^–1. When the melt rate is known, the amount of heat used to melt the ice Q _melt can be calculated using the relation

(12)

where ρ and L _ice are the density and latent heat of ice, respectively. The geothermal heat flux Q _geo is given by the sum of the amount of heat used to melt the ice and the amount of heat conducted through the ice Q _ice :

(13)

Q _ice is determined from the gradient of the observed temperature profile әT/әz at the base at NorthGRIP, i.e.

(14)

where K is the thermal conductivity of ice. Using Q _ice = 70 mWm^–2 (North Greenland Ice Core Project members, 2004) and the basal melt rate found in this study, the geothermal heat flux at NorthGRIP is calculated to be 150±12mWm^–2.

Two-dimensional inversion

In this inversion, the full suite of 32 model parameters was determined. The random walk in the model space was continued until a reasonable statistic was obtained. In the results presented here, 250 000 models were accepted. The distributions of the accepted values for the model parameters are shown in Figures 2 and 3. The three parameters from the accumulation model (c ₁, c ₂ and a ₀) are all well determined by the Monte Carlo inversion (Figs 2a–c), while the distributions for α and h ₀ both show a double peak. The consequence is that the kink height h calculated from Equation (4) is not well determined by the inversion. However, the peaks are close together so the effect on the determination of the basal melt rates is small.

Fig. 2.

Histograms of accepted values for parameters of the two-dimensional model (see Fig. 3 for the melt rates); mean and standard deviation are displayed above each distribution. (a), (b) The relative slopes of the accumulation rate in warm (c ₁) and cold (c ₂) climate, respectively. (c) The present ice equivalent accumulation rate at NorthGRIP a ₀. (d), (e) The parameters h ₀ and α linking the kink height to the melt rate (Equation (4)). (f) The link ß between the fraction of bottom sliding and the basal melt rate (Equation (5)).

Fig. 3.

Histograms of accepted values of the basal melt rate at 4 km intervals along the flow line. The measurements closest to GRIP are displayed at the top left, and those furthest from GRIP are at the lower right. The histogram for the interval containing NorthGRIP is shown in black. The units are mma^–1.

The basal melt rates are well determined for all the 4 km long intervals except the first five (see Fig. 3). The effect of basal melting on the internal layers increases with depth, so the deep layers are very important for the determination of the melt rates. Due to the horizontal movement of the ice, the modelled isochrones have moved out of the first intervals before they have reached great depths. As a consequence, the inversion has not had any constraints in the deep part of the ice for the first part of the line, and the melt rate estimates obtained for that area are badly constrained. The melt rate is seen to vary between 5.3±0.2mma^–1 and 21.2±3.6mma^–1 with the smallest value just upstream from the NorthGRIP drill site. The melt rate at NorthGRIP is found to be 6.1 ±0.2mma^–1. This is considerably lower than the estimate obtained from the one-dimensional model. The higher melt rates upstream from the drill site pull the internal layers down before the ice reaches the NorthGRIP drill site. The one-dimensional model thus compensates for the upstream effect by over-estimating the melt rate.

Figure 4 shows a comparison between observed and modelled isochrones in the lower part of the ice sheet. It can be seen that the modelled isochrones successfully reproduce the large-scale variations of the observed isochrones.

Fig. 4.

(a) RES image collected along the ice ridge between GRIP and NorthGRIP. The surface and bedrock are shown in white. NorthGRIP is indicated by the vertical dotted line. (b) Comparison between observed (dotted) and modelled (solid) isochrones in the lower part of the ice sheet. The shown isochrones have been dated to 28.6, 34.6, 37.6, 44.7, 53.8, 59.7, 75.2 and 79.6 kyr b2k, respectively. The bedrock is shown in the bottom of the plot. The modelled isochrones are seen to reproduce the large-scale undulations of the observed isochrones. Higher accumulation rate and higher melt rate cause the isochrones to be located deeper in the ice sheet upstream (left) from the drill site.

Figure 5 shows a comparison between the shape of the lowest observed isochrone dated to 79.6 kyr b2k and the variation of the melt rate along the line. The two curves show very similar patterns, but the isochrone curve is shifted slightly to the right. The shift is caused by the horizontal flow velocity of the ice. The features created by the melt rate at a given place is carried with the ice along the line. This illustrates the advantage of using a two- to simulate the ice flow.

Fig. 5.

(a) The lowest isochrone from the radar data set used in this study. (b) The melt rates along the line found from the Monte Carlo inversion. Notice the reversed melt rate axis. NorthGRIP is indicated by the dotted line.

Using Equations (12) and (13) and Q _ice = 70 mWm^–2, the geothermal heat flux at NorthGRIP is determined to be 129±2mWm^–2. Both upstream and downstream from the drill site, significantly higher values of the geothermal heat flux are found.

Discussion

The above stated uncertainties are the standard deviations of the histograms of accepted model values. They only reflect the precision with which the Monte Carlo inversion is able to determine the value of the parameters and do not include uncertainties arising from model deficiencies and assumptions. The total uncertainties of the parameters are therefore believed to be larger than the stated standard deviations.

The ratio between the accumulation rate at NorthGRIP and at other locations along the flow line was assumed constant in time. At present the ratio of the accumulation at NorthGRIP to that at GRIP is 83%, but Reference Grinsted and Dahl-JensenGrinsted and Dahl-Jensen (2002) found that this ratio was as low as 66% during the glacial period. This indicates that the accumulation ratio at other places along the line may also have changed in time. Consequently, the assumption of unchanged accumulation pattern along the line with time may be poor. The results from Reference Grinsted and Dahl-JensenGrinsted and Dahl-Jensen (2002) indicate that the accumulation pattern seen today in the area between GRIP and NorthGRIP was more pronounced during the glacial period. If this is the case, the accumulation rates used upstream from NorthGRIP in this model are slightly overestimated for the glacial period, resulting in an underestimation of the melt rates.

The fraction of basal sliding was assumed to be linearly related to the melt rate (Equation (5)). This is based on the premise that a higher melt rate will provide a larger amount of water to lubricate the bed and thus result in a larger sliding velocity. However, this assumption may not hold if the meltwater is drained from the area where it is produced e.g. through valleys or channels. Thus, in assuming Equation (5) is correct, we also assume that the meltwater does not move far from where it is produced.

This study aims to estimate the basal melt rate at NorthGRIP, yet we use a non-thermal model. This can be done because the basal melt rate equals minus the vertical velocity at the base of the ice sheet and thus can be treated as a flow law parameter. However, the melt rate depends on the temperature gradient at the base, which changes with time because the surface climate and therefore the temperature of the ice changes with time. Thus the melt rates found in this study may be considered as average values for the past 79.6 kyr.

Ice core studies have found values of 7mma^–1 and 140 mWm^–2 for the basal melt rate and geothermal heat flux at NorthGRIP (North Greenland Ice Core Project members, 2004). These values fall in between the values found from the two- and one-dimensional models, and considering the assumptions made in the model, the results found in this study do not disagree with those obtained from ice core studies.

Conclusions

The basal melt rate at NorthGRIP is found to be 8.2 mma^–1 using the one-dimensional model and 6.1 mma^–1 using the two-dimensional model. The difference between the two numbers illustrates the importance of using a two-dimensional model even though the computational time is significantly larger.

The basal melt rate is found to vary between 5.3 mma^–1 and 21.2 mma^–1 along the flowline. Assuming the variation is caused by geothermal heat flux variations, Q _geo varies between 121 mWm^–2 and 231 mWm^–2 over scales of 10 km. This requires the sources for the changes in geothermal heat flux to be located near the surface. Large spatial variations in the geothermal heat flux have also been reported by Reference Näslund, Jansson, Fastook, Johnson and Anders-sonNäslund and others (2005). From studies of the Fennoscandian ice sheet during the Last Glacial Maximum they found significant local changes in the geothermal heat flux in Sweden and Finland. The values of the geothermal heat flux found in the present study are, however, quite high.

The drainage system of the meltwater created under the Greenland Ice Sheet is not well known. The water may be transported through small valleys observed in the bedrock topography. The presence of such canals may cause rapid spatial variations in the melt rate and is an alternative way of producing high local melt rates without strong changes in the geothermal heat flux. This is supported by the fact that dips in the isochrones are often observed over the small valleys in the bedrock.