Our surface density dataset consists of 200 point observations, along with the geographic location, annual air temperature and annual accumulation rate for these locations. The oldest surface density data were collected by Benson (1962) in 1954 in Northwest Greenland at latitudes between 70° to 77° N (Appendix C). These measurements include annual accumulation rates and 10 m firn temperatures reported by Mock and Weeks (1965) (Appendix A). Data from both the percolation and ablation areas of the southern and western ice sheet sections near 61.3° N (Nordbo Gletscher) and 69.7° N (Paakitsoq), respectively, were collected by Braithwaite et al. (1982, 1994). Our dataset also includes data from the Program for Arctic Regional Climate Assessment (PARCA) (Mosley-Thompson et al., 2001; Thomas and Investigators, 2001). Further, the SUrface Mass balance and snow depth on sea ice working grouP (SUMuP) provided accumulation rates, snow depths and density values at various sites on the ice sheet (Koenig et al., 2013; Montgomery et al., 2018), including observations from a study of Greenland accumulation (Hawley et al., 2014) and firn aquifers (Forster et al., 2013; Koenig et al., 2014; Miège et al., 2016). We gathered annual air temperatures, accumulation rates, and density observations from snow pits and firn cores from the Greenland Climate Network (GC-Net) (Steffen et al., 1996), the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) (Van As et al., 2016b), and the Arctic Circle Traverses (ACTs) (e.g., Machguth et al., 2016). Lastly, we also included observations by Schaller et al. (2016) from the NEEM to EGRIP traverse, and from the López-Moreno et al. (2016) Greenland circumnavigation. Accumulation rates in the database are not long-term averages, but represent the preceding year's snowfall. Figure 1 provides a map of all measurement locations. All data are available as Supplementary Material. Figure 2a illustrates that 28% of the observations were taken in the mid-1950s, only 2% were taken in the 1980s and 1990s, while 70% were obtained between 1999 and 2016. 94% of the measurements were gathered at elevations exceeding 1,000 m above sea level (Figure 2d).

Defining the surface layer as the upper 0.1 m of snow yields that in most cases the surface layer was deposited in multiple snowfall events, except for areas located at relatively low elevations in the south and southeast of the ice sheet, where individual precipitation events typically produce more than 0.1 m of snow (Burgess et al., 2010). Where possible, the annual air temperature was calculated as the average over the 365 days prior to the date for which the surface snow density was determined. Where air temperature measurements are not available, i.e., for the older data by Benson (1962) and Braithwaite et al. (1994), but firn-core temperatures were, we use 10 m firn temperature as annual air temperature following e.g., Reeh et al. (2005), Polashenski et al. (2014), and Kuipers Munneke et al. (2015). This is a fair approximation since 10 m firn temperatures reflect the conductive temperature wave propagation in places with little or no melt (Benson, 1962). Though valid for the earlier observations in our dataset, recent increases in ice sheet melt area have reduced the dry snow facies of the ice sheet (McGrath et al., 2013) and therefore the applicability of this methodology.

Commonly, snow/firn was sampled in snow pits using a fixed volume cutter at 0.05–0.1 m vertical resolution. These samples were weighed using a variety of scales. When density data were derived from a core, the snow was extracted from the core barrel and typically sub-sampled into 0.1 m sections before being weighed. Conger and McClung (2009) investigated measurement errors of several different density cutters and conclude that measurement accuracy was within 3–12%. They also conclude that the absolute measurement uncertainty is within 11% of true density. A discussion of density cutters by Proksch et al. (2016) reaches a similar uncertainty of 9%. For the data in our database, sampling uncertainty is not documented in any of the field campaigns, however it seems reasonable to assume that surface snow density is known within 10%. Typically, this measurement uncertainty is smaller than the spatial variability in surface snow density in the vicinity of the measurement location (e.g., Proksch et al., 2015). We argue that the point measurements in our dataset do not represent fresh snow, as the persistent katabatic winds in Greenland compact surface snow within days after snowfall (e.g., Liston et al., 2007).

We test a surface snow density parameterization for the Greenland ice sheet that is dependent on temperature, similar to commonly used parameterizations by Kuipers Munneke et al. (2015) and Reeh et al. (2005). We assume a linear dependence of surface snow density (ρ) in kg m⁻³ on annual air temperature (T_a) in °C, in what we refer to as parameterization P1:

