JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 108, NO. B3, 2143, doi:10.1029/2001JB001731, 2003
3.1. Borehole Temperature Profiles and Geothermal Heat Flux
 The deep geothermal heat flux from the Earth into the base of the Greenland ice sheet is largely unknown. The most recent and thorough attempt to develop a global heat flux map is based on a mixture of direct observations and bedrock geology [Pollack et al., 1993]. According to this model, which is illustrated in Figure 6 for the region of Greenland, a strong horizontal gradient should exist across this region. However, the lack of local heat flow measurements in this region leaves the model poorly constrained.
 Determination of background geothermal heat flux through the inversion of observed borehole basal temperatures using thermomechanically coupled ice sheet (and bedrock) models offers the possibility of further constraining the geothermal heat flux for ice-covered regions. Direct input of the heat flux map of Pollack et al.  into our coupled ice sheet model was found to produce excessively warm basal temperatures. Given the few available constraints, we chose to apply a series of bounded linear transformations to this heat flux map in order that the coupled model could be forced to approximately match the present-day observed basal temperatures at the three boreholes through the ice sheet from which such data is available. For the sake of simplicity and also to improve the match of the predictions of the model to RSL observations, we chose to impose much weaker gradients where there were no borehole constraints. The resultant modified heat flux map, show in Figure 6, is of course only reasonably constrained in the vicinity of the three data sites, denoted D, G, and C in Figure 6 (Dye 3, GRIP, and Camp Century, respectively).
 It is clear that borehole temperature data from Renland to the east, and from sites to the north and northeast and one from near the central west margin would allow us to obtain a much more highly constrained heat flux map for Greenland. It should also be mentioned that data sites near to, or in regions of, fast flow (e.g., Camp Century and Dye 3) would have more uncertainty associated with them, both due to the lack of explicit ice stream mechanics in the model and to the more poorly constrained climate forcing used to drive the model than is available from more inland sites such as GRIP.
 As an independent means of assessing the quality of the new heat flux map, it is worthwhile comparing computed and observed present-day surface heat fluxes. The only useful surface measurements available are from a series of boreholes in the Precambrian shield near the southern coastline of Greenland at Ivigtut and Ilimaussaq [Sass et al., 1972] denoted Iv and Il, respectively, in Figure 6. The respective measured values of 43 and 36 mW m-2 are uncorrected for thermal perturbations from past ice cover. For model GrB, the present-day surface heat flows at these respective sites are 13 and 40 mW m-2, a reasonably close match for Ilimaussaq but definitely not for Ivigtut. Understanding the source of such a large difference in modeled present-day surface heat flux between two proximal sites can help elucidate the source of the model-observation misfit at Ivigtut. The input deep geothermal heat flux for the two sites, 41 and 44 mW m-2, respectively, are clearly not the source of the 27 mW m-2 difference between the modeled present-day surface fluxes. As demonstrated by Tarasov and Peltier , who show ice thickness fields for 10 and 9 ka, both sites become ice-free at about the same time in the model (between 10 and 9.5 ka). Furthermore, both sites maintain surface elevations within approximately 200 m of each other, at least post-LGM. Climate forcing is therefore apparently not the source of the difference. Rather, the difference between the predicted heat flows is most probably due to the factor of 2 to 6 contrast in ice velocities between the two sites and the associated >12°C difference in basal temperature that existed prior to onset of the ice-free state (not shown). The higher ice velocity at Ilimaussaq resulted in increased heat advection from upstream basal ice and increased deformation heating. Throughout most of the glacial cycle, basal ice at Ilimaussaq was at the pressure-melting point, thereby allowing basal sliding with resultant surface velocities of up 310 m yr-1 at 10 ka. Ivigtut, on the other hand, never experiences basal sliding in the model, and basal temperatures were approximately -20°C at LGM and -11°C just before deglaciation. Near or above 0°C surface temperatures throughout the Holocene thereby resulted in a significant negative (i.e., downward) surface heat flux during the mid-Holocene.
 These differences between the model and observationally determined present-day heat flux values for the two sites suggest that there is insufficient model ice velocity around the Ivigtut site. This is corroborated by the enhanced Holocene forcing in the whole southwest region of Greenland required to obtain a reasonable match between model and observed RSL histories as discussed by Tarasov and Peltier . Though the forcing employed in this case has enhanced regional Holocene warming, it is clear that part or all of this forcing could be making up for the lack of explicit ice stream mechanics in the model. Increased ice velocities could result from an increased geothermal heat flux, however model tuning favors relatively low geothermal heat flux in the south. Though other factors likely play a role, ice thickness in the southern dome region is generally less than observed and an increase in deep geothermal heat flux would only further exacerbate this discrepancy.
 Partial corroboration of the northeast gradient in the deep geothermal heat flux follows from comparison of the observed and model predicted surface velocity field of the ice sheet. In a comparison with both direct and indirect measurements of surface velocity, Bamber et al.  find that a 3-D thermomechanical ice sheet model, using a spatially constant geothermal heat flux, underestimates velocities near the northern margin and overestimates velocities in the southern margin. Increased geothermal fluxes would tend to increase local velocities and therefore the spatial gradients in the geothermal heat flux map presented herein should allow a better model match to observed velocity fields. In fact, model GrB matches to within one grid point, eight of the 10 direct (i.e., GPS based) surface velocity measurements listed by Bamber et al. . This has been fully discussed by Tarasov and Peltier . The two discrepant sites were situated around Jakobshavns Isfjord, where the very high observed velocity gradients would be difficult to capture in a model of only moderate spatial resolution. Given that model GrB overestimates near margin ice thickness in the north-northwest and east-northeast [Tarasov and Peltier, 2002], it is quite possible that a much stronger northeast gradient is required for the areas covered by our geothermal heat flux map that lacked borehole constraints. However, the model will require better representation of fast flow processes in order to disentangle impacts of geothermal heat flux variations and ice streams.
 Basal temperature at summit is quite sensitive to the deep geothermal heat flux. In the neighborhood of 50 mW m-2, a 1 mW m-2 change in the heat flux will result in a 0.4°C change in present-day modeled basal temperature. We obtained a very close fit to the observed borehole temperature profile at GRIP shown in Figure 7 by using the new heat flux map which has a geothermal heat flux value of 57.56 mW m-2 below summit. The only discrepancies between the prediction of the tuned model and the observed temperature profile occur at the top and bottom of the core. At the top, the model profile has an insufficient and too near surface Little Ice Age cold bump arising at least in part from the lack of accounting for firn densification processes which are observed to be significant for approximately the top 100m. At the bottom, model basal temperature is 0.7°C too cold as compared to observations. This discrepancy was unavoidable in order to best fit the observed temperature profile. The geothermal heat flux required for our model is significantly above the median value of 51.3 mW m-2 required to obtain a close fit in the 1-D Monte Carlo inversion of the observed borehole temperature profile at GRIP by Dahl-Jensen et al. . Furthermore, this latter inversion also used a higher accumulation rate, which would have increased downward cold advection compared to model GrB and thus even further depressed basal temperature.
 This geothermal heat flux discrepancy between 1-D and 3-D modeling is due to a combination of at least four factors. First, the 3-D GrB model is tuned to match the present-day surface elevation at summit, but this does not allow for an exact match of present-day ice thickness. The model underpredicts the ice thickness at summit by 3029 - 2963 = 66 m. Extrapolation of the near basal temperature profile to the observed ice depth at GRIP results in a basal temperature 1.1°C too warm. Second, the vertical resolution of the model introduces some inaccuracies. A doubling of vertical resolution to 130 layers, does result in a slightly warmer base (-9.09°C) due to better resolution of both near-basal strain and the vertical temperature gradient. However, further doubling (-9.01°C) and extrapolation to an infinite number of levels gives -8.9°C. Neglecting the need for retuning of other model parameters, the above two factors account for about a 1.4°C/0.4 mW m-2 per °C or 3.5 mW m-2 reduction in the geothermal heat flux discrepancy.
 The remaining approximately 3 mW m-2 discrepancy between 1-D and 3-D models is likely due to some combination of the lack of accounting for the impact of longitudinal stress deviations and the limitations of 1-D models. The latter must ignore horizontal heat advection, assume a poorly constrained surface elevation chronology, and also must assume negligible impact from the displacement of summit over the glacial cycle. As will be discussed below, the model summit is located within the grid cell corresponding to the present-day summit only during parts of interglacials and never during glacial periods.
 By adjusting glacial, Holocene, and late Holocene values of c (see Table 1), we have obtained a close match to the observed borehole temperature profile at GRIP, with the only discernible discrepancy near the surface, possibly due to the lack of accounting of low-density firn layers. Also shown in Figure 7 is the resultant present-day borehole temperature profile for the untuned model GrB0 (only tuned to maximum surface elevation) and for models with modified cLH (late Holocene) and c(Holocene). It should be clear that the close fit in borehole temperature profile achieved is a nontrivial constraint on the model.
 We have also retuned the model to use the best fit climatic isotopic parameters of Cuffey and Clow  derived by 1-D model-based inversion of the GISPII borehole temperature profile. Following their model tuning, we also added a dynamic precipitation-scale correction (applied everywhere) to ensure agreement between model and inferred accumulation history at the GISPII site. With a 1.5 mW m-2 reduction in the model geothermal heat flux (which gives a value of 56 mW m-2 at GRIP) and an increase in the flow enhancement parameter from 5.1 to 6.0, this new model (“GrC”) also delivers a reasonably close fit to the observed borehole temperature profile at GRIP as shown in Figure 7. Given the significant difference in Holocene climatic isotopic parameters between model GrB and that of Cuffey and Clow  (0.364 and 0.25, respectively), it is therefore clear that the GRIP borehole temperature profile alone provides only a limited constraint on inferred climate.
 As an independent test, it is worth examining other borehole temperature profiles (which were not employed to tune the model). Model GrB somewhat misfits the borehole temperature profile at the GISP II site as shown in Figure 8. Whether this is a consequence of limited model resolution, lack of accounting for the impact of longitudinal stresses on ice flow, inappropriate ice rheology, insufficient surface gradient in climate, or significant discrepancies in the transient evolution of the ice sheet geometry is unclear. This result clearly requires that we acknowledge the limitations of current state-of-the-art 3-D ice sheet models that rely on simplified climate forcing. The predicted borehole temperature profile for the new NGRIP site is also shown in Figure 8. In comparison to the GRIP borehole profile, GrB NGRIP is distinguished by a 63% stronger glacial cold peak and a 30% weaker Holocene warm peak along with a 1.74°C warmer basal temperature (-7.5°). Observationally based inferences indicate that the base of NGRIP is near the pressure melting point, at about -2.5°C (D. Dahl-Jensen, personal communication, 2002) which is thus much warmer than the model prediction. Extrapolation of model temperature to the NGRIP site depth of 3100 m gives a temperature of -2.5 °C for the NGRIP grid cell. Thus the discrepancy in basal temperature is largely if not entirely due to the approximately 170 m thinner ice in the model NGRIP grid cell. The weaker Holocene warm peak is due to the continuous Holocene high surface elevation of the model NGRIP site in contradistinction to the model GRIP site which experiences a mid-Holocene elevation depression of about 140m. During the glacial period, model NGRIP and GRIP site elevation changes were generally quite similar and therefore the colder glacial peak in the NGRIP temperature profile is largely if not solely due to the increasing impact with depth of advected cold ice from higher elevations.
