Introduction

The stratigraphic records in the deep ice cores recovered from central Greenland have provided a continuous detailed record of past variations in the environment around the borehole sites. The stratigraphic record in turn is also altered by dynamical changes in the ice sheet itself, and thus contains an integrated history of the ice-sheet evolution. Here, we use an inverse approach to infer the ice-thickness history, h(t), during the Holocene along the flowline connecting the GISP2 and GRIP drill sites, based on the measured depth–age relation in the GISP2 core. We also obtain error bounds on the derived h(t), which result from the propagation of the data uncertainty into our forward model.

Our derivation of h(t) is possible primarily for two reasons. First, annual layers can be resolved in the GISP2 core well beyond 20 000BP, which enables an accurate depth–age relation to be established (Reference AlleyAlley and others, 1993). Secondly, from mass continuity the ice velocity is a function of the change in ice thickness with time h, the accumulation rate a, and the ice-sheet geometry, and also depends on the ice rheology. Thus the time needed for an ice particle to traverse the trajectory from the surface to a particular depth in the core is a function of h, as well as of the accumulation rate, geometry and rheology. If these variables can be specified as a function of time, we can compute the time needed for an ice particle to travel from the surface to a given depth in the core. With this forward model to compute ages at various depths in the core, we can then use inverse methods to iteratively adjust an initially specified h(t) by comparing computed and measured ages at a number of depths until satisfactory agreement is reached. The resulting h(t) is the derived ice-thickness history.

It is important to note that we make no assumptions concerning the dependence of the ice thickness on the accumulation rate, and, more generally, make minimal assumptions about the physical processes that cause the ice sheet to evolve. Our derived h(t) is simply that function which results in the best agreement between measured ages and those computed by our forward particle-trajectory model.

However, we do make two important assumptions, which we discuss in the next section. One is that the surface topography in the region remains unchanged over the time interval of interest. We make this assumption because while the position of the flow divide and the flow divergence along the flowline can have important effects on the calculated ages, the variation of the divide position and flow divergence with time is not known. The other assumption is that the depth variation of the horizontal velocity (the shape function) as a function of position along the flowline is also time-independent. This results from our inability to compute a solution for the velocity field with a time-varying rheology. Later we will test the sensitivity of our results to both of these assumptions.

Whether or not the Greenland ice sheet actually evolved according to these assumptions, the class of ice-sheet evolution models that includes these assumptions should give results that fall within the error bounds we derive, This is the strength of our approach. The model of Reference Cutler, Raymond, Waddington, Meese and AlleyCutler and others (1995), who derive a thickness history from the layer data by using a specified ice-flow model to relate flow changes to changes in ice-sheet thickness and width, is a member of this class.

Method

We compare measured ages at various depths in the GISP2 core with ages computed by a forward model that calculates particle trajectories. The depths correspond to measured ages at 1000 a intervals from 1000 to 20 000BP. The calculated ages are functions of parameters that describe the change in ice thickness with time, and we take these parameters to be values of h, at 2000 a intervals from 0 to 20 000 BP. This results in a seemingly over-determined inverse problem with 20 data points and 11 parameters.

Changes in the surface topography with time can have a significant effect on calculated ages. For example, as the divide migrates away from the GISP2 site, the length of the flowline increases, which increases the horizontal velocity at a given point. While this effect tends to reduce the travel time, it is compensated for by the increase in the length of the trajectory. Also the rearrangement of the surface elevation contours associated with the divide migration would affect the position of the flowline and the flow divergence, but whether this rearrangement results in the horizontal velocity increasing or decreasing depends on the unknown details of how the surface varies with time.

Some support for the assumption that the current surface topography has been unchanged over the past 20 ka is provided by Reference Anandakrishnan, Alley and WaddingtonAnandakrishnan and others (1994). Using a steady-state, one-dimensional flowline model, they found that the divide position in central Greenland is relatively insensitive to ice-sheet width, with the maximum displacement being about 10 km as the ice margin migrates seaward by 250 km. They also found, as did Reference Letréguilly, Reeh and Huybrechts.Letréguilly and others (1991), that the ice-thickness change in the summit region was probably less than ± 10% over the past 20 ka. These results suggest that changes in surface geometry were likely to have been small, and our assumption of a steady ice-sheet surface is not unreasonable. Later we will examine the effects of changing the flow divergence on the derived ice-thickness changes.

