Introduction
During the 1987 Greenland field season we recovered nine shallow cores to depths of approximately 17 m to determine the spatial variation in the accumulation rate around the Summit site in central Greenland (Fig. 1). This study was part of the GISP2 site-selection activities, which also included an airborne-radar survey to determine the surface and basal topography (Reference Hodge, Wright, Bradley, Jacobel, Skou and VaughnHodge and others, 1990) and the measurement of surface velocities using Doppler satellite-surveying techniques (paper in preparation by J.F. Bolzan).
We established a 150 km × 150 km survey grid centered on the Summit site by a series of surface traverses. Figure 2 shows the location of the grid survey sites and the drill sites in relation to the surface topography. The coordinates of the drill sites are given in Table 1. 2 m pits were excavated at each drill site and one pit wall was sampled every 10 cm for stable oxygen-isotope and gross ß-activity measurements. The shallow cores were recovered by drilling from the floor of each pit using a lightweight hand-auger with a 2 m barrel fabricated by the Polar Ice Coring Office (PICO). The drill was first used at site 57 in the southeast quadrant and was irretrievably stuck at about 10 m depth. Minor modifications to a spare drill barrel subsequently resulted in excellent core recovery, typically about 97%, at the other drill sites, with the cores being about 15 m long. The nine cores were returned frozen to the Byrd Polar Research Center for sampling and analysis. The visual stratigraphy was examined on a light table and mapped.Strata identified as fine-grained and coarse-grained showed no correlation with the δ 18O value. Also, no melt features were seen in any of the cores.
Method
The average accumulation rate,
over the time interval above a specified reference horizon, is
where T is the number of years since the deposition of the horizon, and L(T) is the ice-equivalent length of the core above the horizon. If the core is composed of N(T) segments above the reference horizon, the ice-equivalent length is
where ρk , lk , mk and rk are the density, true length, mass and radius of the kth segment respectively, and ρ ice is the density of ice. Complications arise because of core loss and because of the difficulty in accurately measuring the dimensions of irregularly shaped core segments. Ambiguities in determining the depth of the dated reference horizon also result in uncertainty in the ice-equivalent core length L(T). For this study, the dated reference horizons were based on stable oxygen-isotope measurements which generally showed a well-defined annual signal. Gross ß-activity profiles were also measured to aid in resolving dating ambiguities.
Core loss
To estimate the contribution of core loss to L(T), we used the algorithm developed by Reference Whillans and BolzanWhillans and Bolzan (1988) which determines the core-loss zones by identifying the drill runs where the core break is likely to have occurred at the bottom of the hole. The total core loss between these special runs can then be arbitrarily distributed at the top and the bottom of the corresponding depth interval without violating the requirement that the base of the recovered core be less than or equal to the drill depth. Here, the total core loss between adjacent special runs was distributed equally above and below the total core recovered between the runs. The density for the core-loss region was assumed to be the average of the densities of the core segments immediately above and below. However, because the core loss here amounts to only a few per cent of the total core recovered, our results were not significantly affected by the way in which core loss is distributed between the special runs.
Radius and mass measurements
Measurements of the radius and mass of the core segments were made in the cold laboratory at the Byrd Polar Research Center. The radius of a fim core can vary significantly from segment to segment, and occasionally within a segment, due to wobbling of the drill and abrupt changes in density. Also, handling and shipping causes grains to be shed off the most friable segments, and increases the irregularity in the shapes of these pieces. The error analysis presented in the Appendix shows that the uncertainty in the calculated accumulation rate is by far most sensitive to the uncertainty in the radius of the core segments. Using a vernier caliper with a precision of better than 0.1 mm, we measured the diameter of each core segment at several positions to determine the average radius, and took the uncertainty for that segment to be one-half of the difference of the maximum and minimum radius values. For the lowest-density segments, variations in the measured radius of 1 mm were not uncommon. Deeper in the core, diameter measurements showed little variation and we conservatively took ± 0.5 mm to be the minimum uncertainty in the radius of a core segment. The mass was measured using a triple-beam balance and we again conservatively estimated the uncertainty to be ± 1 g. As a result, even for the segments deep in the core, the contribution to the uncertainty in the ice-equivalent length from the radius errors is about an order of magnitude greater than from the mass uncertainty.
Gross ß-activity measurements
Atmospheric testing of thermonuclear bombs in the 1950s and 1960s resulted in elevated β activity in snow deposited during this time, due primarily to increased concentrations of 90Sr and 137Cs. The first thermonuclear tests in 1952 and 1953 produced a two-fold increase in β activity over the natural background level of ∼ 100 dph kg−1 for snow deposited in late 1953. The Castle series of tests in early 1954 resulted in an activity increase of 5–10 times above the pre-testing background for snow deposited from the summer of 1954 to early 1955. A minimum in β activity apparent in the 1960 layer was due to the testing moratorium in 1959, and the resumption of Soviet tests in 1961 and 1962 caused very high activities in the 1963 layer. This 1963 horizon is very useful for accumulation-rate determinations, as it can be readily identified, even in the absence of other corroborating stratigraphic data (for further details on the variation of gross β activity in polar firn, see Reference Picciotto and WilgainPicciotto and Wilgain (1963); Reference Crozaz, Langway and PicciottoCrozaz and others (1966); Reference Hammer, Clausen,, Dansgaard, Gundestrup, Johnsen and ReehHammer and others (1978); Reference Koide, Goldberg, Langway, Oeschger and DansgaardKoide and Goldberg (1985)).
Samples 10 cm long were prepared for gross ß-activity measurements by pumping the meltwater through ion-exchange filters that collected the strontium and caesium cations. Rather than sampling over an entire core, samples were taken from a zone that bracketed the estimated depth of the 1963 high-activity horizon. Gross β activity was measured on the dried filters using either a Canberra model 2200 or a Tennelec model 1000 gas-flow proportional counter, both with a counting efficiency of about 30%. During the 1 year period over which these measurements were made, the background activity for the Canberra and the Tennelec averaged about 0.65 and 0.45 counts min−1, respectively. However, both counters are located well above ground and are susceptible to significant short-term fluctuations in the counting background. Background measurements were taken before and after counting all of the samples from a given core, mainly to check for instrumental malfunctions, and we made no attempt to determine the short-term background fluctuations that might vary from sample to sample. Thus, in calculating the activity of each sample, we used the average background counted over the duration of all the core measurements. While our β-activity results cannot accurately resolve small differences in the true activity between samples, they are adequate to detect the large changes in activity seen around the most prominent ß-activity horizons.
Oxygen-isotope profiles
Upon return to the Byrd Polar Research Center, the cores were sampled for oxygen-isotope measurements, which were made at the Geophysical Isotope Laboratory at the University of Copenhagen. Samples were 5 cm in length, except where sample size had to be increased or reduced slightly to avoid overlapping a core-loss zone. Measurements of the oxygen-isotope ratio were also made on 10 cm samples from the drill-pit walls. Figures 3 and 4 show the oxygen-isotope values as a function of ice-equivalent depth, along with the the gross ß-activity profiles. The discontinuous lines in the oxygen-isotope profiles correspond to regions of core loss. In all cores, except site 15 and site 57, the 1963 β horizon is readily identified, as is the 1960 minimum; and the 1963 and 1960 summer peaks can be readily assigned. For site 15, the actual accumulation rate was much lower than estimated and, as a result, the β activity was measured below both horizons. However, in this core we see a sharp rise in activity which we attribute to the 1954 Castle tests. For site 57, the drill was lost before the 1963 horizon was reached.
We have dated each isotopic summer peak by taking the first peak below the surface to be 1986, and have used the results of the gross ß-activity measurements to resolve ambiguities. There is Uttle ambiguity in the assignment of years for sites 31, 51, 571 and Summit (see (Fig. 3)), as the number of isotopic peaks in these profiles is consistent with the number of years between the surface and the well-defined 1963 layer. For site 571, we are forced to take two poorly defined peaks assigned to 1967 and 1971, but there are no other candidates for these years. For site 57 (Fig. 3), the summer peaks are well defined and the dating is relatively unambiguous, even without the 1963 β horizon. However, the profiles from sites 13, 15, 37 and 73 show some inconsistency between the gross β and oxygen-isotope stratigraphy (see (Fig. 4)) in that there are too many or too few features that could be identified as peaks between the 1986 and 1963 layers.
Site 73 has the highest accumulation rate of the nine sites, and below the 1963 horizon the isotopic peaks are well-defined. However, above this horizon there are three too many features which could be considered summer peaks. We choose to ignore the features labeled with a “?” on the site 73 profile in Figure 4, as they are two-sample or three-sample peaks which are much less well-defined than the other peaks in the profile.
The site 13 profile contains three extra features which are candidates for summer peaks. The 1984 peak seems to appear both at the bottom of the pit and the top of the core, which may be due to the pit floor at the drill hole being higher than the base of the pit wall where the oxygen-isotope samples were taken. We choose to ignore the two-sample peak seen in 1962–63 and a small bump in 1966–67, which is much less prominent than its neighbours.
The site 37 β profile is consistent with the oxygen-isotope profile below the 1963 horizon. Above this horizon, the isotope profile seems to be missing one well-defined peak. We have assigned 1982 to a single-sample peak lying in the shoulder of the 1981 summer peak. An alternative would be to assign 1972 to a single-sample peak now between 1971 and 1972, but there is some support for our assignment in that high accumulation in 1971–72 is seen also at sites 57 and 13, while no anomalously high accumulation is seen at any sites in 1982–83.
For site 15, if the sharp rise in β activity at 600 cm depth is due to the Ceistle series, then the isotopic profile is missing one peak. We attribute this to core loss, as the core-loss algorithm shows that enough loss could be assigned at about 200 cm depth to account for a five-sample isotopic peak. There is no ambiguity in the dating above this zone, so that the missing peak is assigned to 1976.
There is also some ambiguity as to the time of the year over which an isotopic “summer” peak may occur. Generally, snow deposited from May through August has the least negative δ18O values, and we assign an uncertainty of ± 2 months to the date of the peak’s center within a calendar year.
Results
The average accumulation rate at each of the core sites for the period 1959−87 is given in Table 1, except for site 57 where the record begins at 1964. The uncertainty includes the contribution from core loss, radius and mass measurements, and dating. Using these results and the average accumulation rates from 1943 to 1973 reported by Reference Clausen, Gundestrup, Johnsen, Bindschadler and ZwallyClausen and others (1988), we have plotted constant accumulation-rate contours (isohyets) over the 300 km × 180 km area extending from below Crete to just above the northernmost sites in the grid (Fig. 5). While the values derived here are for a more recent time interval, comparing average accumulation rates over the first 30 years from sites 13, 15, 37 and 571, suggests that the more recent averages are less than the 1943–73 values by about 1–2 cm year−1, which does not differ significantly from our measurement uncertainty.
The isohyets have a spatial gradient that is relatively uniform and a persistent northwest-to-southeast direction that shows little influence due to the surface topography. A similar pattern in surface accumulation for this region was inferred by Reference MockMock (1967) on the basis of a trend-surface analysis, with the multiple-regression coefficients determined from snow-pit studies at 127 sites (including 35 values dating from the work of Koch and Wegener in 1912). Accumulation rates calculated at the grid-coring sites using Mock’s regression model (1) for northern Greenland agree with the values reported here within ±3 cm year−−1 for all sites except 571, for which the predicted value is too large by about 4 cm year−1. Given that the regression coefficients were determined from snow-pit studies that often only spanned a few years at a given site, and from studies conducted more than 50 years apart, the agreement between modeled and measured values is remarkable. We conclude that the large-scale pattern in accumulation rates at least over the northern half of the ice sheet has probably varied on the order of only about 10% since the beginning of the century, which is consistent with the accumulation-rate variations shown in Reference Clausen, Gundestrup, Johnsen, Bindschadler and ZwallyClausen and others (1988) based on shallow ice-core oxygen-isotope profiles.
The northwest-to-southeast trend in the isohyets suggests that the spatial pattern of accumulation is due to air masses that ascend the ice-sheet slope from the southwest (Reference Barry, Kiladis, Radok, Barry, Jenssen, Keen and KiladisBarry and Kiladis, 1982), even though the grid area is nearly equidistant from both the east and west coasts. This is supported by an analysis by Reference KeenKeen (1984), based on the mean vorticity-flux index (VFI) for Greenland from 1946 to 1980. The VFI is a vector quantity defined as the product of the positive vorticity at the 500 m bar level (a scalar) and the long-wave wind field (a vector), and thus is a measure of the advection of cyclonic activity. The VFI is also roughly proportional to the average horizontal flow direction during cyclonic events, which often result in precipitation on the ice sheet. Keen found that upper-level systems approach Greenland from the west and the southwest, and his results show the mean VFI in central Greenland to be nearly perpendicular to the isohyets in Figure 5.
Since there is no longer a moisture source once an air mass ascends the ice-sheet slope, the water vapor should become progressively more depleted in the heavier isotopes as precipitation continues to occur. For moisture flux primarily from the southwest, the isotopic values of the snow should become progressively more negative toward the northeast. Post-depositional processes such as drifting and sublimation alter the isotopic composition of the surface snow, and generate stochastic, high-frequency, isotopic noise. Much of this noise can be eliminated by averaging over a sufficiently long time. In Figure 6, we have plotted contours of average isotopic value based on averages over the 1964–87 interval common to all sides. The isopleths trend nearly in the same direction as the contours of constant accumulation (Fig. 5), consistent with the progressive depletion of the precipitating air masses as they track from the southwest. Note, however, that the average isotopic values for sites to the east of the topographic crest are somewhat enriched, perhaps because these sites receive more moisture from air masses ascending the eastern slope of the ice sheet.
More directly, we can estimate the magnitude of the moisture flux as a function of direction by using automatic weather station (AWS) data. These data have been recorded in central Greenland since 1987 and include air temperature, pressure, wind speed and direction, and relative humiditiy reported at 3 h intervals (Reference Stearns, Weidner, Weiler, Wilson and SeverinStearns and Weidner, 1991). The longest data sets (∼2 years) currently available are from the Cathy (Summit) and Kenton sites. The Cathy AWS operated at Summit from May 1987 through May 1989 and then was re-installed at Kenton in June 1989. Unfortunately, for both sites, there are almost no wind data available between January and May of each year, probably due to the formation of rime ice on the aerovane (Reference Weidner, Stearns, Weiler, Wilson and SeverinWeidner and Stearns, 1991), so that approximately 40% of the wind data are missing for each site. While the missing data are over the coldest time of the year, when the moisture flux may be expected to be small, temperatures in winter and early spring are highly variable with a number of warm events (Reference Weidner, Stearns, Weiler, Wilson and SeverinWeidner and Stearns, 1991), for which the flux could be significant.
It is not the magnitude of the moisture flux per se that is related to the precipitation over a region, but rather the vertically integrated flux divergence. However, here we make the simple assumption that precipitation is more likely when the moisture flux is large. The AWS data enable us to compute only the surface horizontal moisture flux, but as the water-vapor content of the atmosphere decreases rapidly with height (Reference Peixóto, Oort, Street-Perrott1, Beran and RatcliffePeixóto and Oort, 1983), the bulk of the vapor transport is near the surface. Thus, we also assume that the important variations in the vertically integrated horizontal vapor flux are mirrored by the variations in the horizontal surface flux. The latter is calculated for each usable AWS observation as the product of the wind velocity and water-vapor density, ρ v, which is given by:
where e s, is the saturation vapor pressure at temperature T (°K), f is the relative humidity and Rv = 461 J kg−1 K−1 is the gas constant for water vapor (Reference RogersRogers, 1979). Equation (3) follows from the ideal gas law by taking the relative humidity as proportional to the ratio of the vapor pressure to the saturation vapor pressure, an approximation that is better than 1% at the low temperatures measured here. The saturation vapor pressure is computed from the Goff-Gratch formula as given on p. 350 in the Smithsonian meteorological tables (Reference ListList, 1971).
We expect the net contribution to precipitation to be related to the moisture flux and the number of occurrences of wind from a given direction, and in Figure 7 for both the Cathy and Kenton data we have plotted the number of observations and the average horizontal surface-moisture flux as a function of wind direction for 10° angular sectors. Figure 7a shows that the winds most frequently tend to come from the southwest quadrant at Cathy, with 42% of all observations being between 200° and 290°. For the Kenton AWS, which is approximately 40 km west of the topographic crest, Figure 7a shows a broad peak in the wind distribution centered at about 150°, with 58% of all observations being between 110° and 210°. This shift may be due to a change in the wind field with time, and possibly a katabatic effect, which would skew the distribution to the east. Figure 7b shows that at both sides the average flux is largest from the southwest quadrant. This is because temperature, wind speed and relative humidity all generally tend to be higher for winds from these directions. Note that the Cathy data show a secondary peak between 160° and 190° which, however, overlaps a local minimum in the wind frequency. For Kenton, the vapor-flux peak in the southwest quadrant barely overlaps the peak in the wind distribution and, while there is likely to be a significant contribution to the net precipitation from this direction, it is unclear whether the data indicate the primary moisture source is from the southwest. Whether this reflects a secular change between the two data sets, or the effect of a katabatic component at the Kenton site is also uncertain. However, for Cathy we conclude that the AWS data seem to be consistent with the primary moisture flux being from the southwest, in that the average vapor flux and the wind distribution are both peaked in the southwest quadrant.
The question remains as to the initial source areas of the water vapor. Reference Johnsen, Dansgaard and WhiteJohnsen and others (1989) concluded that the high deuterium excess measured at site G was consistent only with a source region with high ocean-surface temperature and mixing ratio such as the subtropical North Atlantic at about 35° N. Their deuterium-excess model also suggested that coastal sites in Greenland receive a significant amount of precipitation from sources with low surface temperatures and mixing ratios, such as Baffin Bay and Denmark Strait. According to Reference Johnsen, Dansgaard and WhiteJohnsen and others (1989), cyclones that result in precipitation on the ice sheet form along the polar front in the North Atlantic south and east of Newfoundland, where the prevailing winds are from the south and west. As the cyclones move northward they become occluded, with the vapor trapped above the occlusion and isolated from the colder ocean surface. This scenario is compatible with a study by Reference Barry, Kiladis, Radok, Barry, Jenssen, Keen and KiladisBarry and Kiladis (1982) on the origin of cyclones affecting Greenland during 1950–65 (i.e. those whose centers were within 200 km of the Greenland coast). They found that about half of the average number of cylcones per year tracking through Baffin Bay and Denmark Strait originated south of 60° N, while only about 30% were of local origin. Barry and Kiladis also noted that accumulation time series from snow-pit data generally show a moderate degree of spatial coherence over a few hundred kilometers, but correlate poorly with coastal precipitation records. While part of the poor correlation may result from the difficulty in measuring snow accumulation at coastal sites, they also concluded that coastal precipitation is strongly influenced by local effects. If the primary moisture source for inland sites is assumed to be the subtropical North Atlantic, then the cyclones originating there and tracking near the Greenland coast must tend to be deeper and more extensive than those that form locally in Baffin Bay.
The well-defined annual signal in the oxygen isotopes enables accumulation-rate time series to be constructed for each core. Using Equations (1) and (2), we have computed the ice-equivalent depth for each isotopic summer peak based on the dating shown in Figures 3 and 4; the results are given in Table 2. In order to identify the recent trends, for each of the cores we have plotted 5 year averages as a function of time in Figure 8. While few of the individual 5 year averages deviate from the mean for the site by more than the measurement error (which includes a ±2 month uncertainty for the summer peak at the beginning and end of each 5 year interval), every site shows an increase for the interval 1982–87 over the previous 5 year period. Also, except for sites 15 and 51, the average over the 5 year interval 1977–82 is less than the previous interval. However, except for this recent feature, the data do not show any persistent longer-term trends.
The lack of a persistent trend in the individual accumulation-rate time series implies that there has been no recent trend in the total annual accumulation over the study region. To check this, we computed the total accumulation for each year during 1964–86 over the 180 km × 180 km region shown in Figure 2. We used a gridding algorithm to interpolate the accumulation rates for the nine sites over a regular grid centered on Summit, and then integrated the interpolated values. Here, the annual accumulation rates are referred to accumulation years, spanned by the isotopic peaks. The results are shown in Figure 9. Because the total accumulation for each year is based on only nine measured values, the detailed structure displayed here in the total-accumulation time series may be questionable. However, it is apparent that there is significant short-term variability, with the standard deviation being about 11 % of the mean value. This amounts to a significant reduction in the variance of the individual accumulation-rate time series, for which the standard deviation varies from about 15% of the mean (sites 31, 51, 73 and Summit) to 20–25% of the mean (sites 13, 15, 37,57 and 571). A linear regression gives an extremely small positive slope (0.01 km year) that is not different from zero given the likely uncertainty (∼1 km3 year−1) in the annual values, and indicates no trend over the period 1964–86 in the net annual accumulation in the central Greenland survey area.
Summary
The spatial variation in average accumulation rates determined from the oxygen-isotope and gross ß-activity stratigraphy in nine firn cores suggests a moisture flux source to the southwest of the central Greenland ice sheet. Progressive depletion in the average stable oxygen-isotope values along the southwest to the northeast direction is consistent with this picture. Automatic weather-station data over 2 years from the Summit site near the topographic crest indicate that the primary moisture flux is from the southwest, while data from the Kenton AWS for the two following years are more ambiguous.
Accumulation-rate time series constructed from the dated isotopic summer peaks show no long-term trend, except for a recent 5 year increase, which seems to be present in all the cores. Similarly, the total accumulation over the survey grid as a function of time is highly variable but shows no trend. These results suggest that the precipitation regime in central Greenland has not changed significantly over the past several decades.
Acknowledgements
We should like to acknowledge the efforts of the other members of our field party which included M. Jackson, K. van der Veen and M. Mabin. W. Boiler of the Polar Ice Coring Office was responsible for the excellent core recovery. We also appreciate the exceptional efforts of K. Swanson of the Polar Ice Coring Office, who coordinated the logistics. We are grateful to W. Dansgaard, N. Gundestrup and others at the Geophysical Isotope Laboratory in Copenhagen for making the oxygen-isotope measurements. We also thank E. Mosley-Thompson, D. Bromwich and F. Robasky for providing helpful suggestions to improve this manuscript. Finally, we should Uke to thank Dr J. Dionne, whose energy and persistence were critical ingredients in the successful completion of this project. This work was supported by the U.S. National Science Foundation Division of Polar Programs grants DPP-8520855 and DPP-8822081. This is Byrd Polar Research Center contribution C-805, and contribution number GISP 91-14 of the Greenland Ice Sheet Project (GISP2).
The accuracy of references in the text and in this list is the responsibility of the authors, to whom queries should be addressed.
Appendix Error Analysis
Equation (2) is valid for 100% core recovery up to the true surface. For this study where the core was recovered from the bottom of a 2 m pit, and some core loss occurred, Equation (2) can be re-written as:
where m k and r k are the mass and average radius of the kth core segment, respectively, L pit is the ice-equivalent length of the 2 m pit wall, and L(T)loss is the ice-equivalent length of the missing core above the reference horizon. To compute L pit density profiles along a wall of each drill pit were measured by sampling each prominent visible stratigraphic layer. If D is the diameter of the density tube, V is its volume and M is the number of density measurements, we can express the ice-equivalent pit depth as:
where mj is the mass of firn in the tube, dj is the depth at the top of the tube and again we have approximated the density of the unmeasured zones as the average of the densities above and below. To estimate the uncertainty in L pit, we assign an uncertainty of ± 1 g in measuring the mass on a triple-beam balance in the pit, and assume the volume of firn in the 500 cm3 density tube is known to within ± 10 cm3. This results in a fractional error of about 3%. Thus:
We find δL pit varies from about 3.5 to 6.5 cm. The uncertainty in L(T)loss can be treated in the same way, except that we must consider the two extreme cases for assigning core loss between the special drill runs identified by the core-loss algorithm (Reference Whillans and BolzanWhillans and Bolzan, 1988). In one extreme, all the core loss in the depth interval between two consecutive special runs is placed on the top of the depth interval; in the other extreme all the core loss is placed at the bottom of the interval. The uncertainty will differ for each extreme depending on the difference in densities between the core segments bracketing the core-loss zone, and we use the larger value of the two. We find that for all the sites the maximum uncertainty varies from about 1 to 3 cm. The uncertainty in the ice-equivalent length above the horizon T is then:
where the uncertainty δmk is estimated to be 1 g and the uncertainty δrk is one-half the range of the radius values measured for the fcth segment. Then from Equation (1):
where δT is the uncertainty in the age of the horizon. For time intervals of 20 years or so, this contribution from δT amounts to less than 1 % of the total uncertainty in the accumulation rate. However, note that, if the interval is defined by two isotopic peaks instead of the surface and a peak, δT is at least ±4 months, not including other ambiguities in the dating of the peaks. If accumulation rates are to be computed over short time intervals, it is more useful to reference the values to accumulation years, defined as the time interval between adjacent isotopic peaks or troughs.