2.1. Datasets

Reference Ohmura and ReehOhmura and Reeh (1991) compiled a Greenland accumulation dataset (hereafter OR), combining direct snow-gauge observations with estimates from glaciological measurements (Fig. 2a). They found the average accumulation rate over the ice sheet to be 0.31 m a⁻¹, with a maximum of around 1.5 m a⁻¹ on the southeastern flank. Broadly speaking, the accumulation rate decreases with elevation and with latitude: the colder the temperature of the air column, the less moisture it is able to hold. Secondly, the local maxima and minima in the accumulation rate are attributable to the atmospheric circulation over Greenland (Reference Ohmura and ReehOhmura and Reeh, 1991). The low-level winds are generally anticyclonic, with some variation in strength and position over the course of the year. Where the prevailing winds are forced upslope on the ice sheet, there is enhanced precipitation because the rising air column cools and becomes saturated. Conversely, where prevailing winds are down-slope, precipitation is decreased: the so-called “rain shadow” effect observed in the lee of mountain ranges (e.g. Reference SmithSmith, 1979). There are also seasonal and regional variations in storminess over the ice sheet. During winter, for example, many storms impinge on the southeastern flank of the ice sheet, contributing to the high precipitation rates there (e.g. Reference Chen, Bromwich and BaiChen and others, 1997).

Fig 2.

Different accumulation and precipitation datasets for Greenland. (a) Annual accumulation ratefrom Reference Ohmura and ReehOhmura and Reeh (1991); (b) ERA mean annual precipitation for 1985–93; and (c) Reference Chen, Bromwich and BaiChen and others (1997) annual precipitation rate diagnosed from MCEP re-analysis winds and temperatures (CBB). For clear comparison and to show interior accumulation detail, the color scale in each figure is truncated at 1 m a⁻¹.

Ohmura and Reeh associate an accuracy of ±20% with the individual measurements used to produce Figure 2a. Because of this uncertainty, as well as the poor spatial and temporal coverage of the observations, there have also been attempts at some indirect, model-based estimates of the precipitation. Figure 2b shows the mean precipitation for the years 1985–93 from the European Centre for Medium-Range Weather Forecasts re-analysis project (hereafter ERA; Reference Gibson, Hernandez, Kållberg, Nomura, Serrano and UppalaGibson and others, 1996). The ERA re-analysis is a product of a numerical weather-prediction model which has meteorological observations assimilated into it as part of its operational cycle. The precipitation field is determined mainly by the numerical model rather than the direct observations. Since temperatures generally rise above zero only in summer, and then only over a relatively small fraction of the ice sheet, the precipitation rate in Figure 2b may be compared to the accumulation rate in Figure 2a. While the broad features of the distribution are the same, the ERA precipitation rate over the southeast is significantly higher than the OR precipitation rate, and the interior accumulation rate is much lower than observed — a recognized deficiency in the ERA data (Reference Bromwich, Cullather, Chen and CsathóBromwich and others, 1998).

Reference Chen, Bromwich and BaiChen and others (1997) made an independent estimate for precipitation over Greenland by diagnosing the vertical velocity in σ-coordinates using the U.S. National Centers for Environmental Prediction (NCEP) re-analysis dataset for horizontal winds, temperatures and relative humidity. Figure 2c is their estimate for the years 1985–95 (referred to hereafter as CBB). It is similar to the ERA data in that interior values are significantly underestimated. It also fails to show a local maximum in west central Greenland that appears in both the OR observations and the ERA re-analysis. Figure 2c does, however, show another observed local maximum in precipitation in northwest Greenland (close to Thule), not seen in the ERA.

2.2. Simple parameterization of precipitation

Here we outline a parameterization for precipitation that was used in Reference Roe and LindzenRoe and Lindzen (2001) for ice-age studies using a simple climate model. It was designed to represent the processes causing precipitation that were outlined in section 2.1. In section 3 it will be used to generate accumulation patterns to compare with the datasets already presented, and will be used to identify features of the distribution that are important for the response of an ice sheet.

The parameterization determines precipitation rates from time-mean (e.g. seasonal or monthly) wind and temperature fields, such as might be generated from a simple climate model. We assume that these fields are averaged over a time period that is much longer than the day-to-day synoptic variability. This variability is accounted for by including a probability distribution of vertical velocities, centered around the time-mean vertical velocity (method of calculation is described below). For the fraction of time for which the vertical velocity, w, is between w′ and w′ + dw′, the precipitation, dP, is taken to be

(1)