We make use of the shape-function approach in calculating the ice velocity as a function of time. The essence of this method is that the horizontal velocity is written as the product of two terms. One is the mean or depth-averaged velocity, which is computed from mass continuity and is a function of horizontal position and time; the other is a shape function which describes flow the horizontal velocity varies with depth. In principle, the shape function may also be time-dependent, although the greatest utility of the method is when the shape function can be taken to be time-invariant. Here, we infer the shape function using the two-dimensional finite-element ice-flow model described in Schott and others (1992), which computes a steady-state velocity field along the flowline connecting the two boreholes.

While the use of a shape function derived from Holocene conditions may be appropriate for the past 10 ka, the shape function prior to the Wisconsin–Holocene warming may have been significantly different due to colder temperatures, particularly in the upper part of the ice sheet, which make the ice stiffer through the Arrhenius term in the flow law. Also, during the late Wisconsin the upper zone was impurity-laden, which tends to soften the ice (Reference PatersonPaterson, 1991). However, whether the temperature and impurity effects compensate each oilier is unclear. As a result, here we determine shape functions for late-Wisconsin and Holocene conditions and derive ice-thickness changes for each shape function, rather than attempt to interpolate between the two shape functions during the Wisconsin–Holocene transition. We find no appreciable difference in the resulting ice-thickness changes.

The Particle-Trajectory Model

We compute the age of a layer at a specified depth in the core by determining the path traveled by the layer since its deposition. The computation of particle trajectories is based on the equation of continuity and the assumptions discussed above. We also assume dial the flow is planar. The shape-function approach we use here is similar to that of Reference Kostecka and WhillansKostecka and Whillans (1988) and Reference ReehReeh (1988), except that we take the shape function to be known and allow it to vary with horizontal position. The horizontal velocity components u_x and u_y are given by

where z is positive upward, ϕ is the shape function, ū, is the depth-averaged velocity component, and i = x,y. The depth-averaged components are obtained from the equation of mass continuity,

where h(x,y) is the ice thickness, a(x,y,t) is the accumulation rate in m ice a⁻¹. and the time t is positive toward the future.

Following Reference Kostecka and WhillansKostecka and Whillans (1988), we simplify the continuity equation by defining the coordinate system so that the x axis follows the flowline, with x = 0 at the GRIP site, and write the y component of mean velocity as:

If flow is perpendicular to the elevation contours, the spreading parameter R is the radius of curvature of the elevation contours (Reference Whillans and CassidyWhillans and Cassidy, 1983). The continuity equation then becomes

which is solved with the boundary condition of zero velocity at the ice divide. The solution provides the mean horizontal velocity along the trajectory at any time for a given ice-sheet geometry {h(x,t), R(x,t)} and accumulation rate a(x,t). We take R(x) = R ₀ x and, using the shape function discussed below, we find that R ₀ = 1 makes the surface velocity at GISP2 about 1.3 ma⁻¹, in good agreement with the result measured by Waddington and Morse (Reference WaddingtonWaddington and others, in press). Also, R ₀ = 1 implies that the elevation contours are concentric circles, a reasonable approximation to the current surface topography between the two drill sites (Reference Hodge, Wright, Bradley, Jacobel, Skou and Vaughn.Hodge and others, 1990). If we take a and h constant along the flowline, Equation 4 can be solved in closed form by use of an integrating factor:

As discussed above, we determine the shape function from a finite-element simulation of the velocity field between the boreholes for both late-Wisconsin and Holocene conditions. For the Holocene simulation the surface temperature is −30°C, and the softness parameter in the flow law is enhanced by a factor of three below the Wisconsin–Holocene transition at about 1700 m depth. For the late-Wisconsin ice sheet the surface temperature is −40°C, and the enhancement factor is three throughout the entire ice thickness. Here the basal topography between the drill sites is based on the radio-echo data reported in Reference Hempel and Thyssen.Hempel and Thyssen (1992), while the surface topography agrees to within 5 (10) m with the airborne radio-echo survey of Reference Hodge, Wright, Bradley, Jacobel, Skou and Vaughn.Hodge and others (1990) for the Holocene (Wisconsin) simulation. The finite-element grid used here is shown in Figure 1.

Fig. 1. A depiction of the finite-element grid used to generate the shape function, with the nodes shown by the dots. The basal topography between the GISP2 and GRIP drill sites is taken from Reference Hempel and Thyssen.Hempel and Thyssen (1992), while the basal topography on either side of this zone is based on Reference Hodge, Wright, Bradley, Jacobel, Skou and Vaughn.Hodge and others (1990). Also shown are the particle trajectories corresponding to ages of 5, 10 and 20 ka.