We determine the fit coefficients by orthogonal linear regression to all available T_a values in our dataset, and find a best fit for A = 362.1 and B = 2.78 (Table 1) for the top 0.1 m of snow. Kuipers Munneke et al. (2015) determined the coefficients of Equation (1) using annual surface temperature T_s in °C simulated by RACMO2.3 and the average density of the uppermost 1 m of snow/firn, and found what we here refer to as parameterization P2:

Reeh et al. (2005) derived surface snow density as a function of the 10-m firn temperature (T_f) from the near-surface part of their depth-density profiles by determining the load at 5 m depth, as calculated by their model, so that it fits the corresponding load derived from the measured depth-density profiles (parameterization P3):

There is a ca. 40% overlap between our dataset and the data feeding into the Kuipers Munneke et al. (2015) and Reeh et al. (2005) parameterizations that stems from them also using the Benson (1962), Braithwaite et al. (1994), and PARCA (Mosley-Thompson et al., 2001) datasets.

To highlight the importance of sampling depth ranges in producing an observationally-based boundary condition for firn models in Greenland, we also test P1 (Equation 1) using the average density of the top 0.2 and 0.5 m of snow/firn in our analysis (Table 1). We theorize that, by using density data obtained as close to the surface as possible, we avoid introducing a systematic bias due to compaction. Yet by focusing only on the top layer of snow/firn, we likely introduce more scatter in our results due to additional variability by single weather events. We investigate such considerations below.

We use a smaller depth range to better represent surface snow density than previous studies. Assessing density closer to the surface is important for producing a more accurate upper boundary condition to be used in firn evolution models that would produce too high firn densities along the entire depth profile. Figure 4, Table 3 confirm that using relatively large depth ranges in determining a surface snow density parameterization results in overestimated values by Kuipers Munneke et al. (2015) and Reeh et al. (2005). Our smallest tested depth range of 0–0.1 m reveals larger natural variability, but would not introduce a considerable systematic bias in firn evolution modeling even if a vertical grid resolution finer than 0.1 m is used. In surface mass balance modeling, the choice of vertical resolution of the subsurface directly influences the calculation of key variables, such as the meltwater retention capacity of the snow/firn column. The more variable density in the top 0.1 m of snow compared to the top 0.2 m (factor of 1.4 more variable) or 0.5 m (1.8), is due to the influence of single precipitation events and subsequent weather forcing. We contend that this increased variability is preferable over the introduction of a systematic bias in surface mass balance modeling.

The top of the snowpack compacts rapidly after snowfall (e.g., Brun et al., 1997; Liston et al., 2007), as the crystal structure of freshly deposited snow breaks down within days due to wind and redistribution of drifting snow (e.g., Kotlyakov, 1961; Kojima, 1967; Pahaut, 1976). Surface snow densification by wind, which generally only influences the top 0.1 m, becomes insignificant after a few days (Brun et al., 1997). For most or all observations in our dataset, we can safely assume that wind compaction has occurred already. Therefore, our dataset and resulting products should not be used in models to prescribe or validate fresh snow densities (e.g., Vionnet et al., 2012), but rather to define the upper boundary condition (i.e., minimum density) in firn evolution models that do not calculate micro-scale snow physics and densification by wind, snow drift and redistribution.

In regions where large snowfall events occur, such as in south Greenland, density measurements of the top 0.1 m of snow may reflect the conditions during one snowfall event and subsequent weather-dependent densification prior to measurement. All of the snow-density measurements in our database were taken in spring and summer, meaning that our average and parameterization may be seasonally biased. Dibb and Fahnestock (2004) investigated the seasonality of the surface snow density at Summit in Greenland, and found a seasonal standard deviation of 30% in density in the top 0.03 m of snow as determined from 22 measurements during a two-year period. However, seasonal variation in the surface snow density is likely to increase with elevation (Brun et al., 1997) with standard deviation values lower than 30% in regions away from the three dome sites in Greenland where persistent katabatic winds and their influence on snow compaction do not occur (e.g., Noël et al., 2014). In general, katabatic winds are strongest in winter due to surface radiative cooling, and at lower elevations (below 2,000 m above sea level) due to larger surface slopes (e.g., Van As et al., 2013; Noël et al., 2014), resulting in wind-packing of fresh surface snow within days.