f(w′) is the probability of w being between w′ and w′ + dw′, and is defined below. a and b are constants. e _sat is the saturation vapor pressure at the surface. It is given by the Clausius–Clapeyron equation (e.g. Reference HoltonHolton, 1992), and is a good proxy for the moisture content of the overlying air column. It is an exponential function of the surface temperature, T _s : e _Sat(T _s)= e ₀ exp[c ₁ T _s/(c ₂ + T _s)], where e ₀ = 6.112 mbar, c _i = 17.67 and c ₂ = 243.5, and T _s has units of °C. We take f (w′) as a Gaussian distribution centered on the time-mean vertical velocity, w ₀:

(2)

N is a normalization factor to ensure that . The measure of variability of the vertical velocity, α, was taken to be 1.15 cm s⁻¹, which was the average value over the region at 850 mbar, determined from daily NCEP data. A possible refinement of the parameterization would be to allow for α to vary with position (e.g. with proximity to the storm track). However, the results of section 3 show that this would probably have little impact on the resulting shape of the modeled ice sheet.

The precipitation rate for a given surface temperature is found by integrating Equation (1). The integral can be reduced to an expression involving error functions for which accurate asymptotic expansions are available (see Reference Roe and LindzenRoe and Lindzen, 2001).

Equation (1) can be related to the equation for conservation of water mass (Reference RoeRoe, 1999). This analysis shows that the parameter a roughly represents a background precipitation rate that would occur in the absence of time-mean vertical motion. A value for b can also be derived from this analysis, from which we take b = 5.9 × 10⁻⁹ s⁻³ kg⁻¹. This value is similar to that used in a related parameterization (Reference Sanberg and OerlemansSanberg and Oerlemans, 1983; see also Reference RoeRoe, 1999), which was determined by tuning to observations over Scandinavia. That both approaches lead to similar numbers gives confidence that the parameterization yields physically plausible values.

The vertical velocity was determined by two separate methods. The first was to use the NCEP re-analysis dataset of pressure vertical velocity, ω, available on constant pressure surfaces on a 2.5° by 2.5° horizontal grid. On daily time-scales and longer, a good approximation is w ₀ = − (H/p)ω, where p is the pressure and H is the atmospheric scale height (∼8.5 km). The second method was to calculate the vertical velocity using the relationship , where z _s is the surface topography and is the horizontal wind (also taken from the NCEP re-analysis). This treatment assumes that all of the vertical motion is forced by flow over the topography. Over an ice sheet, however, there is a significant component which is due to air sinking over a cold surface. This second formulation therefore assumes that this sinking motion can be represented by its areal average and subsumed into the parameter a of Equation (1). The atmospheric data were interpolated onto the 40 km horizontal grid of the ice-sheet model (see section 3).

In the model integrations presented below we use atmospheric data from the 850 mbar pressure surface, which has a mean height of about 1400 m. As seen from Figure 1, this level coincides with the steepest slopes on the ice sheet, and is therefore the most appropriate for the orographically induced accumulation. Where this pressure level lies below the physical surface of the ice sheet, the NCEP re-analysis gives values of and ω for directly above the surface. Data from the 1000, 925 and 700 mbar surfaces were also used to generate accumulation fields, and the differences were small. As long as the atmospheric fields change only slowly with height this will be true. The exception to this was that when data from the 700 mbar surface (average height 2800 m, generally well above the ice-sheet surface) were used, the accumulation rate on the southwestern flank of Greenland was excessive (compared to OR) due to the larger westerly winds at higher altitudes.

The accumulation rate for Greenland was determined from Equation (1) using the observed topography and monthly-mean NCEP re-analysis data, averaged between 1985 and 1995 (for comparison with CBB), and using the NCEP 2 m temperature for T_s . Precipitation was deemed to fall as snow if the monthly-mean T _s was < 0°C, although the results were not sensitive to this choice. Any rainfall was assumed to run off the ice sheet. Monthly accumulation rates were summed to give the annual accumulation rate. Note that there is an implicit smoothing associated with the 40 km grid spacing, and also from using the 2.5° by 2.5° NCEP re-analysis data.

The focus of this work is to compare an ice sheet’s response to different patterns of accumulation, and so the accumulation rate was calculated as described above, adjusting the value of a in Equation (1) until the mean accumulation rate over the ice sheet agreed with the OR observations (0.31 m a⁻¹).

3.1. Model integrations using accumulation datasets

For all of the runs presented, the model is initiated with the current Greenland topography and integrated for 100 kyr, which is sufficient time for the model to come to equilibrium with the forcing climate. A summary of all the results is presented in Table 1.

Table 1.

Summary of ice-sheet model integrations

The model integration using the OR-observed accumulation produces a good simulation of the observed ice sheet (Fig. 3a). Differences are that the maximum height is slightly overestimated (by 40 m), and the modeled ice sheet is overextended slightly in the northwest of Greenland. Another minor discrepancy is that the axis of the southern dome has a slightly different orientation. Overall, though, the area, volume and maximum height of the ice sheet are well simulated.

Fig. 3.