The ice-flow model generates velocity components at the nodes of the finite-element quadrilaterals. We convert the discrete horizontal-velocity values to shape-function values by dividing u at fixed x by the corresponding surface value. We then expand the discrete set of shape-function values in Chebyshev polynomials to obtain an interpolating polynomial representation for ϕ(x,z), which can be readily integrated and differentiated.

The vertical velocity is computed by integrating the vertical strain rate, which is obtained from Equation (1), Equation (3) and Equation (5), and the condition that ice be incompressible. The result is

and integrating Equation 6 from the base to height z gives the vertical velocity w(x,z):

where we have taken the ice sheet to be frozen to the bed (Reference Firestone, Waddington and Cunningham.Firestone and others, 1990) and have ignored the effects of isostatic adjustment.

Computation of the Particle Trajectory

Let x ₀ be the distance of the borehole from the flow divide, and from now on let u = u_x denote the horizontal velocity component. To generate a particle trajectory, a starting-point (x ₀,z ₀) corresponding to depth ζ in the borehole is given, with t ₀ = 0 the initial time. The calculation proceeds backward in time. The components u ₀ = u(x ₀,z ₀,t ₀) are computed, and a first estimate for the change in vertical position resulting from a horizontal step Δx (negative in our coordinate system) is

The next point on the trajectory would be

with the corresponding time t ₁ at (x ₁,z ₁) computed from

The vertical position and time are further refined by iteratively averaging the velocity components at the beginning and end of the horizontal step, so that after the ith iteration:

where and similarly for . The iteration continues until the times computed in two successive iterations converge satisfactorily. Then the program moves to the next horizontal position up the flowline, x ₂ = x ₁ + Δx, and the process is repeated. This backward-stepping calculation is continued until the trajectory comes within 0.1 m of the surface. The trajectory origin, time and velocity are found by linearly extrapolating the appropriate values at the two adjacent points just below the ice-sheet surface, The corresponding time is then the age of the layer at depth ζ in the core.

The size of the distance step has to be chosen so as to maximize accuracy within a feasible computational time. After some numerical experimentation, we choose Δx for each step such that the estimated time to take the step is less than some value dT, which varies linearly from 30 years for the 1 ka trajectory to 150 years for the 20 ka trajectory. Between 11 and 15 years, when the accumulation rate varies rapidly, dT is reduced to 40 years.

Finally, note that the change in position over one step depends only on the velocities at the beginning and end of the step. If the velocity at (x ₂,z ₂) is computed from the accumulation rate at t ₂ alone, small changes in convergence parameters or the length of the distance step can result in large changes in travel time when the accumulation rate is highly variable with time. This is because the time t ₂ computed at each iteration is used to compute the velocity for the next iteration, and the velocity depends on the accumulation rate at t ₂. Thus if t ₂ differs by even a few years from the previous iteration, the velocities computed from one iteration to the next can vary wildly if the accumulation rate varies greatly between the t ₂ computed for the two iterations. To remedy this, the accumulation rate used to compute (u ₂,w ₂) is not the value at t ₂, but the average value over the interval t ₂ – t ₁.

Accumulation Rate

To generate the accumulation-rate time series, we determine how measured layer thicknesses are related to the original thicknesses at the surface. If λ(ζ) is the layer thickness at depth ζ, and λ₀ is the original thickness, then the two are related by:

where the summation is over the number of distance steps needed to compute the trajectory, and the subscripts on and u refer to the trajectory coordinates at the beginning and end of each step. We refer to D(ζ) as the dilation factor, which is the ratio of the original surface layer thickness to the thickness measured in the core. Note from Equation (5) and Equation (7) that since a and ḣ are independent of x, the ratio is independent of a – ḣ. Thus the dilation factor is independent of the accumulation rate history a(t), so that the reconstructed a(t) depends only on the ice-sheet geometry, the shape function and the layer data. We can assume a simple accumulation-rate history to compute the dilation factor and generate a(t), as was done by Reference Schøtt, Waddington and RaymondSchøtt and others (1992) and Reference AlleyAlley and others (1993). The a(t) computed by our model is essentially identical to that of Reference AlleyAlley and others (1993) because (1) we use nearly the same velocity held in our calculations, (2) the accumulation rate varies by less than 10% along the GISP2–GRIP flowline (Reference Bolzan and Strobel.Bolzan and Strobel, 1994), and (3) we assume that ḣ is independent of position along the flowline.

