Introduction
Glaciers and ice sheets continuously adjust their dimensions in response to current and past climatic change. Long-term (>103 years) changes of the surface elevation are primarily due to changes in ice dynamics reflecting past mass-balance states of the ice sheet. Short-term decadal or seasonal changes are caused by interannual or seasonal variations of weather conditions such as snow accumulation and surface air temperature (e.g. Reference Braithwaite, Laternser and PfefferBraithwaite and others, 1994). Changes in surface elevation of polar ice sheets are of great importance to mass-balance studies since changes in ice volume may be determined from changes in ice thickness obtained from surveys of changes in surface elevation by satellite altimetry (e.g. Reference Davis, Kluever and HainesDavis and others, 1998; Reference Shepherd, Wingham, Mansley and CorrShepherd and others, 2001; Reference Zwally, Brenner and FuZwally and Brenner, 2001). Surface elevation changes are equivalent to ice-thickness changes minus the vertical motion of the bedrock, which is generally smaller and can be separately estimated. However, short-term changes caused by variations in rates of near-surface firn compaction must also be accounted for (Reference Arthern and WinghamArthern and Wingham, 1998).
In dry-snow zones, seasonal and interannual changes in surface elevation may be caused by variations in the rates of accumulation and firn densification. The rate of densification is also affected by the firn temperature and rate of accumulation. Field observations in West Greenland show that variations of near-surface firn density can cause annual surface elevation changes of the order of ± 10–20 cm (Reference Braithwaite, Laternser and PfefferBraithwaite and others, 1994). Reference McConnell, Mosley-Thompson, Bromwich, Bales and KyneMcConnell and others (2000) have recently modeled accumulation rates based on data from 11 ice cores located near the 2000 m contour of the Greenland ice sheet. They suggest that the ice-sheet elevation varies by tens of centimeters from year to year simply because of changing snow accumulation.
Ice-sheet surface elevations derived from satellite radar altimetry have exhibited seasonal variations (Reference Zwally, Brenner, Major, Bindschadler and MarshZwally and others, 1989; Reference Yi, Bentley and StenoienYi and others, 1997), but the specific cause of these seasonal variations in dry-snow zones has been uncertain (Reference Yi, Bentley and StenoienYi and others, 1997). In the percolation areas of Greenland, an observed seasonal cycle (2 m peak-to-peak at 1200–1700 m elevation) with a maximum in mid-April just before the melt season and a minimum in mid-October is consistent with a lowering caused by a combination of summer melting and melt-induced firn compaction (Reference Zwally, Brenner and FuZwally and Brenner, 2001). Also, the observed elevations in the ablation zones show similar summer minima that can be associated with summer melting. However, in the dry-snow zones, observed amplitudes of several tens of centimeters, with a minimum typically occurring in summer, are generally larger than expected from seasonal variations in snowfall alone. Other potential causes for seasonal variations in Geosat altimetry data discussed by Reference Yi, Bentley and StenoienYi and others (1997) include unresolved satellite orbit errors, variations in the subsurface vs surface radar backscatter and variations in firn densification. However, examination of European Remote-sensing Satellite (ERS) orbit quality over oceans eliminates orbit errors as a cause. Also, the observed seasonal cycles with a winter maximum have not been explained by seasonal variations in the effective radar-back-scatter depth. If anything, new fine-grained winter snow should allow more radar penetration and cause detection of a lower surface, giving a winter elevation minimum in contradiction to the observations. Furthermore, calculations of a seasonal cycle in firn compaction using a densification model and similar to Reference Herron and LangwayHerron and Langway (1980) or Reference Arthern and WinghamArthern and Wingham (1998) gave a small seasonal amplitude of only a few centimeters (personal communication from M. Spencer and R. Alley, 2000).
The temperature dependence of densification rate from the constitutive relation used in previous models is based on the Arrhenius equation, with an exponential function of both activation energy and temperature, and a multiplicative rate constant that is independent of temperature. Although the Arrhenius formulation is widely used to describe various processes in snow and ice, including ice deformation (Reference PatersonPaterson, 1994) and grain growth (e.g. Reference GowGow, 1969), different values of the activation-energy “constant” are typically used for different processes or for different regimes of the same process. For example, Reference Herron and LangwayHerron and Langway (1980) used different rate constants and activation energies for density less than or greater than 550 kg m−3 in their densification model. Reference Arthern and WinghamArthern and Wingham (1998) used different formulations and constants at different depths in their densification model to describe the dominant deformation processes such as grain-boundary sliding (e.g. Reference AlleyAlley, 1987a) or dislocation creep (Reference Wilkinson and AshbyWilkinson and Ashby, 1975). The relevant parameters were determined by fitting field density profiles assuming accumulation rate and temperature are both constant (Reference Herron and LangwayHerron and Langway, 1980; Reference AlleyAlley, 1987a). However, this is not the case in practice, because snow precipitation varies in time, and temperatures in upper firn layers vary greatly during the year. Also, measured density profiles exhibit large variations around a generally increasing density with depth. Some of these variations appear to be seasonal and have been attributed to variations in deposition or various processes in the firn (e.g. Reference GowGow 1968; Reference Koerner and CraryKoerner, 1971; Reference AlleyAlley, 1987b; Reference Jiankang and YangHan and Yang, 1988; Reference Dahe and YoungQin and others, 1988; Reference Gerland, Oerter, Kipfstuhl, Wilhelms, Miller and MinersGerland and others 1999).
In our model, we take into account the seasonal cycle of temperature in the upper firn, which is important because of the non-linear dependence of temperature even in the standard Arrhenius-type formulation. More importantly, however, we use a temperature-dependent activation energy and a temperature-dependent rate constant based on the crystal-growth studies of Reference Jacka and LiJacka and Li (1994). Using the detailed records of accumulation and temperature from an automatic weather station (AWS) at the Greenland summit (Reference Steffen, Box and EstupinenSteffen and others, 1999), and atmospheric model results on precipitation variations (D. H. Bromwich and others, unpublished information), we model the densification rates of firn and obtain elevation changes for comparison with the satellite altimeter measurements. Our reference location for this study is 72°34′ N, 38°27′ W (Summit), which is the approximate location of Greenland Ice Sheet Project 2 (GISP2) and the AWS about 25 km west of the actual summit. The radar altimeter measurements are averages over an area of 100 km radius around this location.
Surface Conditions
A continuous surface daily mean air-temperature record for the Greenland summit region over the period 1987–99 has been constructed by Reference Shuman, Steffen, Box and StearnsShuman and others (2001) from a composite of data from three AWSs since 1995 and satellite passive-microwave data from the same locations. Figure 1 shows their smoothed daily temperature curve. The surface air-temperature data over the period April 1992–April 1999 show an average maximum of about −6°C, with an annual mean of about −29°C. Both mean summer (June–August) and mean annual temperatures show a general increasing trend and an interannual variability of about 5°C, with a minimum in 1992 and a maximum in 1995. Hourly AWS data for 1996–99 (Reference Steffen, Box and EstupinenSteffen and others, 1999) show an average diurnal amplitude of about 9°C.
Variations of surface air temperature, accumulation rate and snow surface height at the Greenland summit ice sheet, April 1992–April 1999. The temperature data were compiled from AWS and passive-microwave measurements by Reference Shuman, Steffen, Box and StearnsShuman and others (2001). Accumulation rates are from atmospheric (P − E) modeling results from D. H. Bromwich (unpublished information) normalized to the mean accumulation rate at the summit, and from AWS sonic measurements of surface height by Reference Steffen, Box and EstupinenSteffen and others (1999). Surface elevation time series is derived from ERS-1/-2 satellite radar altimeter measurements with multivariate linear and seasonal sinusoidal fit.
Snow accumulation data have been collected from AWSs for the Greenland summit between May 1995 and August 1998 using the sonic technique (Reference Steffen, Box and EstupinenSteffen and others, 1999). A sonic sensor continuously measures the distance between the snow surface and the sensor on the AWS tower. Decreases in distance are caused mainly by snow accumulation (snowfall and drift), and increases in distance by evaporation and wind erosion. The regression line through the measured distance vs time gives a surface rise of 62.5 cm a−1, which, multiplied by the near-surface firn density, gives the accumulation rate. Averaging over the firn density profile to 62.5 cm depth, which corresponds to the firn added in 1 year, gives an average firn density of 400 kg m−3 (compared to a surface density of 332 kg m−3) and an annual accumulation rate of 250 kg m−2 a−1. This density estimate attempts to account for the effects of new-snow compaction for an AWS pole tower inserted a few meters in the firn. The accumulation rate obtained is close to the 255 kg m−2 a−1 from Reference Zwally and GiovinettoZwally and Giovinetto (2000) and slightly higher than the 220 kg m−2 a−1 from Reference Ohmura and ReehOhmura and Reeh (1991). The accumulation time series derived from the AWS in Figure 1 shows two sudden increases of snow surface height of approximately 40 and 20 cm in January 1995 and July 1997 caused by large snowstorms.
Recently, D. H. Bromwich and others (unpublished information) re-evaluated the precipitation–evaporation (P – E) distribution over the Greenland ice sheet over the period 1985–99, using the re-analysis of atmospheric model data. A pronounced downward trend in annual precipitation for the period 1985–95 was followed by a marked increase for the period 1995–99. However, the model-derived P − E for the Greenland summit is about 50% smaller than the accumulation rate derived from the AWS sonic measurements. Therefore, we normalize their modeled snow precipitation for the Greenland summit by using the mean accumulation rate (250 kg m−2 a−1) derived from AWS measurements over the period April 1992–April 1999, as shown in Figure 1. Interestingly, the modeled precipitation also shows significant peaks in the fall of 1995 and the spring of 1997 corresponding to the storms observed in the AWS data, albeit with smaller amplitudes.
The summit region of the Greenland ice sheet is in the dry-snow zone due to its negligible summer melting (Reference BensonBenson, 1962). The density profile (Reference GowGow and others, 1997) for the summit region shows a firn–ice transition depth around 75–77 m where the firn density reaches approximately 830 kg m−3. During May 1987, a number of shallow cores (to approximately 17 m depth) were recovered over a 150 × 150 km survey grid in central Greenland (Reference Bolzan and StrobelBolzan and Strobel, 1994). Detailed density measurements were made in depth intervals of < 0.1 m (personal communication from J. F. Bolzan, 1999). The regression line through the data from eight cores gives a surface snow density of 332 kg m−3, increasing to 550 kg m−3 by 15 m depth.
Surface H(t) from Radar Altimetry Data
The observed elevation variations are derived from ERS-1 and -2 radar altimeter data for the period April 1992–April 1999 using procedures described in Reference Zwally, Brenner and FuZwally and Brenner (2001). Elevation changes are derived from surface elevation differences, dH 21 = H 2 − H 1, measured at crossover locations where sub-satellite paths intersect at successive times t 2 and t 1. Sets of N values of (dH 21) i are averaged over selected areas to reduce the error of the mean. Time series of surface elevations, H(t), having sufficient resolution to show seasonal changes are created by the sequence of average crossover differences between the first 90 day interval and each of the successive 90 day intervals, combined with the sequence from the second interval crossed with each successive interval, and so forth for the sequences for the third and greater intervals. Typically, crossovers within a 100 km radius and ±250 m elevation of the central point are included. The data are also corrected for an unexpected inter-satellite bias that was determined by Reference Brenner, Zwally, Cornejo, Saba and Sawaya-LacosteBrenner and others (2000) from analysis of crossover differences acquired during the 12 months of overlapping operation of ERS-1 and -2. At Summit, the bias lowers ERS-2 elevations by 10.6 cm relative to ERS-1. A multivariate linear/sine function is then fitted to the H(t) series giving a linear trend, amplitude of the seasonal cycle in the data, and its phase.
The H(t) elevation series for the summit vicinity for the period April 1992–April 1999 is shown in Figure 1 with the fitted linear trend and seasonal cycle. The regression line shows that the winter snow surface is higher than the summer one, with an average amplitude of 25 cm peak to peak and a minimum in mid-July. The average increase of the surface height is 4.2 cm a−1 during this period. Note that a significant (∼0.5 m) downward trend in the surface elevation occurs between 1992 and 1995, followed by a pronounced increase from 1995 to 1999.
Surface-Elevation and Firn-Densification Model
Surface elevation change
The elevation change of the ice-sheet snow surface is the consequence of several vertical velocity components as illustrated in Figure 2. Snow accumulation increases the surface height at the rate A(t)/ρ 0, where A(t) is the accumulation rate and ρ 0 is firn density at the surface (≈300 kg m−3). Firn compaction, ice flow and surface ablation reduce the surface elevation (H),
The velocity components contributing to surface elevation change at a given location are accumulation, ablation, firn compression, and the vertical motion due to ice flow.
where V fc (z 0, t), V ice and B(t) are the components of vertical velocity of the surface due to firn compaction, ice flow and surface ablation, respectively. All components in Equation (1) are a function of time (t), but V ice changes on time-scales much longer than the seasonal changes of the other components. In dry-snow zones such as the Greenland summit, B(t) may be neglected. On short time-scales (e.g. decadal), V ice is approximately constant. Therefore, variations of the surface height H(t) on seasonal to interannual time-scales are determined by variations in the rates of snow accumulation and densification.
Velocity of firn compression
In steady state, the mass flux, V(z)ρ(z), through any firn depth is a constant that is equal to the mean accumulation rate, 〈A〉 = ρ i V ice, where 〈A〉 is independent of time and ρ i is the density of solid ice (917 kg m−3), so V(z) = (ρ i/ρ(z))V ice. Since the total velocity V(z) is equal to V fc(z) + V ice, the velocity of firn compaction, V fc(z), due to densification of firn below depth z is given by
For example, near the surface, V fc(0) ≈ 2.1 V ice, and the total downward velocity of firn near the surface is V(0) ≈ 3.1 V ice. Therefore, during short intervals when A(t) = 0, the rate of surface lowering is about three times faster than the downward motion of ice at depth below the firn, even though the system is in steady state and dH(t)/dt = 0 on longer time-scales.
In non-steady state, V fc(z, t) depends on the integral of the rate of densification below depth z,
where z i is the depth at which ρ(z) is ρ i.
Firn densification rate
The densification rate at depth in dry firn is determined by the overlying pressure, firn temperature and firn density (cf. Reference BaderBader, 1962; Reference Shapiro, Johnson, Sturm and BlaisdellShapiro and others, 1997). In our model, we use a constitutive equation for densification based on the semi-empirical equation derived by Reference Herron and LangwayHerron and Langway (1980):
where K is a “rate constant”, which is solely dependent on temperature, and A is the mean accumulation rate. The constitutive relation given by Equation (4) accounts for effects of stress and temperature on the densification rate, without attempting to describe the specific processes in the pressure sintering of ice such as grain-boundary sliding (Reference AlleyAlley, 1987a), dislocation creep (Reference Wilkinson and AshbyWilkinson and Ashby, 1975) and boundary and lattice diffusion (Reference CobleCoble, 1970) as used in the treatment by Reference Arthern and WinghamArthern and Wingham (1998).
Equation (4) was based on the suggestion that the proportional change in air space during densification is linearly related to the change in stress due to the weight (p) of overlying snow (Reference RobinRobin, 1958), i.e. dρ/dt ∝ (ρ i − ρ) dρ/dt. The weight change is represented by the accumulation rate (A) which is normally a function of time (t), i.e. dρ/dt = A(t)/ρ i, and K(T) represents the constant of proportionality. The exponent α is an empirical constant introduced by Reference Herron and LangwayHerron and Langway (1980) to account for differing stages of densification or densification mechanisms. Using data collected from 17 sites in Greenland and Antarctica, Reference Herron and LangwayHerron and Langway (1980) derived a value of α ≈ 1 for densities less than 550 kg m−3, and we use α ≈ 1 for all densities.
The temperature dependence of the densification rate has been taken to follow the Arrhenius relation
similar to formulations for the processes of ice creep and grain growth. Although both K 0 and the activation energy E have usually been taken to be constants independent of temperature, Reference BaderBader (1962) found that E varies from 41.8 to 100.4 kJ mol−1 for snow with different densities, with a “reasonable mean” of 58.6 kJ mol−1. Reference Herron and LangwayHerron and Langway (1980) empirically derived average E values of 10.16 and 21.4 kJ mol−1 for densities below and above 550 kg m−3, respectively. Reference AlleyAlley (1987a) obtained a value of 41 kJ mol−1 from modeled density profiles for four Antarctic sites, a value which is rather close to the one often used for grain growth (42 kJ mol−1; Reference PatersonPaterson, 1994). However, previous studies on ice creep and grain growth (e.g. Reference Barnes, Tabor and WalkerBarnes and others, 1971; Reference Budd and JackaBudd and Jacka, 1989) suggested that the activation energy is a function of temperature. Based on laboratory data on grain-growth and ice-creep rates, Reference Jacka and LiJacka and Li (1994) examined the temperature dependence of the activation energy (Fig. 3a) by applying Equation (5) to data for small increments of the temperature. They found that E(T) increased significantly with temperature, especially above −10°C. Values of E for grain growth were similar to the values for ice creep.
Field and laboratory experiments of ice-grain growth data show that both activation energy E (a) and the rate “constant” K 0 (b) are strong functions of temperature. E(T) data and K(T) values obtained from Reference Jacka and LiJacka and Li (1994).
Increasing values of E(T) with increasing temperature cause the rate K(T) in Equation (5) to decrease with temperature, if K 0 is taken to be a constant. However, a decrease would be in marked contradiction to the measurements of K(T) that show a strong increase with temperature (Reference Jacka and LiJacka and Li, 1994). To increase with temperature strongly enough given the data on activation energy, the factor K 0(T) in Equation (5) must increase with temperature. Using the K(T) data for grain growth given by Reference Jacka and LiJacka and Li (1994, table 2) and their derived E(T) values, we derived the corresponding K 0G(T) as shown in Figure 3b. The best-fit power-law curves through the activation energy E(T) and the corresponding values of the K 0G(T) factor are shown in Figure 3a and b. Both curves have correlation coefficients greater than 0.99. The stronger temperature dependence of K 0G(T) compared to E(T) is shown by the larger exponent on the inverse temperature.
To account for differences between the rates for the processes of densification and grain growth, we introduce an empirical constant β in the rate factor K 0(T):
where K 0(T) is the rate factor for densification and K 0G(T) is the derived factor for grain growth. Changing the value of β varies the rate of densification, and should therefore change the gradient of density with depth in the modeled density profiles. Figure 4 shows a numerical experiment with our densification model (described below) for three values of β using the smoothed density profile from the Greenland summit as an initial condition in the model. As expected, the gradient of the density changes with β. In addition, the modeled amplitude of the seasonal cycle of density also varies with β. The chosen value of β = 8.0 produces a modeled density profile that matches the smoothed density profile of the initial condition (also shown in Fig. 4). Although the measured profile appears to lie notably below the modeled profile in the figure, the smoothed value is actually very close to the measured density (difference is <10 kg m−3), because the modeled seasonal variation in density is not symmetrical.
Comparison of density profile computed at three different values of empirical constant β as indicated beside each curve.
Firn temperature profile
Li and others (in press) have analyzed firn- and ice-temperature evolution driven by surface air temperature during the period 1982–99 for the Greenland summit. The temperature analysis follows the standard one-dimensional time-dependent heat-transfer equation (Reference PatersonPaterson, 1994) written by
where c = c(T) is the heat capacity and k = k(ρ(z), T) is the thermal conductivity, v is vertical velocity, and f is internal heating. For the advection term, which is small, we use v(z) = 〈A〉/ρ(z) based on steady-state conditions and a constant 〈A〉 = 250 kg m−2 a−1. The small f term is also neglected.
Figure 5a presents computed firn-temperature variations at several selected depths over the period April 1992–April 1999. We use daily mean temperature data (Fig. 1) and an approximation for the diurnal temperature cycle of about 9°C. The approximation in our model takes the temperature for half of each day to be the daily mean plus 4°C, and for the other half-day the daily mean minus 4°C. The temperature analysis shows that in addition to the seasonal cycles of firn temperature, significant interannual variations of the mean summer (June–August) temperatures occur at various depths. The interannual variability of summer temperature is > 5°C at the surface, diminishing to about 1°C at 10–15 m, as shown in Figure 5b. The trend in surface summer temperature is positive at 0.3°C a–1, and the trend in mean annual temperature for this period is also positive at 0.4°C a−1.
(a) Variations of firn temperature calculated at several depths above 15 m for the Greenland summit region, April 1992–April 1999. The number beside each curve indicates the corresponding depth. (b) Interannual variations of summer (June–August) mean firn temperatures at several depths, April 1992–April 1999. The trend in surface summer temperature is +0.3°C a−1.
Multi-layer numerical model
The surface height variations with time (Equation (1)) depend on the variations of new accumulation at the surface, A(t), and variations in the velocity of firn compaction, V fc(z, t), as given by Equations (3–6). The rate of firn densification at time t depends on the existing density profile ρ(z, t), as well as the temperature profile, T(z, t), which also depends on the density profile according to Equation (7). Therefore, V fc(z, t) depends on the densification and temperature history, and solution of the time-dependent equations must be coupled in order to calculate V fc (0, t) and dH/dt.
The equations are solved using a multi-layer numerical model. Taking the density within each layer n of thickness hn to be uniform, the rate of change (vnj ) of the layer thickness that is caused by the density change at depth zn for one time interval (dt) at time j is given by
We consider densification in the firn between the surface and depth z d to be the seasonal part, where z d ≥ 15 m, because of seasonal variations of firn temperature in this range, and below z d as the non-seasonal part. Below z d, the vertical velocity of firn compression is calculated based on Equation (2) under the steady-state conditions from
which is therefore independent of the small temperature variations below z d. Above z d, Equation (7) is used to calculate firn temperatures.
After substituting V ice = 〈A〉/ρ i, neglecting B(t), and splitting V fc(z) into three terms, Equation (1) for the height change dH for each time-step, dt, becomes
The model is initialized at t = 0 with the measured density profile measured at the summit (personal communication from J. F. Bolzan, 1999). The vertical velocity from additional accumulation at each time-step is upward by A(t j)dt/ρ 0, followed by a decrease on subsequent days as the new precipitation compresses. The compression of the new layers is given by the second term, which sums from 1 to M − 1 the change in thickness of these layers, as given by Equation (8), where M is the new surface layer formed at t j, and M = t j/dt is the total number of these layers. The firn from the initial surface to z d, which is initially set to 15 m, is divided into N layers, each with initial thickness of 0.1 m. Near the surface, these layers correspond to about 1/7 of the annual deposition. The third term sums the change in thickness of the N layers, also as given by Equation (8). The fourth term is V fc(z d) for the compression below z d, which is given by Equation (9). The surface elevation H(t) is obtained by numerical integration of Equation (10).
Results and Discussion
Seasonal elevation and densification variations
To examine seasonal variations, the densification/elevation model is driven by a constant accumulation rate of 〈A〉 = 250 kg m−2 a−1 and a fixed seasonal cycle of surface air temperature. The temperature cycle has an annual mean of −29°C and a seasonal amplitude of 27°C, with a maximum In July. The model time-steps are dt = 12 hours, and for each time-step A(tj ) · dt = 〈A〉/730. Diurnal temperature variations, which can be significant in the top meter of firn, are included in the model because of the strong non-linear dependence on temperature. As noted above, the diurnal cycle is approximated by adding 4°C to the seasonal cycle for half of each day and subtracting 4°C for the other half. This treatment is rather close to field observations. Reference Alley, Saltzman, Cuffey and FitzpatrickAlley and others (1990) monitored hourly surface snow-temperature variation at the GISP2 site. Their data indicate that 12 hour mean surface temperature is about 5°C above and below the daily mean.
Figure 6a shows the seasonal variations in the surface elevation, which are expanded in Figure 6b, as well as the sinking of the initial surface as it is buried with additional accumulation. April is chosen as the starting month because that is the start of the altimeter elevation series and corresponding data in Figure 1. In this case, the driving of surface increases is constant in time, and the resulting seasonal variations in the surface elevation are driven solely by the temperature-dependent variations in densification. The modeled peak-to-peak amplitude of the seasonal variation of surface height is 18 cm, which is close to the 25 cm mean amplitude shown in the radar altimetry measurements. The minima occur in late August which is also close to the mid-July average minima in the altimeter elevation data.
Modeled steady-state seasonal variations of surface elevation for the Greenland summit (a, b) together with corresponding density changes (c). Solid and dashed lines in (c) show the modeled density variation and thefield density measurements (personal communication from J.F. Bolzan, 1999), respectively. The variations are driven by surface air temperature with a regular annual cycle with amplitude 27° and mean ∼29°C at constant accumulation rate of 250 kg m−2 a−1 The seasonal variation in surface elevation is 18 cm, but the seasonal variation in the height of the initial surface decreased quickly as it was buried by new accumulation. The seasonal variation in density is about 140 kg m−3, which decreases with depth.
Much of the firn densification and consequent surface lowering occurs during about 3 months of the late-spring to early-summer season when the upper firn is warmest. After the maximum in summer temperature in July, the upper firn cools and the rate of densification decreases. After the period of maximum compaction ends in late summer, the accumulation rate dominates the compaction rate through winter until the next spring and the surface rises. The enhanced firn compaction driven by the higher temperatures of summer is most striking in the near-surface firn within 1 m, and quickly diminishes with increasing depth, as indicated by the small variations in elevation of the initial surface curve in Figure 6a.
The diurnal temperature cycle is also significant in warming the upper part of the firn and increasing the rate of summer densification. When the model is run with a fixed daily temperature, without the approximation for the diurnal cycle, a slightly large value of β = 11 is required to match the modeled density profile to the measured profile. The corresponding seasonal amplitude is also somewhat smaller at 11 cm instead of 18 cm.
The seasonality of the compaction rate is also reflected in the seasonal variations in the density profile (Fig. 6c). Corresponding density variations with depth are shown in the density measurements made at 10 cm depth intervals. The density-profile data from the summit core (personal communication from J. F. Bolzan, 1999) for the top 4 m are also shown in Figure 6c. The model and the core data are reasonably similar considering that the point measurements of the core also include spatial variability in deposition processes not included in the model. The magnitude of the density variations, which is about 100 kg m−3 in amplitude, reduces with increasing depth in both the model and the measurements, but not as much in the field data below 4 m as would be implied by the small amplitudes at 2–4 m depth. Below 4 m, the core data also show approximately one cycle in density per year. A similar seasonal density cycle has been observed in high-resolution density measurements made on firn/ice cores in Antarctica (Reference Gerland, Oerter, Kipfstuhl, Wilhelms, Miller and MinersGerland and others, 1999), and may be reflected in the typically wide scatter of densities in other less detailed density profiles as noted in the introduction.
Interannual elevation variations
As shown in Figures 1 and 5, both the temperature and accumulation datasets used to drive the densification/elevation model have significant interannual variations during the period April 1992–April 1999. The modeled H(t) elevation series from driving the model with the temperature and accumulation data is shown in Figure 7. For comparison, the observed H(t) elevation series derived from radar altimetry measurements (overlain in Fig. 7) closely follows the modeled elevation, although the observed series shows several larger spikes. Both curves indicate that surface elevation at the Greenland summit decreased about 40–60 cm between 1992 and 1995 followed by a marked increase of 50–80 cm. Along with these variations, both curves also show significant seasonal peaks, caused by the seasonal cycles of the surface temperature in the model.
Comparison of modeled and observed interannual changes of the snow surface elevation, 1992–99, at the Greenland summit.
To separate the effects of temperature and accumulation variations that are input to the model, we first compute the elevation change with the observed temperature variations (seasonal and interannual), but with a fixed accumulation rate (curve a in Fig. 8). Curve a shows the seasonal cycle in surface elevation with a sharp lowering during late spring to early summer. However, it also shows an interannual variability in both the seasonal amplitude and the year-to-year changes. For example, the seasonal amplitude is large during the warmer summer of 1995 (cf. temperature in Fig. 5b) and small during the colder summer of 1996. There is also a general decrease in elevation, which is due to the increasing trend in temperatures during this period, particularly in summer. A modeled elevation decrease of about 20 cm is the response to the increase in summer temperatures during this period.
Modeled elevation change caused by (a) observed temperature variations with constant accumulation and (b) observed/modeled accumulation variations with a seasonal variation in temperature without interannual variations.
To examine the effect of accumulation variations (interannual and seasonal), we compute the elevation change with the observed variable accumulation rate and a fixed seasonal temperature cycle, but without any interannual variations in temperature (curve b in Fig. 8). Analysis of the accumulation data in Figure 1 using the multivariate fit (linear with a sinusoidal seasonal cycle) gives an increasing trend of 10 kg m−2 a−1 a−1, a seasonal amplitude of 119 kg m−2 a−1, and a seasonal maximum at the end of May. For the first 3 years, in particular, the seasonal cycle of accumulation is out of phase with the cycle of firn compaction, which is shown by the reduced amplitude of the seasonal cycle of curve b vs curve a.
Curve b shows that the interannual variations of accumulation dominate the interannual changes of the surface elevation over this period. Although the 7 year trend in accumulation is positive, there is a significant decrease through mid-1995 followed by a marked increase. In response to accumulation changes, the surface elevation decreases notably between 1992 and 1995, followed by a significant increase. Although the decrease from 1992 to 1995 is primarily driven by accumulation, the warmer summer temperature of 1995 also contributes to the minimum shown in both the model and the altimeter observations. Over the 7 years, the effect of warmer temperatures lowers the surface by about 3 cm a−1, but the effect of increased precipitation dominates to bring the overall increase to 4.2 cm a−1 as obtained from the multivariate fit to the altimeter data. Supporting information for the variation in accumulation rates is provided by the analysis of shallow ice cores, which shows an anomalously low accumulation in 1995 followed by marked increase in 1996 (Reference McConnell, Mosley-Thompson, Bromwich, Bales and KyneMcConnell and others, 2000).
Summary and Conclusions
A principal feature of our densification model is the stronger temperature dependence of the constitutive relation governing firn densification. This stronger temperature dependence is based on laboratory measurements of rates of grain growth and ice creep. Firn densification, grain growth and ice creep are processes that have been previously treated similarly to rates described by the Arrhenius equation, but with activation energies (E) and rate constants (K 0) that were independent of temperature. We apply the temperature dependence of E and K 0 found in the laboratory measurements of grain growth to the densification process as well. Although previous treatments took E and K to be independent, they usually used different values of E and K for different processes or stages of processes, such as the densification in the upper and lower levels of firn. Our treatment uses a multiplicative factor β with K 0 to account for differences between the processes of densification and grain growth. The value of β is empirically selected to match the modeled density profile to measurements over the depth range 0–40 m. A key feature of this result is the matching of the density profile with one set of parameters. Secondly, after matching the density profile, the same set of parameters produces the seasonal cycle in compaction and elevation change.
A principal result of the model is that most of the firn compaction in the upper firn occurs in late spring and early summer, when temperatures are warmest. This seasonal variability in compaction rate is sufficient to cause a cycle in the density profile that is gradually dampened in amplitude as the firn continues to compress with depth. This provides an explanation of observed density variations in firn, which have previously been attributed mostly to depositional processes rather than to post-depositional variations in compaction. The model shows that snow deposited from winter to summer is compacted to the higher densities, and snow deposited from summer to winter is compacted to the lower densities. This result is in contrast to depositional density associations that tend to be characterized in the literature as due to either summer or winter processes.
The agreement between the modeled surface elevation, which is driven by observed temperatures and accumulation variations, and the elevation time series derived from radar altimetry provides an explanation for the seasonal and interannual elevation variations observed in the altimetry data. The modeled seasonal variations are largely due to variations in compaction, but are also influenced by seasonal variations in precipitation, particularly from precipitation events associated with large storms in a few years. Although the agreement between the modeled and observed amplitude and phase of the seasonal cycle is good, there are small residual differences. Therefore, other possible contributions such as variations in radar backscatter and penetration depth may also have some secondary effect on the altimeter results.
The interannual variations produced by the model are driven in part by interannual variations in temperature and compaction, but are mainly caused by the interannual variations in precipitation. In particular, although the marked decrease in elevation in summer 1995, followed by a marked increase, is partially driven by warmer temperatures in 1995, the predominant effect is the interannual variation in precipitation. The precipitation variation is supported by independent evidence on accumulation from ice cores. The observed elevation trend for 1992–99 is 4.2 cm a−1, which is slightly larger than the 3.6 ± 2.1 cm a obtained by Reference Hamilton and WhillansHamilton and Whillans (2000) by the global positioning system (GPS) coffee-can method at the GISP2 summit and different from the −3 ± 4 cm a obtained by Reference Hvidberg, Keller, Gundestrup, Tscherning and ForsbergHvidberg and others (1997) nearer the actual summit about 25 km to the east.
At any time, the rate of densification depends on the temperature and accumulation history, which was simplified in our present model by assuming the firn was in steady state and using a measured mean density profile on the initial day (14 April 1992). However, longer-term temperature and precipitation records can be used to calculate an initial state of temperature and density profile with seasonal variations. Future efforts should also include laboratory measurements of the temperature dependence of firn-densification rates to refine the coefficients used in the rate equation, direct measurements of densification rates vs depth by automated instrumentation placed in the firn, and extension of the elevation model to firn percolation zones and cold, low-accumulation regions. Finally, other forcings such as drift and wind compaction to the snow may need to be considered. Strong wind may enhance snow compaction producing higher-density firn layers (Reference GoodwinGoodwin, 1991), while in other cases it also brings in large amounts of low-density snowfall in very short periods causing larger fluctuations of surface height. Some of the departures between modeled and observed results shown in Figure 7 may be due to such causes.
The effect of interannual variations in densification (caused by changes in temperature, density and accumulation rate at the surface) on the use of altimeter observations for studies of mass balance was addressed by Reference Arthern and WinghamArthern and Wingham (1998) and Reference WinghamWingham (2000). While a greater sensitivity to temperature would seem to complicate the interpretation of satellite altimeter data, it may only change the response time to variations in accumulation. Since the response time of the surface elevation to accumulation anomalies is governed in effect by how fast the anomaly moves down the density profile, on average the greater sensitivity to temperature should not make much difference. However, accumulation anomalies that occur during the warmer seasons when near-surface densification is faster should persist for a shorter time than anomalies occurring during the colder seasons. Satellite observations of surface temperature (Reference ComisoComiso, 2000), calibrated with data from AWS stations on the ice sheets, can be used to model the elevation changes expected from temperature variations, so that elevation changes caused by precipitation variations and ice flow can be evaluated in satellite altimeter data.
Acknowledgements
We appreciate very much discussions with R. B. Alley and M. K. Spencer on their work on densification, the assistance of W. Abdalati and K. Steffen in providing and helping us to interpret the AWS data, the cooperation of D. H. Bromwich in providing the precipitation modeling results, the help of W. L. Wang in the analysis of the temperature profiles, and the suggestions of reviewer R. B. Alley.