Result of integrating SICOPOLIS model using OR accumulation data, (a) Surface topography, contours every 200 m; (b) contours and vectors of vertically integrated mass flux. Contour interval 60 km m a⁻¹

The vertically integrated mass flux (Fig. 3b) also reflects the overlying accumulation distribution, although not directly. In steady state, the divergence of the mass flux is equal to the local net mass balance (i.e. including surficial and basal melting). Furthermore, probably the most important determinant of the mass flux is the basal topography: high mass flux occurs over low basal elevations and within bedrock channels (see, e.g., Reference Letréguilly, Huybrechts and ReehLetréguilly and others (1991) and Fig. 5b (shown later), a model integration using uniform accumulation). However, if all other factors stay the same then a higher accumulation rate drives a larger mass-flux divergence, and leads to a larger mass flux. This relationship between ice flux and accumulation has important consequences for the internal ice dynamics. Regions of high mass flux produce greater internal frictional heating which can raise the basal temperature to the melting point. When this occurs, the ice is no longer frozen to the bedrock but becomes free to slide. Areas of the ice sheet with temperate ice at the bottom coincide closely with the contoured regions of high mass flux in Figure 3b.

The CBB dataset was available as monthly precipitation rates. To convert it to an annual accumulation rate it was assumed that precipitation fell as snow if the NCEP re-analysis mean monthly 2 m temperature was < 0°C, and the monthly rates were summed over the annual cycle. This processing only affects the accumulation rate at low elevations in southern Greenland, and has little impact on the modeled ice sheet. The resulting pattern was also multiplied by 0.97 so that the average accumulation rate was the same as the OR observations. The equilibrium ice sheet obtained by integrating the ice-sheet model with this accumulation pattern is shown in Figure 4a. The most obvious difference from Figure 3a (the OR accumulation) is that the lower accumulation rate in the northern interior causes an underestimate of the volume and maximum height. Moreover, the peak of the northern dome sits much closer to the coast at a location where there is only a secondary maximum in the observations. Differences in the accumulation pattern are also reflected in the mass flux. The weaker accumulation rate in North Greenland produces a much weaker mass flux there (Fig. 3b). In southern Greenland the mass flux agrees fairly well in the two integrations.

Fig. 4.

As for Figure 3, but for SICOPOLIS integration using CBB accumulation data.

A model integration was performed using the ERA precipitation dataset, processed to give an accumulation rate in the way described above. It produced an ice sheet that looked very similar to that using the CBB dataset, and the results are summarized in Table 1.

Caution is required in comparing these model integrations to the observed ice sheet. The dynamical adjustment time-scale for an ice sheet the size of Greenland is several thousand years, and in general, therefore, it is not going to be in equilibrium with the climate forcing at any given time. However, paleoclimate records in the region indicate that the climate has been comparatively stable in the region over the last 10 kyr (e.g. Reference DansgaardDansgaard and others, 1993), so we might expect the Greenland ice sheet to be fairly close to equilibrium. Addressing this question, Reference Ritz, Fabre and LetréguillyRitz and others (1997) performed a time-dependent simulation of the Greenland ice sheet using a plausible reconstruction of the climate forcing over the last glacial cycle. They compared the resulting ice sheet for the current time to an ice sheet that was in steady state with the current climate forcing. For a range of parameters used in the ice model, the maximum height and volume in the steady-state integration were about 5% and 10% smaller, respectively, than in the transient run. Reference GreveGreve (1997) found similar results.

Because the CBB and ERA distributions have known deficiencies (Reference Bromwich, Cullather, Chen and CsathóBromwich and others, 1998), and because of the good agreement between the observed ice sheet and that modeled using the OR accumulation (especially their overall shape), the model integrations in this section support the OR observations as being closer to the actual accumulation distribution than either the CBB or ERA datasets. Therefore we shall compare the model integrations produced using the accumulation parameterization presented in section 2.2 to Figure 3a and b.

3.2. Model integrations using accumulation parameterization

We now examine what features of the accumulation pattern are necessary to model the Greenland ice sheet. The simplest possible pattern is a uniform accumulation field equal to the average of the OR observations (0.31 m a⁻¹). The equilibrium ice sheet from an integration using this forcing is shown in Figure 5a. While the area is close to the observed value, both the volume and maximum height are significantly overestimated (Table 1): too much snow accumulating in the interior leads to an excessive build-up of the ice there. Comparing to the integration using OR observations, the differences in ice mass flux reflect the differences in accumulation rates (Fig. 5b). The mass flux is significantly over-predicted in northern and central Greenland, while it is under-predicted along the southern coast.

Fig. 5.

As for Figure 3, but for integration using uniform accumulation field.