Maximum-Likelihood Inverse Method

We use the maximum-likelihood inverse method to infer the set of model parameters from the measured ages in the core. In what follows we briefly discuss some of the important features of this method; more complete discussions can be found in Reference Tarantola and Valette.Tarantola and Valette (1982) and Reference MenkeMenke (1989).

We pose our problem by using 20 data values (t_i ) to determine 11 model parameters (h_j ), so that the problem appears solvable by non-linear least-squares methods, which derive the model parameters that minimize the (weighted) residuals between the measured and modeled data values. However, the solution depends critically on the a priori information incorporated into the problem. For example, assume we know nothing about the model parameters, so that our a priori information is total ignorance. Then the least-squares method is free to generate model parameters that give very small residuals but may result in derived physical variables (such as h(t) in our case) which oscillate wildly and have little physical significance. This happens because the data are not exact but contaminated with random errors, so that the lack of smoothness in the solution reflects the least-squares attempt to fit the random error in the data.

Often we do have some a priori information about the model parameters. As discussed above, we have reason to believe that the change in the ice-sheet thickness over the past 20 ka was probably < 200 m. One way to incorporate this information into the inversion process is by assuming an a priori Gaussian joint probability distribution for the model parameters, characterized by a mean and standard deviation. Thus if each model parameter is believed to lie within some range, the mean and standard deviation could be taken as the center of the range and one-half the range, respectively, while for the state of total ignorance the standard deviation would be infinite.

Now consider the N + M-dimensional space spanned by the N parameters and M data, and suppose the data d and model parameters m are related by d = g(m). Given the a priori probability distribution for the model parameters and the probability distribution that describes the data, the probability of observing a given point in this space is the product of the data and parameter probability distributions. The non-linear model which relates the parameters and data forms a hyperplane in this space. The maximum-likelihood method then looks for the point on the hyperplane for which the joint probability distribution is maximized.

II the computed data values are non-linear functions of the parameters, then the search for the maximum-likelihood point is done iteratively. After the kth iteration the new parameters are (Reference Tarantola and Valette.Tarantola and Valette, 1982):

where Λ_d , and Λ_m are the data and a priori model-parameter covariance matrices, respectively; m⁰ are the a priori parameter values; Ḡ is the matrix transpose; and G^k has the matrix elements

Here, we compute the matrix elements by central differencing:

where g_i (m) is the computed value of the ith datum, and after some numerical experimentation we take δ = 10⁻¹⁰ ma⁻¹.

In practice, the iterative process is stopped when the solution for successive iterations changes by a negligible amount or when x ² is minimized, where for the kth iteration:

with d_i , the ith datum and σ_i , the corresponding standard deviation.

Recall that here the parameters are the rate of change of thickness with time at 2000 a intervals (h_j, j = 1...11). We start the inversion by choosing an initial point in parameter space. If our parameterization of the problem is adequate, the inversion results must be independent of the starting-point. To check this we compute thickness changes with the initial ḣ_j = 0.00,0.01 and −0.01 ma⁻¹. We choose the a priori parameter standard deviation, σ_p, to correspond to an uncertainty in the thickness of ±200m over the past 20 ka; this gives σ_p = 0.01 ma⁻¹. For the data, the estimated error in the measured ages varies from ±70 a at 3300BP to ±520 a at 17380BP, to ±3000 a at 40 500BP (Reference AlleyAlley and others, 1993). We take the standard deviation in the measured ages to vary linearly with age between these values.

Results

The change in thickness with time using the Holocene rheology is shown in Figure 2 for the three different sets of a priori parameter values. All three results are in good agreement over the past 10 ka, but prior to that the computed values for ḣ are nearly equal to the a priori values. This suggests that these parameters have not been determined by our method. This is confirmed by noting that for the ḣ prior to about 10 000BP, the estimated a posteriori parameter standard deviations are nearly identical to the assumed a priori values, so that little information is being propagated from the data to these parameters. Thus, the ice-sheet history prior to about 10 000BP remains unresolved by our method.

Fig. 2. The derived ice-thickness history for assumed initial parameter values of 0.00 ma⁻¹ (solid line), +0.01 ma⁻¹ (dashed line) and –0.01 ma⁻¹ (dot-dash line). The assumed a priori parameter standard deviation in all cases is 0.01 ma⁻¹. Also shown are the ± 1 σ bounds (light solid line) on the thickness history obtained from the initial parameter values of 0.00 ma⁻¹.