Higher air temperatures result in higher snow and firn densities through increased compaction (Zwally and Li, 2002). It is therefore desirable to ensure that parameterizations of surface snow density remain appropriate even as the climate changes (e.g., Reeh et al., 2005; Morris and Wingham, 2014). Studies of Greenland accumulation rates and firn properties document a recent densification in the overall firn column and attribute it to climate warming (e.g., de la Peña et al., 2015; Charalampidis et al., 2016a; Machguth et al., 2016; Overly et al., 2016). If we assume that temperature-dependent densification processes are responsible for the transformation of freshly-fallen snow to the surface snow densities of our dataset, the inclusion of temperature as a variable in a parameterization (Equation 1) explicitly accounts for atmospheric warming. In this case, a parameterization is better capable of representing changing surface conditions due to climate variability. For instance, our temperature-dependent parameterization suggests that the observed 2.7°C warming at Summit over the period 1982–2011 (McGrath et al., 2013) had lead to a local surface snow density increase of 8 kg m⁻³. But even a large temperature increase of 10°C anywhere in Greenland would only cause a densification of 28 kg m⁻³ in the top 0.1 m of snow, which is smaller than the 32 kg m⁻³ measurement uncertainty and 44 kg m⁻³ standard deviation of the dataset (Table 2). For larger depth ranges the temperature sensitivity (B-values in P1, see Table 2) is considerably smaller and thus more insignificant given the measurement uncertainty and natural variability. The insignificant densification as a result of warming supports the notion that temperature is a poor predictor of the variability of surface snow density in the top 0.1, 0.2, and 0.5 m of snow/firn, and that using a constant value may be preferable in some applications.

The choice of a surface snow density boundary condition influences calculations of available pore space by models simulating the surface mass budget of the Greenland ice sheet. Steger et al. (2017) discussed the limitations and inaccuracies of their Greenland ice sheet surface mass budget simulations by regional climate model RACMO2.3, and conclude that the Kuipers Munneke et al. (2015) parameterization systematically overestimates surface snow density, impacting pore space available for refreezing at depth. Langen et al. (2015) and Charalampidis et al. (2016b) applied a constant surface snow density value of 330 kg m⁻³ in regional climate model HIRHAM5, while Langen et al. (2017) applied a parameterization depending on latitude, longitude and elevation, derived from our dataset, in a new model version. The latter study found that the parameterization yields an ice-sheet-wide average of surface snow density that is 7% lower than using a constant density value of 330 kg m⁻³, signifying a higher meltwater retention capacity in the snow and firn. Langen et al. (2017) also documented that the firn density profiles simulated by HIRHAM5 using their parameterization satisfactorily resemble measured profiles. Yet based on our own findings we suspect that using a constant surface snow density value of 315 ± 32 kg m⁻³ as boundary condition, a value 5% lower than that used by Langen et al. (2015), should perform equally well in Greenland-wide surface mass budget simulations.

Using our dataset for the top 0.1 m of snow, as opposed to using those for larger depth ranges, comes at the cost of a higher variability (standard deviation in Table 2) due to a larger influence of meteorology-dependent processes like snow drift (Brun et al., 1997; Hörhold et al., 2011; Koenig et al., 2016). The larger variability could also stem from snowfall events depositing more than 0.1 m of snow. Layers below the top 0.1 m are not much influenced by wind compaction, snow drift and redistribution, and will primarily be subject to the less efficient densification through rounding or settling of snow grains from vapor fluxes in the subsurface layers (e.g., Albert and Shultz, 2002). Surface mass budget models using our constant value of 315 kg m⁻³ for the top 0.1 m of snow may therefore misrepresent relatively low-density layers below 0.1 m depth deposited during large snowfall events. Regions where snowfall may exceed 0.1 m in single events are typically located at lower elevations on the southern and southeastern parts of the ice sheet (e.g., Burgess et al., 2010).

Our dataset has sparse coverage in the northern and eastern sectors of the ice sheet, possibly introducing a spatial bias in our results. Figure 5 illustrates that the elevation distribution of our measurement locations broadly reflects the overall area-elevation distribution of the ice sheet as determined from the GIMP digital elevation model (Howat et al., 2014). Some lower elevation ranges (1,000–1,750 m above sea level) are relatively underrepresented in our dataset, while some higher elevation ranges are comparatively overrepresented. Our parameterization could benefit from acquiring additional measurements from elevations between 1,000 and 1,750 m above sea level, i.e., in the lower percolation area of the ice sheet (Benson, 1962). The lower percolation area is considered crucial for properly determining the surface mass budget, as firn properties influencing meltwater retention capacity vary substantially across the ice sheet (Van As et al., 2016a; Langen et al., 2017).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer, CM, declared a past collaboration with two of the authors, JB and LK, to the handling editor.