Next an accumulation pattern is generated only from the surface temperature observations. The procedure outlined in section 2.2 followed, but with the parameter b set to zero. Thus the winds do not influence the accumulation distribution, which includes only the effect that decreasing temperature with elevation and latitude reduces the ability of an air column to hold moisture. The parameter a in Equation (1) is adjusted so that the average accumulation rate over the ice sheet is the same as in the OR observations, and the resulting pattern is shown in Figure 6a. The accumulation rate in the Greenland interior is slightly overestimated, and while there is a band of maximum accumulation around the southern coast, it is not as large as in the OR observations. Figure 6a confirms, as was argued in section 2.1, that it is the atmospheric circulation that generates these features.

Fig. 6.

(a) Accumulation derived using simplified parameterization with surface temperature dependence T _s only; (b) surface height, contours every 200 m; and (c) contours and vectors of vertically integrated mass flux. Contour interval 60 km m a⁻¹.

Using this accumulation forcing, the SICOPOLIS model gives the equilibrium ice sheet shown in Figure 6b. The area, maximum height and volume all agree quite closely with observations (Table 1). The main discrepancy with the integration performed using the OR observations is the ice mass flux, and again there is a clear connection to the accumulation-rate differences. The mass flux in the northern half of the ice sheet is overestimated, and around the southern margin it is underestimated (Fig. 6c).

The next step is to include the effect of the atmospheric circulation on the accumulation distribution. The accumulation field determined from Equation (1), using and data from 850 mbar, is shown in Figure 7a. The gridscale variability is attributable to having interpolated the coarse-resolution atmospheric data onto the finer grid on which the topography is calculated. Comparing with the OR observations (Fig. 2a), the localized maximum on the southeast coast is seen (although substantially over-predicted), as is the minimum in the southwest. The local maximum near Thule also appears, although again over-predicted. In the interior, the accumulation rate is about equal to the OR observations. Overall, while the large-scale features are fairly well captured, there is a lot of small-scale variation, and exaggeration of localized maxima and minima.

Fig. 7

As for Figure 6, but including dependence on atmospheric circulation (i.e. using ).

The integration of the SICOPOLIS model with this accumulation is shown in Figure 7b. It is notable that even though the accumulation field has substantial gridscale noise, the resulting ice sheet is much smoother. The diffusive nature of ice flow is highly effective in smoothing out small-scale variations in the accumulation rate. The overall shape of the ice sheet resembles the observed, except in southwest Greenland where the underestimation of accumulation rate has caused the ice sheet to thin and recede.

The regions of high mass flux within the ice sheet (Fig. 7c) agree well with those obtained using the OR accumulation, although there is some overestimate of the ice flux in the southeast where there is excessive accumulation. That the mass flux agrees more closely than either the uniform or temperature-dependent accumulation patterns demonstrates that capturing the interaction between the prevailing winds and the surface topography is necessary to successfully model the underlying ice dynamics.

An accumulation pattern was also generated using the NCEP re-analysis ω to calculate ω ₀ in Equation (1). This approach produced an accumulation pattern that picked up the circulation features seen in the OR observations, although smeared out. When it was used as climate forcing, the SICOPOLIS model reproduced the shape of the ice sheet (Table 1). The modeled mass flux, however, did not do as well as that seen in Figure 7c in capturing the variations seen in the integration using the OR accumulation.

In summary, the overall shape of the Greenland ice sheet is quite robust to different patterns of accumulation. Only in the model integration forced by uniform accumulation did the error exceed 10% (Table 1), and then only in volume. So the results suggest that provided the accumulation rate decreases with surface temperature, following the Clausius–Clapeyron equation, the model integration will produce an ice sheet resembling the modern one. However, the mass flux within the ice sheet is a more sensitive measure of the local accumulation distribution, and reflects quite strongly the differences in accumulation pattern caused by the prevailing atmospheric circulation.

As already mentioned, the ice mass flux is closely connected with an ice sheet’s basal boundary condition, changes to which have been invoked as an explanation of the periodic discharge of large parts of the Laurentide ice sheet during Heinrich events (Reference MacAyealMacAyeal, 1993). An implication of the results here therefore is that the successful modeling of such events in climate models likely requires knowledge of the atmospheric circulation controlling the accumulation distribution over the ice sheet.

3.3. Sensitivity to mean accumulation rate

The accumulation patterns used in section 3.2 all had the same mean value. To test the sensitivity of the ice sheet to changes in this value, a series of integrations were performed using the OR accumulation multiplied by a factor varying from 0.5 to 2.0. For this range, the area and maximum height of the equilibrium ice sheet varied by ±10%, and the volume by ±20%. This relative insensitivity of the ice sheet to mean accumulation rate is consistent with simple models of ice flow (e.g. Reference PatersonPaterson, 1994). A coastal margin is fixed by calving processes, and a land-based margin is constrained because ablation is a strong function of temperature: an increase in accumulation rate causes only a small advance in the ice-sheet margin before it encounters warm enough temperatures for the ablation to balance the increased accumulation.