JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 108, NO. B3, 2143, doi:10.1029/2001JB001731, 2003
[10] The 3-D ice sheet model (ISM) that we have developed and employ herein has full thermomechanical coupling (including bed thermodynamics). It is coupled asynchronously to a physically based viscoelastic model of the glacial isostatic adjustment process to describe the time-dependent displacement of the surface of the solid earth due to surface loading and to a positive degree-day surface mass balance model with temperature-dependent degree-day coefficients and physically based refreezing model. Model parameters that appear in the discussion to follow are summarized in Table 1.
2.1. Thermomechanical Ice Sheet and Bedrock Response Models
[11] The base thermomechanically coupled model is that originally described by Tarasov and Peltier [1999] with more recent improvements and Greenland specific details provided by Tarasov and Peltier [2002]. Only a brief review of these primary components will be provided here. The ice dynamics component of the model is based upon the vertically integrated form of the equation for the conservation of mass, as
in which H is local ice thickness and G is the net surface and basal mass balance. The standard Glen flow law for ice rheology is employed to compute the horizontal ice velocity V(z) with a factor 5.1 flow enhancement (unless otherwise stated). The temperature dependence of the ice rheology is as per the European Ice Sheet Modeling Initiative (EISMINT) II intercomparison project specifications (available at http://gopher.ulb.ac.be/phuybrec/eismint.html; also see C. Ritz et al., manuscript in preparation, 2002). A factor 2.5 reduction of the flow parameter is applied to interglacial ice (i.e., ice that formed during either the Holocene or Eemian periods) in accord with observations [Dahl-Jensen and Gundestrup, 1987].
[12] The computation of the ice temperature field (T) is based on the conservation of internal energy and takes into account advection, vertical diffusion, and heat generated by deformation heating (Qd) as represented by the following partial differential equation:
Heating due to sliding is accounted for in the basal thermal boundary conditions as discussed by Tarasov and Peltier [1999]. The thermodynamic solver is fully coupled to the ice dynamics module, and in the base configuration includes 65 levels in the vertical scaled to the local ice thickness. The model for the internal temperature field is also implicitly coupled to a 5 level thermodynamic (vertical diffusion of heat only) bedrock model that spans a depth of 2 km. The deep geothermal heat flux employed herein as the lower boundary condition for the thermodynamic bedrock model is described in a subsequent section.
[13] Bedrock response is computed on the basis of a complete linear viscoelastic field theory for a spherically symmetric Maxwell model of the Earth [Peltier, 1974, 1976]. For an arbitrary surface load per unit area L(
, 
, t), the bedrock displacement R(, 
, t) is governed by the following space-time convolution:
in which 
(
, t - t
) is the radial displacement Green function [Peltier, 1974]. The radial response is evaluated using a spectral representation [Peltier, 1976] of the convolution integral (3) truncated at degree and order 256.The radial viscosity profile is represented by that of the VM2 model [Peltier, 1996; Peltier and Jiang, 1996] with a 90 km thick lithosphere and the PREM model [Dziewonski and Anderson, 1981] is assumed to describe the elastic structure. Bedrock response is computed every 100 years and is asynchronously coupled to the ice sheet model. The model is initialized with the present-day observed surface and bedrock topography. It is first brought to thermomechanical equilibrium (under the assumption of isostatic equilibrium) using the climate forcing inferred for 250 ka. The model is then run with full dynamical coupling from 250 ka to present using a climate forcing based primarily upon the summit 
18O record as discussed below. A more complete description of the isostatic adjustment model is provided by Tarasov and Peltier [2002] while the relative sea level solver has been reviewed by Peltier [1998].
[14] The surface mass balance model is based upon a positive degree-day representation of the ablation process with temperature-dependent positive degree-day coefficients derived from the analyses of Braithwaite [1995]. The influence of surface refreezing is included, taking into account both capillary retention and latent heating following the methodology of Janssens and Huybrechts [2000].
|
|
[15] As has become standard for ice sheet models that are not explicitly coupled to General Circulation Models of climate evolution, precipitation is computed by assuming a surface temperature-dependent perturbation to the present-day observed precipitation climatology P(0, x, y), as follows:
with a base value of 0.062 for the parameter 
, chosen to provide a good fit to both the inferred age profile of the GRIP ice core (see Figure 1) and to the observed borehole temperature profile.
[16] We have also found it necessary to explicitly allow for some regional sensitivity of the local precipitation rate upon temperature (represented by the parameter xy(x, y) in equation (4)) in order to better fit relative sea level observations from the coast of Greenland as well as additional topographic constraints. Briefly, xy(x, y) is linearly reduced to 0.5 from 68°N to 60°N (then held at 2.0 south of 60°N), is increased to 2.0 in the northeastern sector of the ice sheet and is elsewhere equal to 1.0. The southern zone reduction in xy(x, y) offsets the decrease in glacial precipitation due to the stronger climate forcing imposed in that region as described in section 2.3.
|
|
[17] A slope-dependent factor is also included to capture orographic impacts on precipitation and follows the form employed by Ritz et al. [1997] though with an added weakening of the slope-dependence with elevation so as to better match observed borehole temperature and age profiles and inferred accumulation rate history at GRIP. Specifically, is defined as
where h is the surface elevation. The resulting accumulation rate history that is predicted by the model for the GRIP site, shown in Figure 2, is at almost all ages less than that inferred for GRIP by Johnsen et al. [1995] but generally lies above that inferred for the nearby GISP II site by Cuffey and Clow [1997]. The age profile for our base model (denoted GrB) at both the GRIP and GISP II sites, as shown in Figure 1, matches the corresponding inferred age profiles to within the well-constrained GISP II age profile uncertainties discussed by Meese et al. [1997], at least for the top 2860 m (corresponding to approximately the last 150 kyr). At the GRIP site, a 60 m model misfit in ice thickness results in discrepancies for the bottom of the core.
[18] The present-day precipitation climatology utilized in the model was presented by Tarasov and Peltier [2002] and is a combination of results obtained from the inversion of a new digital accumulation map for higher elevations [Ohmura et al., 1991] and an older precipitation map for lower elevations[Ohmura and Reeh, 1991] with present-day seasonal variability taken from Legates and Willmott [1990]. Snow fraction is computed assuming a normal distribution for temperatures around monthly means with solid fraction assumed for temperatures below 2°C.
[19] Calving of ice is a key component of Greenland mass balance. The representation of calving processes is problematic in light of the subgrid-scale processes involved. Given the lack of a consensus alternative, we simply calve ice in proportion to the excess buoyancy of the marginal ice with the time-independent proportionality factor hand-tuned to provide a match between model and observed relative sea level chronologies as detailed by Tarasov and Peltier [2002].
2.3. Temperature Forcing: Determination of the Proxy for Climate Variability
[20] We use the 18O GRIP record as a proxy for regional climate variations and constrain the climatic isotopic sensitivity 
c = d18O/dT to be less than the present-day spatial sensitivity of approximately 0.67‰ °C-1. Holocene and glacial values for c are obtained so as to provide a close match between the model and observed borehole temperature profile at GRIP. Following Cuffey [2000], we will take into account the effect of the observed isotopic lapse rate in central Greenland on the temporal transfer function from 18O to temperature forcing. Specifically, using the observed value for the 18O lapse rate in central Greenland of 
= -6.2‰ km-1 [Johnsen et al., 1989], the regional climate forcing 
Tc will be assumed to be given by
in which ZG(t) is the change in surface elevation relative to present-day at the GRIP site. Surface temperature change is then computed from Tc using a constant environmental lapse rate of 7.5°C km-1. Present and past surface temperature gradients are computed with the widely used parameterization of Huybrechts et al. [1991] based on the surface temperature maps of Ohmura [1987].
[21] The largest single uncertainty in glacial cycle modeling of large-scale ice sheets is that of the accuracy of the climate forcing. Given the variability of midlatitude storm tracks and their impact on regional Greenland climate, especially during glacial periods, the use of a single proxy to drive climate changes over all of Greenland is clearly suspect. Using constraints from RSL observations [Tarasov and Peltier, 2002] and observed borehole temperature profiles, we have found it necessary to include three further modifications to the regional temperature forcing. First, we have imposed at 5% per degree latitude increase in the amplitude of the 18O derived temperature forcing (relative to present-day) for all areas south of summit. This is in approximate accord with the results of paleothermometry for the Dye 3 borehole which indicates a 50% larger amplitude in the temperature signal there [Dahl-Jensen et al., 1998], which is most probably due to the closer proximity to the North Atlantic storm track which is thought to have migrated northward during glacial periods [e.g., Fawcett et al., 1997]. Second, we have imposed a set of ad hoc regional modifications to the Holocene temperature forcing in order to obtain a close match with observed relative sea level (RSL) histories and present-day net mass balance estimates for the whole ice sheet [Tarasov and Peltier, 2002]. These further modifications are applied only to the near marginal regions(ablation zone) and therefore have no direct thermal impact on the borehole temperature profile for the summit region ice cores. Third, we have also added a slight 0.3°C warming and cooling during the last 1.2 kyr in order to obtain a better match to the observed borehole temperature profile. In detail, the temperature forcing is lowered by -0.3°C linearly from 1.2 ka to 900 years before present and then elevated linearly to the unmodified level by the present time.
2.4. Semi-Lagrangian Tracer Module
[22] Passive tracer fields (F) such as ice formation date and 18O ratios for ice parcels are conservatively advected through the ice. Their transport is therefore described by
with appropriate boundary conditions at the external surfaces of the ice sheet.
[23] Tarasov and Peltier [2002] described the implementation of a finite difference tracer module for transporting the ice parcel formation date through the ice. Application of this module to track 18O ratios was found to be problematic due to the action of excessive dissipation. Dissipation can be reduced through the use of higher-order finite difference methods and increased resolution, but only at great computational cost. Numerical representations of advective transport using finite difference methods are subject to the CFL criterion that constrains advective transport to advance no more than a single grid cell per computational time step. The semi-Lagrangian methodology avoids this limit and easily extends to higher-order form. Semi-Lagrangian methods assign transported field values based on the value of the field at the precursor point from the previous time step [e.g., Staniforth and Cote, 1991]. This is the point (x - U*t) from which the contemporaneous velocity field U* would have transported the tracer field F(t, x) to the current grid point over one tracer time step, i.e.,
[24] Application of this methodology requires two computations every tracer time step. First, an appropriate mean velocity field U* must be computed so that a precursor point can be determined. Second, the value of the field at the precursor point (which is generally not a grid point) must be computed using some form of interpolation. Previous theoretical and analytical studies have generally found that while at least a third-order interpolation scheme is required for the interpolation of the field value to the precursor point, linear interpolation is adequate for determining the midpoint velocity value that is used to represent an effective mean time step velocity [Staniforth and Cote, 1991].
[25] The midpoint mean velocity U* is obtained from the condition
and could be interpreted as implementation of the midpoint form of numerical integration. As a test of the stability of the method, we have also implemented an alternative expression for the midpoint field that corresponds to the trapezoidal form of numerical integration, namely:
Use of this form for U* resulted in no discernible difference for a sine wave test field as compared to the prior U* form. The above recursive equation for U* is iterated up to a maximum of 5 times or until the weighted point-wise difference in computed midpoint position between iterations is <3 m. The explicit mathematical form of this condition is then
For each tracer time step, the previous U* field is used to initiate the iterations.
[26] Field values for precursor points situated above the ice sheet surface are set to the contemporaneous surface value of 18O(Tsurface, h), as inverted on the basis of the temperature forcing (determined by equation (6)), and zero age (or, numerically, the ice formation date is set to the contemporaneous time). To allow finer resolution, the tracer module is implemented on a subgrid that does not cover the whole ice sheet. The resolution of the subgrid is given in Table 2. Aside from investigations incorporating the NGRIP site, the subgrid was bounded by 40.25° and 34.25° west longitude and by 71.1° and 74.1° north latitude so as to incorporate all contemporaneous summit positions that obtain during the transient run. Experiments with an expanded horizontal subgrid boundary found no discernible impact on near-GRIP site tracer field results.
[27] It should be noted that since an age field is only computed for the subgrid region, the interglacial flow parameter reduction described above is only applied in the region of the subgrid for model runs using the SL tracer module. Model runs without the SL tracer module use a coarse resolution finite difference tracer calculation to compute age fields for the whole ice sheet. Sensitivity analyses (not shown) indicate that this geographic restriction of the flow enhancement has a relatively minor impact on the modeled ice sheet, at least for the analyses presented herein.
[28] We have found tricubic spline interpolation [e.g., Cheney and Kincaid, 1985] to be the most accurate given available computational resources. It also has the useful feature of conserving mass for divergence-free flows [Bermejo, 1990]. In the tricubic spline case, the interpolation is centered around the point in question. For precursor points near the horizontal boundary of the subgrid, the tracer model reverts to lowest-order bilinear interpolation horizontally but still remains cubic spline for the vertical interpolation. In order to ensure stability, points with ice thickness less than 100 m are assigned surface boundary values throughout the depth.
|
|
[29] Three-dimensional interpolation is very intensive numerically and generally involves a Cartesian product of 1-D interpolations that in this case first involves interpolation onto the z plane of the point in question, followed by interpolation within this horizontal plane. For this reason, a fourth-order 3-D interpolation will require a minimum of 64 linear interpolations for each point. In order to preserve the phase and low- to midfrequency amplitude components of the Eemian segment of the tracer 18O chronology in the model and remain within available computational resources, we have found the best trade-off between accuracy and resources to be obtained with a 1/32° by 1/64° horizontal (longitude, latitude) resolution and 1025 equidistant layers in the vertical for the tracer subgrid. This vertical resolution corresponds to a temporal resolution at the GRIP site of about 1200 years at 2750 m which corresponds to a time approximately 100 kyr ago. Theoretical analyses indicate that amplitude and phase errors of SL tracer results are affected by both the relative wavelength (wavelength/grid size) and by the Courant number (ut/x) [McDonald, 1984]. Phase error decreases with the Courant number while amplitude errors are minimized for near integer values of the Courant number. For deep ice (corresponding to Eemian or pre-Eemian periods) in the summit region, the Courant number for the coarse resolution ISM grid is less than 0.03 for the horizontal projection. The presence of significant dissipation with tracer models using this coarse horizontal resolution makes it clear that the Courant number is not the critical factor here. Increased relative wavelength will decrease both phase and amplitude errors. Analyses suggest that a relative wavelength of order 10 or greater is generally required to avoid significant dispersion and dissipation over long-term integrations [McCalpin, 1988]. This obviously justifies the need for high vertical resolution. Given the significant improvement of the tracer signal with higher horizontal resolution that is clearly evident in Figure 3, it is evident that even the horizontal projection of the 18O signal has significant high-frequency variability. This is especially true near the summit, where the very low horizontal ice velocities and diverging flow results in short relative wavelengths.
[30] It is generally observed that dissipation decreases with increasing time step [McCalpin, 1988] due to the decrease in the total number of interpolations involved. In fact, if the precursor point could be determined exactly, then dissipation and dispersion error would be minimized by using a single time step [McDonald, 1984]. There is a trade-off with shorter time steps needed to reduce uncertainty in the precursor point and to maintain temporal resolution of the tracer field near the boundaries. As such, we vary the tracer time step from 500 to 20000 years, with shorter time steps reserved for periods of interest (e.g., the Eemian). Tests with different tracer time steps had no impact on tracer signal phases. Tracer amplitudes for, say, the Eemian segments of model cores were generally improved using large post-Eemian time steps for the tracer module.
2.5. Intercomparison of Different Semi-Lagrangian Tracer Methodologies
|
|
[31] An extensive comparison of different semi-Lagrangian (SL) models was undertaken using the ISM horizontal resolution of 0.5 by 0.25 degrees longitude/latitude. Both linear and quadratic interpolations for the midpoint velocity determination were compared, as were linear, cubic spline, and order 3 to 7 Newton-Lagrange interpolation methods (a comparison of a number of these methodologies for a 10 kyr sine wave test forcing is shown in Figure 4). This has led us to the following conclusions that may be valuable to others involved in similar work.
[32] First, even finite difference methods work reasonably well for smooth monotonic fields such as ice age, though discernible improvements are obtained with higher-order SL tracers.
[33] Second, for fields with strong high-frequency variance, such as 18O, tricubic spline interpolation combined with high horizontal and vertical resolution are required to preserve signal phase and amplitude into the Eemian period as demonstrated in Figure 3. For more recent periods, horizontal resolution can be degraded without significant impact. Trilinear interpolation is not appropriate for signals with high-frequency variance.
[34] Third, in accord with experience in the numerical weather prediction field, higher-order interpolation for the midpoint evaluation did not produce any significant improvements in the quality of the tracer signal recovered.
[35] Fourth, tricubic spline interpolation for the determination of the tracer field value at the precursor point preserves tracer signal quality better than even eighth-order Newton-Lagrange interpolation with quasi-monotonic constraints. Other studies [Reames and Zapotocny, 1999] find no significant improvement in tracer field conservation with Newton-Lagrange interpolation using orders higher than 8.
[36] We have also explored the impact of employing quasi-monotonic constraints on the precursor field value interpolation. This constraint enforces the bounding of interpolation results by the extremal values at the grid box vertices enclosing the point in question [Bermejo and Stanforth, 1992; Williamson and Rasch, 1989]. Newton-Lagrange interpolation requires such a constraint to avoid significant spurious overshoots of the signal. Spline interpolations can also overshoot, but for the 18O signal the overshoots were relatively insignificant. Furthermore, for the sine wave test field, the imposition of quasi-monotonic constraints does slightly increase signal dispersion in accord with previous studies [Bermejo and Stanforth, 1992].
[37] In summary, the high-resolution tricubic spline method appears to be the most computationally efficient SL tracer methodology to employ with ice sheet models for signals with high-frequency variance. We are able to preserve phase and a significant fraction of amplitude at the GRIP site reaching back to the Eemian interglacial. For the top section of the model core corresponding to the last 60 kyr, we have near complete capture of the 18O signal at the computational temporal tracer resolution.
[38] As detailed by Tarasov and Peltier [2002], our model version GrB is tuned to five primary constraints: surface elevation at summit (GRIP), borehole age, and temperature profiles at GRIP, a set of 14C dated RSL records, and approximate basal temperatures at Camp Century, Dye 3, and GRIP. Secondary constraints that were not necessarily fully compatible with the primary constraints include GRIP accumulation chronology between that inferred for GRIP and GISPII, minimized RMS difference with the observed surface topography, minimized discrepancy with the observed ice volume, and closeness of fit to observed ice thickness at GRIP, Camp Century, and Dye 3. Tuning parameters for these constraints are summarized in Table 3.
|
|
[39] With respect to observations, discrepancies in ice thickness (secondary constraint) range from about 65m at GRIP to about 250 m at Dye 3. It is clear from Table 4 that many of these differences are attributable to errors and resolution limitations associated with the input topography, though other factors such as the simplicity of the climate forcing used to drive the model, the assumption of contemporaneous isostatic equilibrium for model initialization at 250 ka, and the lack of explicit ice stream dynamics must also play a role. The overall RMS topographic difference of the ice sheet with respect to the digitized observational topography that was used to initialize the model is a respectable 70 m. Model GrB also has the advantage of a perfect (to grid resolution) match to present-day summit location which provides some minimum confidence in our subsequent analyses of summit migration. Figure 5 shows a sequence of snapshots of the predicted evolutionary history of Greenland, in terms of topography, from the baseline GrB reconstruction developed by Tarasov and Peltier [2002]. It is notable that for this model the Eemian configuration of the ice sheet at 121 ka is not significantly different from the present configuration of the ice sheet.
Citation: Greenland glacial history, borehole constraints, and Eemian extent, J. Geophys. Res., 108(B3), 2143, doi:10.1029/2001JB001731, 2003.