Over the most recent 10 ka, our results show that the ice thickness has varied by less than ± 10 m about the current value, with an estimated a posteriori uncertainty in the thickness of about ±85m at 10 000BP. The computed current rate of thickening varies from about 4 ± 9 mm ⁻¹ to −3±9 mm a⁻¹ for a priori parameter values of +0.01 and −0.01 ma⁻¹, respectively. These results are consistent with the regional mass balance of 17±17 mm a⁻¹ measured by Reference BolzanBolzan (1992), which was based on measured surface velocities and accumulation rates over a 150km× 150 km grid centered on the Summit site.

We assumed here that the ice rheology and surface topography were time-invariant. To test the sensitivity of our results on the constant-rheology assumption, we have computed thickness changes by substituting the shape function inferred from the finite-element model of the late-Wisconsinan flowline described above. We find that h(t) differs by less than 0.6 m over the past 20 years from the result shown in Figure 2, This is because most of the trajectories are within the upper half of the ice sheet, where differences in the glacial and Holocene shape functions are small.

A varying surface topography may affect the divergence along the flowline (and also the direction of the flowline) as discussed in a previous section. We estimate the effects of a varying flow divergence by computing the thickness change with R ₀ = 10. This causes the horizontal velocity to nearly double compared to the model with R ₀ = 1 (see Equation 5), but again we find that h(t) is within 0.6 m of the result in Figure 2. This is because even though we expect computed ages to be reduced due to the increase in the horizontal velocity, the shapes of particle trajectories also change so that the trajectory lengths are increased. For example, the 20 000BP trajectory intersects the surface at x = 14.4km for R ₀ = 1, but at 8.3 km for R ₀ = 10. The result is that the increase in velocity is offset by the increase in distance traveled, and computed ages are nearly identical to those for R ₀ = 1.

We could reduce the uncertainty in the results and improve the parameter resolution in a number of ways. One would be to increase the number of data used in the inversion. This could be readily done with only a small increase in the computational cost for the more recent part of the record, but it would be much more expensive to more tightly constrain the late-Glacial ice-sheet history. We could also re-parameterize h(t) in the forward model prior to 10 000BP by specifying ḣ and ḧ at only a few points in time. This would have the effect of making more computed data values in this time interval sensitive to these parameters. Finally, a seemingly obvious way to reduce the a posteriori parameter uncertainly is to reduce σ_p, the a priori parameter uncertainty, but this would be unjustified as it implies more precise knowledge of the ice-sheet history than we actually possess.

Conclusions

We have investigated the use of a formal inverse approach to derive the ice-sheet thickness changes and the associated uncertainties over the Holocene from the measured age–depth data in the GISP2 ice core. This method strives to generate robust results, i.e., we make only a few very general assumptions and in return expect to obtain results that are applicable to all members of the much wider class of ice-dynamics models that incorporates our assumptions as a sub-set of its own. This generality is possible because our method estimates the uncertainty in the ice-sheet thickness history directly from the data, and is relatively independent of the details of the ice-sheet evolution model employed.

Our results suggest that the ice sheet has been in steady state over the past 10 ka. However, the estimated uncertainty in the ice thickness at 10 000BP is ±85m, so that at the onset of the Holocene we can say with confidence only that the thickness was within about 100 m of the current value. We find the current rate of thickening to be within ±0.01 ma⁻¹. This result is consistent with that of Reference Cutler, Raymond, Waddington, Meese and AlleyCutler and others (1995), obtained using an independent technique, and with a regional mass-balance determination by Reference BolzanBolzan (1992), based on measured surface velocities and accumulation rates over a 150 km × 150 km grid.

Prior to about 10 000BP, h(t) remains undetermined by our approach. We believe this is largely due to the parameterization of h(t) we adopted which, in retrospect, sought more resolution of the thickness history than was justified. A more robust result could be obtained by incorporating more Holocene data values and reducing the number of parameters that describe h(t) during the last glacial.

Adding more information, in the form of using more data values in the inversion or reducing the a priori parameter standard deviation by incorporating results of other determinations of h(t), can only reduce the uncertainty that we compute here. As a result, we believe the results derived here provide an upper bound to the likely uncertainty in the thickness history over the past 10 ka.