3.2. Constraints on the Eemian Evolution of the Greenland Ice Sheet
 Having described the model and the quality of the fit to the imposed constraints, we will next proceed to use the model to investigate Summit migration, 18O tracer predictions, and the constraint on Eemian Greenland ice sheet evolution.
3.2.1. Summit Migration and Ice Source Elevation
 As 1-D inversions of both the borehole temperature and age profiles at GRIP rely on the assumption of a stationary summit, it is an important check to consider the summit location chronology of tuned 3-D models. As shown in Figure 9, model GrB summit position matches the current location only during the latter parts of interglacials. During glacial periods, the model summit is located a quarter to half a degree to the south and oscillates a half degree (single grid point distance) to the west. This single grid point oscillation may arise solely from the limited numerical resolution. Marshall and Cuffey  find a somewhat similar pattern of summit migration, although with approximately twice the migration distance. As noted by Marshall and Cuffey , such summit displacement can explain the absence of an observed bump in isochronal layers (“Raymond bump” [Raymond, 1983]) under summit which is expected for the case of a stationary summit [e.g., Schott-Hvidberg et al., 1997].
 The cyclical displacement of summit also implies that much of the ice in the GRIP ice core did not originate from that location. Given the 18O to temperature dependence on elevation, ice source elevation needs to be explicitly taken into account in deriving a climate forcing from the GRIP 18O record. We have traced ice source elevation back to the Eemian optimal. Furthermore, to eliminate the influence of the assumption that the ice was sourced from the time-dependent elevation of the GRIP site, we have twice iterated model runs with ice source elevation chronologies from previous runs (we employed the ratio of the source elevation minus GRIP elevation to the difference between the contemporaneous summit elevation and the GRIP elevation to fix the source elevation history for the next iteration). From approximately 70 ka until present, ice originated from very close to the GRIP site as shown in Figure 10. Further back in time, source elevation varies away from that of the contemporaneous GRIP site, bounded of course by the elevation of summit and of the GRIP site (though numerical dissipation in the tracer model occasionally oversteps these bounds). The difference in elevation between summit and the GRIP site remains small until the Eemian. Most important for constraining the Eemian minimum of the Greenland ice mass, source elevation during the Eemian minimum is effectively that of GRIP which is significantly below the contemporaneous modeled summit. This will have a significant impact on the magnitude of the contribution inferred for the Greenland ice sheet to the Eemian sea level highstand (see below).
3.2.2. 18O Tracer Analyses
 It has remained somewhat contentious concerning the extent to which the Eemian section of the GRIP ice core has suffered flow/fabric disturbances and therefore disturbances to the inferred 18O chronologies [Johnsen et al., 1997]. The profiles of certain chemical tracers across the sudden apparent cooling events of the Eemian are difficult to explain if these events were the result of stratigraphic disturbances [Steffensen et al., 1997]. Records from the subpolar North Atlantic suggest that coupled surface-deepwater oscillations occurred just prior to the Eemian interglacial [Oppo et al., 2001]. There is also some far-field evidence from Lake Baikal biogenic silica and microfossil abundance records for a mid-Eemian cooling event [Karabanov et al., 2000]. On the other hand, the methane variations across the Eemian segment of the GRIP and GISPII cores has no counterpart in the Vostok record, suggesting that there is stratigraphic disturbance [Chappellaz et al., 1997]. On the basis of the apparently much more stable Eemian Antarctic climate, Cuffey and Marshall  argue that the Vostok Deuterium record provides a more realistic chronology for the period prior to 100 ka. We have followed their suggestion in splicing the Vostok chronology onto the Greenland record with extremal amplitude adjusted to that of the GRIP record. However, this is problematic in that matching of the extremal amplitude from the last 250 kyr (see the reconstruction labeled DVw in Figure 11) results in a significantly different reconstruction as opposed to matching for the period prior to 105 ka (reconstruction DV). Unless otherwise specified, we will therefore use the DV chronology for subsequent analyses.
 The synthetic GRIP 18O profiles obtained in this way have an excellent match with observations for the glacial period back to just prior to the Eemian as is evident in Figure 12, further validating the quality of the tuned model. However, during the Eemian and also for the period after last glacial maximum, phase differences indicate limitations in the model ice sheet chronology. For this latter period, downward displacement of the phase profile indicates excessive ice accumulation at the start of the Holocene period. This is likely a result of the computationally convenient assumption of thermodynamic control on temporal precipitation change which holds reasonably well for glacial periods. This assumption has been shown not to hold for the early Holocene for which changes in atmospheric circulation appear to dominate temporal changes in precipitation [Cuffey and Clow, 1997]. For the mid-Eemian period, the phase error of the model 18O profile is opposite to that of the early Holocene (Figure 12). This is again likely due to the impact of atmospheric circulation changes on the regional accumulation rate. However, at this depth, we cannot rule out that some phase error may also be arising from limitations of the ice flow chronology, possibly due to the lack of accounting for longitudinal stresses or differences in the evolved geometry of the ice sheet.
 The impact of the assumed 18O lapse rate on the inferred time series of surface 18O variations that is required to drive the tracer model is small relative to other sources of uncertainty and tends to consist of a slight shift to higher values of the downcore 18O values relative to those predicted without adjustment of the lapse rate effect. It must be remembered that the source elevations employed for the 18O inversion are obtained by iterated tracing of ice source location for the model GRIP core. The impact of the assumed 18O lapse rate on sites further from summit will likely be larger.
 Validation of the tuned coupled model and of the tracer module allows a clear test of a possible source of the high Eemian variability of the GRIP 18O record. Comparison of the Eemian profiles for models driven with the raw GRIP 18O chronology and the DV 18O chronology with much lower Eemian variability in Figure 12 shows that with an isotropic flow law model that ignores longitudinal stresses, there is no discernible contribution to 18O Eemian variability arising from model summit migration. Therefore, it is likely that the main source of flow disturbance to the Eemian segment of the GRIP 18O record is due to dynamically induced folding of the stratigraphy arising from anisotropic components of the ice rheology [Dahl-Jensen et al., 1997].
 We have also examined whether the 18O record could be used to constrain the Eemian 18O sensitivity to temperature and thereby indirectly Eemian extent. However, we find no significant difference in the GRIP site 18O tracer record for models using 0.312 and 0.5 values of c for the Eemian period.
 While the tracer 18O profile for the GRIP site has significant phase errors only during the interglacial and terminal glacial periods, the tracer profile for the GISP II site (Figure 13) displays excessive downward phase displacement throughout the core. The GISPII record is indeed displaced downward during the glacial period relative to the GRIP core, but not as much as the tracer model predicts. The tracer 18O chronologies for the GRIP and GISPII sites are very similar and the GISPII site tracer is closer to the inferred GRIP record than to the inferred GISPII record (not shown), indicating the spatial persistence of the forcing chronology when using such a simple 18O forcing for the ice sheet.
 The 18O tracer chronology for the NGRIP site is displaced upward by about 2.1 per mil relative to that for the GRIP site (Figure 14). This is a direct result of the 344 m difference in present-day surface elevation between the two model sites and the 18O lapse rate used to compute the upper boundary condition for the tracer. Aside from this bias, the 18O tracer for NGRIP closely follows the GRIP record with a cold (increasing isotopic depletion) drift toward the Eemian arising from the nearer to summit sourcing of the older (i.e., deeper) ice. It is also worth noting the significantly reduced signal dissipation at the NGRIP site that most likely arises from a longer projected horizontal relative wavelength due to a stronger and more horizontal flow. Further improvements in signal quality will likely arise further downstream until the vertical component of the ice flow starts to increase near the margins. This would benefit possible future studies combining advanced 18O deposition models with coupled ice sheet tracer models. Comparison of model and far from summit ice core 18O profiles could offer powerful constraints on the transient history of existing ice sheets.
3.2.3. Eemian Sea Level Contribution
 Crucial to constraining the minimum volume of the Greenland ice sheet during the Eemian interglacial is the representation of climate forcing assumed in the analysis, or in this case the 18O paleothermometric calibration. If we were to assume that c(Holocene) is close to c(Eemian), the difference in c(Holocene) between models GrB and GrC (both of which fit the GRIP borehole temperature and inferred age profiles) already implies significant uncertainty in the inferred Eemian climate. Furthermore, a more complete probe of parameter space would likely offer an even larger range in possible values of c(Holocene). Other independent constraints on c would therefore be useful.
 A number of processes are likely involved in determining the value of c [Cuffey, 2000]. One-dimensional model-based analyses, for instance, suggest that changes in source temperature and regional evaporative recharge rates can alone explain the difference between observed spatial gradients of 18O (with respect to temperature) and inferred temporal gradients for the Antarctic [Hendricks et al., 2000]. A different physical basis for the small value of cG (during the glacial period) in comparison to the present-day value inferred on the basis of observed spatial and temporal isotopic gradients has been suggested by AGCM-based analyses [Fawcett et al., 1997; Krinner et al., 1997; Werner et al., 2000]. A reduction in c from the warm Holocene to the cold glacial is attributed to the increased seasonality of precipitation during glacial periods with little accompanying change in the temporal form of the seasonal cycle of either the surface temperature or 18O. The 18O record thus becomes a more warm-biased temperature proxy during glacial times. Other possible contributing factors, such as changes in the origin of precipitation and changes in tropical sea surface temperatures, have been found to be much less important in AGCM studies [Werner et al., 2000]. While an increase in present-day mean temperature could further reduce the seasonality of precipitation, the rate at which this takes place as a function of temperature change is likely to be much less than for colder climates. For glacial climate, GCM simulations indicate a much more zonal winter circulation which drastically reduces the transport of moisture to the ice sheet [Werner et al., 2000]. There is no direct evidence nor any apparent physical basis for suggesting a similar change in atmospheric circulation with warmer climates typical of the mid-Holocene or Eemian. The impact of changes to ocean water 18O will also be much smaller for warmer than present as compared to colder than present climates. In contradiction therefore to recent 1-D model-based borehole temperature inversions for the Holocene period, one would expect c for periods warmer than present to be no smaller than that for present-day and to also be larger than cG. Therefore, we base our upper bound for c(Eemian) on the present-day bound for c. Shuman et al.  have measured the present-day (seasonally based) value of c to be most likely between 0.4 and 0.5, with 95% confidence intervals of 0.34 to 0.68 [Shuman et al., 1995]. The 1-D model-based borehole temperature inversion of Cuffey et al.  also required a value of c of 0.47 for the last 500 years. The older and more poorly fitting paleothermometric calibration of Johnsen et al.  has an c value of 0.6 for high 18O values that characterize the Holocene. Given the ranges of these previous analyses, we choose a value of 0.6 as a defensible upper bound for c(Eemian).
 As a lower bound for plausible values of c(Eemian), we initially choose the value of c(Holocene) = 0.25 from the analyses of Cuffey and Clow . Rerunning the constrained model using this range of Eemian c, we obtain a range of minimum ice volumes from 2.23 to 1.12 × 1015 m3 for model GrB when the Vostok DV 18O chronology is assumed. Extracting sea level contributions is further complicated by the 7% excess present-day ice volume in the model. Figure 15 presents estimated eustatic sea level contributions reduced by the ratio of the present-day excess model volume (i.e., by 1/1.07).
 On the basis of the total gas content record of the GRIP core [Raynaud et al., 1997], Cuffey and Marshall  argue that the minimum Eemian elevation of the GRIP core ice was no lower than about 2900 m. As indicated by the boxes in Figure 15, this provides an upper bound to maximum Eemian sea level contributions of about 4.4 m and 3.9 m, respectively, for model GrB with accounting for GRIP ice source elevation and model GrC (without such accounting). Lower bounds for the models are near 2 m. Different model precipitation sensitivities to temperature change likely accounts for the model tuned with the low value of c(Holocene) (=0.25, GrC) also having the lowest sea level contribution for the lower boundary value of c(Eemian). The significantly different slope of the GrC sea level sensitivity curve suggests even further caution in interpreting model-based constraints upon the volume of the Greenland ice sheet during the Eemian interglacial.
 One way in which we might better quantify the impact of the entire set of constraints employed to construct model GrB is to compare the sea level sensitivity curve for the untuned version of GrB (Figure 15, “untuned model”, with DV 18O chronology) to that for the fully tuned model. This “untuned” model was only tuned to approximate present-day ice volume and elevation at summit (using only flow enhancement and sliding parameters). RSL dates, present-day topography, and borehole temperature and age profile constraints were thereby ignored. This relatively unconstrained model predicts approximately an extra 0.5 m of excess Eemian sea level contribution arising from the prolonged period of high climate variability.
 While the predicted Eemian sea level excess for Greenland is most sensitive to the assumed value for Eemian isotopic sensitivity other factors are not insignificant. Use of the 18O chronology with high Eemian variability (i.e., the unmodified GRIP record) for the climate forcing results in an extra 0.6 to 0.7 m sea level contribution. Given its dependence upon the vagaries of midlatitude storm track migration, it is likely that the true Eemian Greenland climate resides somewhere between the extremes represented by the GRIP and Vostok DV chronologies. The assumed source elevation of the GRIP site ice is also important. The assumption of GRIP ice sourced to the contemporaneous summit can result in more than a meter of extra sea level contribution as compared to the assumption of locally sourced ice and more than 0.7 m excess relative to models that take into account the actual (model) source elevation of the GRIP site ice.
 One additional uncertainty requires consideration. As detailed by Tarasov and Peltier , our model of Greenland ice sheet evolution required significant modifications to the near coastal Holocene temperature forcing in order to match computed and observed relative sea level observations. We suspect that a significant fraction of this extra forcing is accounting for the lack of explicit fast flow mechanics in the ice sheet model. Whatever this extra forcing represents, the apparent need for this enhanced forcing during the Holocene would suggest that the deglaciation event predicted under the assumption of a simple climate forcing based upon a single proxy may well underrepresent the diminution of ice volume that actually occurred during the Eemian. During the Holocene minimum at 8.5 ka, differences in ice volume between the fully tuned model and the model lacking the modified Holocene regional temperature forcing were 3.205 × 1015 - 2.925 × 1015 = 0.28 ×1015 m3. Subsequent to this, differences diminished. If we simply take this as an additional uncertainty in our estimation of the maximal excess Eemian sea level contribution, this adds 0.7 m of possible additional sea level contribution.
 Most problematic is the reliance on a single climate proxy from the summit of Greenland to derive a climate chronology for the whole of Greenland. Recent observations find anticorrelations between local climatic responses for the east and west coastal regions [White et al., 1997]. Yet it is precisely the near-margin climate that will largely determine the extent of Greenland ice. The apparent absence of deep water formation in the Labrador Sea during the Eemian interglacial [Hillaire-Marcel et al., 2001] would have limited oceanic heat transfer to the region, suggesting limited regional warming. On the other hand, this situation might have been caused by a large freshwater flux from Greenland, suggesting fast and significant deglaciation. It is therefore important to investigate other possible constraints on the extent of Eemian Greenland deglaciation.
 Inferences on the marginal extent of the Eemian Greenland ice sheet might be obtained from ice cores located in regions that were potentially deglaciated during the Eemian. Considering Camp Century and Dye 3, for instance, Camp Century retains ice even for c(Eemian) = 0.312 (Figure 16), while Dye 3 is barely covered by the margin of a residual dome for c(Eemian) = 0.6 and is fully deglaciated with c(Eemian) = 0.5. For all cases, the southern dome is cutoff from the rest of the ice sheet at the time of minimal Eemian extent, in contrast to the Holocene minimal extent at -8 kyr. In the model, the ice divide upon which Camp Century resides is much more robust against Eemian warming than the southern dome upon which Dye 3 resides. Koerner  presents some evidence for Eemian ice-free conditions at both Dye 3 and Camp Century but no definitive conclusion. If observations could provide strong independent evidence for bedrock exposure at Camp Century during the Eemian, this would tend to imply a very strong Eemian deglaciation and a much warmer climate. On the other hand, our analyses do suggest that Dye 3 was ice-free during the Eemian.
 Independent observational constraints on Eemian climate are generally lacking. However, macrofossil analyses of a till-covered sequence in Washington Land, northwestern Greenland [Bennike and Jepsen, 2000], suggests that the Eemian interglacial climate was not too different from modern. Given that the regional sea level adjusted warming during the peak of the Eemian is a significant 6.8°C relative to present-day for c(Eemian) = 0.4, this data favors a larger value of c(Eemian) (and therefore limited warming).
 The above model results and observations suggest that a reasonable upper bound for Eemian Greenland eustatic sea level contribution is about 5.2 m (2900 m GRIP Eemian elevation limit for model GrB, with contemporaneous ice source inversion plus 0.7 m uncertainty arising from Holocene model tuning as discussed above), while a lower bound is no more than 2 m. However, given the numerous sources of uncertainty, even wider bounds are possible. More likely values are arguably bounded by c(Eemian) = c(Holocene) (or c(Eemian) corresponding to 2900 m minimum GRIP site elevation, whichever constraint is stronger) and by the present-day observed value of c of approximately 0.5. This corresponds to a range of about 2.7 m to 4.5 m for the set of models GrB and GrC (with 0.6 extra meters added to GrC to allow for source elevation discrepancies and regional variations in climate). This contrasts with the 4 to 5.5 m likely range estimate of Cuffey and Marshall , who use a stricter upper bound for c(Eemian) of 0.43 based on central Greenland borehole temperature analyses from the last 5 kyr and who set their lower bound to the 2900 m elevation constraint. Some of the range difference is also due to different model sensitivities and assumed contemporaneous summit sourcing of GRIP ice (though Cuffey has indicated that their model has low sensitivity to assumed ice source location (K. M. Cuffey, personal communication, 2001)). For instance, at c(Eemian) = 0.36, model GrB (DV 18O chronology) has Eemian sea level contributions of 3.8 m and 4.5 m for respectively source-traced and summit-traced ice, while Cuffey and Marshall  have almost a 5 m contribution using summit-traced ice. This increased Eemian deglaciation is likely due to the stronger slope-dependent feedback in the precipitation factor that Cuffey and Marshall  employed which tends to reduce summit elevation changes and thereby allow stronger inferred climate changes for a given c(Eemian) (refer to equation (6)). In order to match borehole temperature and age profiles for the GRIP site, we found it necessary to weaken the slope-dependent factor for higher elevations (equation (5)). Given the more complete set of constraints imposed upon the model, our results suggest that the 5–6 m excess eustatic rise generally assumed to characterize the penultimate interglacial [e.g., see Rostami et al., 2000] likely also involved mass loss from the south polar cryosphere.
Citation: Greenland glacial history, borehole constraints, and Eemian extent, J. Geophys. Res., 108(B3), 2143, doi:10.1029/2001JB001731, 2003.
Copyright 2003 by the American Geophysical Union.