Brief communication: Firn data compilation reveals the evolution of
the firn air content on the Greenland ice sheet

Abstract. The firn covering the Greenland ice sheet interior can retain part of the surface melt, buffering the ice sheet’s contribution to sea level, but its characteristics are still little known. Using remote-sensing observations from 2000–2017, we estimate that firn covers 1,405,500 ± 17,250 km2 of the ice sheet. We present 344 firn-core-derived observations of the top 10 m firn air content (FAC10), indicative of the firn’s meltwater retention capacity. FAC10 remained stable in the coldest 74 % of the firn area during 1953–2017, while FAC10 decreased in the warmest and driest 12 % of the firn area between 1997–2008 and 2011–2017, resulting in a loss of 180 ± 78 km3 (−26 ± 11 %) of air from the near-surface firn.



Introduction 30
More than half of the Greenland ice sheet's current contribution to sea level rise stems from surface melt and subsequent runoff (van den Broeke et al., 2016). During summer, surface melt occurs across a large area of the ice sheet, even reaching the highest ice sheet elevations during extremely warm summers like 2012 (Nghiem et al., 2012). Most of the surface meltwater produced in the firn-covered interior of the ice sheet is refrozen into the snow and firn and does not immediately contribute to sea-level rise . This retention capacity of the firn is controlled by (i) the areal extent of the 35 firn; (ii) the firn air content (FAC); iii) firn temperature; and iv) firn permeability. In this study, we use in-situ and remotelysensed observations to estimate the firn's extent and its air content in the upper 10 m.

2
The first attempt to delineate the ice-sheet firn area into characteristic zones dates back to the 1950s (Benson, 1962). Later studies delineated the firn area using satellite remote sensing (e.g. Nolin and Payne, 2007) but had limited spatial and temporal coverage. Recently Fausto et al. (2018) published maps of remotely-sensed end-of-summer snowline but did not discuss the simple implication of these snowlines for the firn area extent.
The FAC is the volume of air contained in the firn from the surface to a certain depth per unit area. It indicates for that depth 5 range the maximum volume available to store percolating meltwater either in liquid or refrozen form . Harper et al., (2012Harper et al., ( ) used two years (2007Harper et al., ( and 2008 of observations from 15 sites along the western slope of the ice sheet to estimate the spatially integrated FAC of the entire percolation area in spite of the diversity of firn structures across the ice sheet (Forster et al. 2014;Machguth et al., 2016). Estimation of the FAC was also made using firn models forced by regional-climate-model (RCM) output (e.g. Ligtenberg et al., 2018) but model results differ significantly from each other in 10 terms of firn density and therefore of FAC (Steger et al., 2016).
The depth to which meltwater may percolate, and therefore to which FAC must be calculated, varies with melt intensity and firn permeability. Heilig et al. (2018) did not observe percolation deeper than 2.3 m below the surface at 2300 m a.s.l. in west Greenland indicating that, at their study site, only the near-surface FAC was being used to store meltwater. In a warmer region, ~400 km to the north and at 1555 m a.s.l., Humphrey et al. (2012) observed percolation below 10 m meaning that the 15 FAC of the whole firn column, from the surface to pore-close-off depth, might be used for meltwater retention. Nevertheless, Machguth et al. (2016) showed that percolation depth may not increase linearly with melt intensity and that refrozen meltwater can reduce percolation, subsequently preventing meltwater from accessing part of the FAC. Due to the absence of an ice-sheet-wide estimation of percolation depth and to the scarcity of firn observations covering the whole firn column, we here focus on the top 10 m of firn, for which numerous observations cores are available, and assume that the firn air content 20 in these top 10 m (FAC 10 ) represents the maximum volume of meltwater that can be stored in the firn.
In this study, we first estimate the firn area extent using remotely-sensed end-of-summer snow extent maps from Fausto et al., (2018). We then calculate the FAC 10 using a set of 344 firn cores collected between 1953 and 2017. We finally present the spatial distribution and where possible the temporal evolution of FAC 10 . Fausto et al. (2018) used surface radiance remotely sensed by the MODIS Terra satellite between 2000 and 2017 along with in-situ measurements of albedo at PROMICE automatic weather stations (Ahlstrøm et al., 2008) to estimate the end-ofsummer snow-covered and bare ice areas. Using these data, we determine the firn area, defined as the region where only snow has been detected during the entire 2000-2017 period. 30

Firn cores dataset
We gathered from the literature 324 firn-density profiles from cores that were at least 5 m long (Table 1 and S1 of the Supplementary Material) supplemented by 20 firn cores extracted in April-May 2016 and 2017 as part of the FirnCover campaigns and for which the density was measured at 10 cm resolution. Potential gaps in the density profile were filled using a logarithmic function of depth fitted to the available densities. When near-surface snow densities were missing, we 5 assigned a density of 315 kg/m 3  to the top centimetre of snow before the gap-filling.

Calculation of the FAC 10
For a discrete density profile, the FAC 10 in m 3 m -2 is calculated to depth as: where, for each section k, is the density and is the mass. The density of the ice formed after meltwater infiltration and 10 refreezing, , is set to 873 kg m -3 after Machguth et al. (2016). The value for was seen to vary within ±25 kg m -3 by Machguth et al., 2016 andHarper et al., 2012 used a value of 843 kg m -3 . Changing by ±25 kg m -3 leads to a variation of ±1 to 10% for FAC 10 values ranging from 5 to 1 m 3 m -2 . Addressing the variability of and its potential drivers is beyond the scope of this study.
With 121 cores shorter than 10 m in our dataset, we need to extrapolate shallow measurements to a depth of 10 m. For that 15 we find the 10+ m core that best matches the FAC vs. depth profile of the shallow core, with lowest Root Mean Squared Difference (RMSD) amongst all available cores, and append the bottom section of this 'twin' core to the FAC profile of the shallow core (see Figure S1 of the Supplementary Material). When testing this methodology on the available 10+ m long cores from which we remove the deepest 3 m of the FAC profile, find a mean difference between extrapolated and real FAC 10 inferior to 1% and a RMSD of 3-15% for FAC 10 values ranging from 5 to 1 m 3 m -2 . 20 The accuracy of the firn and infiltration ice densities as well as the effect of spatial heterogeneity can be assessed by comparing FAC 10 measurements located within 1 km ( Figure S2 of the Supplementary Material). The standard deviation of co-located FAC 10 observations is below 0.15 m 3 m -2 for the majority of sites (20 of 27). We therefore attach to any FAC 10 measurement an uncertainty of ±0.3 m 3 m -2 , twice the standard deviation, ±6 to 30% for FAC 10 values ranging from 5 to 1 m 3 m -2 . 25

Zonation of the firn air content
FAC 10 depends on the site's temperature and accumulation history and on meltwater refreezing at depth (Reeh, 2008). We extract each core site's long-term   temperature ̅̅̅ respectively from Box (2013) and Box et al. (2013) (Figure 1a). We also find ̅ and ̅̅̅ for all locations that occur within the firn area derived in this study (outlined area in Figure 1b).
Borrowing the terminology from Benson (1962), we define three regions where FAC 10 shows markedly different behaviours: (1) the Dry Snow Area (DSA, amber area in Figure 1a); (3) the Low Accumulation Wet Snow Area (LAWSA, red area in Figure 1a); (3) the High Accumulation Wet Snow Area (HAWSA, green area in Figure 1a). The DSA encompasses low 5 temperature regions of high altitude or latitude where melt is rare and where FAC 10 can be explained by well-known dry-firn compaction equations dependent on ̅ and ̅̅̅ (amber area in Figure 1a). Toward higher ̅̅̅ (lower altitude and/or latitude) two patterns are evident. At lower ̅ sites in the LAWSA, more scatter appears in FAC 10 , and a slope change occurs in the FAC 10 's temperature dependency (Figure 1c). At higher ̅ (in the HAWSA), FAC 10 remains higher than at locations with similar ̅̅̅ in the LAWSA and the slope of the temperature dependency is different from the one found in the DSA or 10 LAWSA (green area in Figure 1d).
The boundary between the colder (DSA) and warmer regions (LAWSA and HAWSA) can be defined as the ̅̅̅ where the dry firn densification cannot explain FAC 10 variations and increasing scatter appears in the FAC 10 values (Figure 1c). We set this threshold to ̅̅̅ = -16 o C as it is the temperature for which the standard deviation of FAC 10 within 1 o C-wide bins first exceeds 0.3 m 3 m -2 ( Figure 1c). Few firn cores are available in the transition zone from LAWSA to HAWSA and a boundary 15 could be anywhere between 543 mm w.eq./yr (core with highest accumulation in the LAWSA, Figure 1a) and 762 mm w.eq. yr -1 (core with lowest accumulation in HAWSA, Figure 1a). We chose the rounded value of 600 mm w.eq. yr -1 to separate LAWSA from HAWSA. The geographical delineations of the DSA, LAWSA and HAWSA are presented in Figure 1a.

Firn air content mapping
To map FAC 10 over the entire firn area from our collection of observations, we fit empirical functions of ̅ and ̅̅̅ to FAC 10 20 measurements and use these functions to predict FAC 10 anywhere in the firn area. We prefer this empirical approach to purely statistical approaches (e.g. kriging) or to the use of RCM and firn models which still do not accurately reproduce observations of firn densities and thus FAC 10 in the lower accumulation area (Steger et al., 2016;Langen et al., 2017). The form of these empirical functions is arbitrary but is tightly constrained by the numerous FAC 10 observations that they should fit. In this section, we briefly explain how we build these functions while further details and illustrations are available in 25 Figure S3 of the Supplementary Material.

Dry Snow Area
In the DSA, we use the steady-state firn densification model by Arthern et al. (2010) and we tune the surface snow density to 302 kg m -3 through least squares method to match the 209 available FAC 10 observations ( Figure S3

5
Material). We also investigated other densification laws and the one from Arthern et al. (2010) gave the best match with FAC 10 observations. We investigate the robustness of our method using a sensitivity analysis for each period separately. For 1000 repetitions, we randomly exclude four observations (respectively 9% and 11% of the observations in 1997-2008 and 2011-2017) and fit our 15 empirical function FAC 10 (̅ , ̅ ) to the remaining measurements. We then calculate the standard deviation of all possible FAC 10 at each (̅ , ̅̅̅ ) location and double it to quantify the 95% envelope of uncertainty that applies to any predicted FAC 10 in the LAWSA depending on (̅ , ̅̅̅ ). Since we do not consider that prediction can be more precise than field measurements, we set 0.3 m 3 m -2 as the minimum possible uncertainty.

High Accumulation Wet Snow Area 20
The HAWSA is described by only 15 firn cores drilled between 2010 and 2017 at 7 sites in the interior of the HAWSA, more than 20 km from the firn line. To overcome the scarcity of observations we use our remotely-sensed firn line as additional measurements where FAC 10 = 0 m 3 m -2 . We then fit our empirical function FAC 10 (̅ , ̅ ) as follows. In the surroundings of the core sites, we use a linear function of ̅̅̅ fitted through least squares to the FAC 10 observations: FAC 10 ( ̅̅̅ ) = −0.07 * ̅̅̅ + 3.4 (RMSD = 0.25 m 3 m -2 , Figure S3 in the Supplementary Material). The residuals of this fit did not show 25 any correlation with ̅ meaning that not enough measurements are available to disentangle the control of ̅ and ̅̅̅ on FAC 10 in the HAWSA. We then need to describe how the FAC 10 decreases from the core sites down to 0 m 3 m -2 at the remotelysensed firn line. We can make three estimates: i) a mid-range estimate where FAC 10 is bilinearly interpolated between the available firn cores and the firn line; ii) an upper-range estimate where FAC 10 follows the temperature dependency presented The Cryosphere Discuss., https://doi.org /10.5194/tc-2018-172 Manuscript under review for journal The Cryosphere Discussion started: 28 August 2018 c Author(s) 2018. CC BY 4.0 License. 6 above until the firn line where it drops abruptly to 0 m 3 m -2 ; and iii) a lower-range estimate where the FAC 10 drops to 0 m 3 m -2 shortly after the observations. These three surfaces are presented in more detail in Figure S3 in the Supplementary Material.
The mid-range estimate is considered as our most realistic estimation and is the empirical function FAC 10 (̅ , ̅ ) we use for the mapping of FAC 10 in the HAWSA while half the spread between the upper-and lower-range estimates quantifies the uncertainty applying to our FAC 10 map. 5

Delineation of the firn area
The firn area, illustrated in Figure 1b, covers at least an area of 1,405,500 km 2 or 78.5 % of the ice sheet (when compared to the contiguous ice extent from Citterio and Ahlstrøm (2013)). Moving the firn line 1 km in-or outwards (the resolution of the MODIS surface radiance product) suggest an uncertainty of ±17,250 km 2 (~1%). However, one should keep in mind that 10 the firn line is an idealized view of a patchy and gradual transition from bare ice to snow and firn that survived the summer melt in previous years (Benson, 1962). Spatial heterogeneity in melt and snowfall leave isolated snow and firn patches at the end of the melt season which might be missed by the method of Fausto et al. (2018). Using the average from 18 years of data gives a robust firn line and reduces the noise from local heterogeneities but also provides the absolute minimum extent of the firn area during 2000-2017. Ephemeral or thinner firn patches may exist outside of our strict delineation but we do not 15 believe they play an important role in the overall retention capacity of the firn.

Stable FAC 10 in the Dry Snow Area
The spatial distribution of FAC 10 in the DSA is shown in Figure 2. In that region representing 74% of the firn area, we find a spatiotemporally average FAC 10 value of 4.9 m 3 m -2 . We find a RMSD between predicted and observed FAC 10 (Figures 3b and 3c), confirms that the decrease in FAC 10 between the 1997-2008 and 2011-2017 is greater than the uncertainty associated to the two FAC 10 maps. The FAC 10 loss in the LAWSA is concentrated in a 60 km wide band above the firn line (Figure 2). Summing the FAC 10  Recent studies attributed the increasing near-surface firn densities and subsequent loss of FAC to increasing surface melt and 10 meltwater refreezing (de la Peña et al., 2015;Machguth et al., 2016;Graeter et al., 2018;Charalampidis et al., 2015).
However, firn density and FAC are also greatly dependant on annual snowfall (Herron and Langway, 1980) and a decrease in snowfall could also trigger a decrease in FAC 10 . Nevertheless, the lack of distributed observation of snow accumulation for the 1997-2017 period and the contradicting trends in precipitation given by the RCMs (Lucas-Picher et al., 2012;van den Broeke et al., 2016;Fettweis et al., 2017) make it impossible to quantify the contribution of snowfall into change in FAC 10 . 15 In our regional approach, we grouped measurements from different years, which made it impossible to define the precise years that were responsible for this loss of FAC 10 . However, repeated observations in western Greenland (Machguth et al., 2016;Charalampidis et al., 2015) indicate that 2010 and 2012 extreme melt seasons are responsible for the greatest changes in near-surface firn density and FAC 10 during the 1997-2017 period.
Finally, the loss of FAC 10 near the firn line likely turned firn areas into bare ice and triggered the upward migration of the 20 firn line. Unfortunately the relative uncertainty of our FAC 10 maps is greatest near the firn line (Figure 2 and 3) making it difficult to infer firn line migration from FAC 10 changes. Using our remotely-sensed firn line for all periods was a simplification but discussing its temporal evolution requires additional data and is beyond the scope of this study.

Potential deep meltwater percolation in the High Accumulation Wet Snow Area
In the HAWSA, comprising 14 % of the firn area, we calculate an average FAC 10 value of 2.4 m 3 m -2 . The spatial 25 distribution of FAC 10 is shown in Figure 2 while our estimated uncertainty map is presented in Figure 3. Summing over the HAWSA, we calculate that 560 ±154 km 3 of air exists within the top 10 m of firn.
The markedly (up to 90%) higher FAC 10 in the HAWSA compared to the LAWSA for any given ̅ (Figure 1) indicates that for similar surface melt (of which ̅̅̅ is a proxy) the firn in the HAWSA was less dense and potentially less saturated by ice than in the LAWSA. In line with observations at aquifer sites located in the HAWSA (Forster et al., 2014), we hypothesize 30 that in this region meltwater percolation may exceed 10 m, facilitated by high snow accumulation insulating the firn from the The Cryosphere Discuss., https://doi.org /10.5194/tc-2018-172 Manuscript under review for journal The Cryosphere Discussion started: 28 August 2018 c Author(s) 2018. CC BY 4.0 License. 8 winter cold. This deep percolation also implies that FAC 10 may not be the most adequate indicator of the firn's retention capacity in the HAWSA as additional retention may occur at deeper than 10 m. The fate of this deep meltwater whether it refreezes below 10 m, flows to the nearest aquifer, or reaches the bed, remains unknown and requires further investigation.

Total firn meltwater storage capacity
The total air content in the DSA during 1953-2017, in the LAWSA during 2011-2017 and in the HAWSA during 2010-2017 5 can potentially accommodate 5,480 ± 420 Gt of meltwater (15 ± 1.2 mm sea level equivalent). Harper et al. (2012) estimated that the firn located in the long-term percolation area (as modelled by a RCM) could potentially store between 322 ±44 Gt (when considering the top 10 m of firn) and 1 289 −252 +388 Gt (if considering the FAC to pore close-off depth). The smaller retention capacity found by Harper et al. (2012) mainly owes to the smaller area considered in their study. With our distributed approach, it is now possible to determine each year which areas of the firn are available to store meltwater and 10 the available FAC 10 given the extent of surface melt.

Uncertainty of long-term average climatic conditions
To investigate how uncertainties in ̅ and ̅̅̅ impacts our FAC 10 maps, we repeat our procedure using the 1979-2014 ̅ and ̅̅̅ predicted by the Modèle Atmosphérique Régional (MAR; Fettweis et al., 2017), as illustrated in Figures S4 to S7 of the Supplementary Material. The MAR-derived FAC 10 fits equally well (within measurements uncertainty, RMSD < 0.03 m 3 m -15 2 ) the FAC 10 observations. Since Box et al. (2013) gives 2 m air temperature and MAR gives surface temperature the value of ̅̅̅ threshold between DSA and LAWSA/HAWSA is adjusted to -20 o C. Nonetheless differences between the two FAC 10 predictions exist in areas where the spatial pattern of temperature and accumulation in Box (2013) and MAR differ, especially in the southern and eastern firn area. In these regions fewer cores are available to constrain the FAC 10 estimates.
More observations in these sparsely observed regions would therefore not only improve FAC10 estimates, but also elucidate 20 which ̅ and ̅̅̅ source best describes the spatial pattern in FAC 10 .

Isolated FAC 10 measurements
In the LAWSA and HAWSA, our dataset also contains 40 firn cores that were too isolated in space and time to be used FAC 10 in mapping. Isolated measurements nevertheless provide insight on how the FAC 10 might have been at a given place and time. Renaud (1959) reported a measurement in the LAWSA with a FAC 10 ~30% higher (+ 1 m 3 m -2 ) than the one 25 measured in 2006-2007. Conclusions from a single measurement are dubious, but it still indicates that the FAC 10 may have The Cryosphere Discuss., https://doi.org /10.5194/tc-2018-172 Manuscript under review for journal The Cryosphere Discussion started: 28 August 2018 c Author(s) 2018. CC BY 4.0 License. 9 been significantly higher in the 1950's. Ten observations in the LAWSA in the 1980's by van der Veen, et al. (2001), did not appear systematically different than our calculated FAC 10 for the 1997-2008 period, suggesting that FAC 10 was similar in the LAWSA during these two periods. A last measurement raises questions in the HAWSA: the core 6348 by Mosley-  drilled in 1998 indicates a FAC 10 of 3.9 m 3 m -2 . It is a factor of 2.2 higher than our mid-range FAC 10 estimate based on measurements from 2010-2017 and indicates either a drastic decrease in FAC 10 between 1998 and 2010-5 2017 or inaccuracies of our methodology for that location. A re-survey of that location would be of great interest even though the 1998 core, being a single point in space, can always be suspected to be anomalous or not representative.

Conclusions
Our study provides, for the first time, a delineation of the firn area of the Greenland ice sheet over the 2000 -2017 period.
Using remote-sensing observations from 2000 to 2017, we estimate that firn covers 1,405,500 ± 17,250 km 2 or 78.5% of the 10 ice sheet. This result allows further study of possible migration of this boundary in the past or in the future. Additionally, we present a collection of 344 firn cores spanning 65 years from which the firn air content from the surface to 10 m depth (FAC 10 ) could be calculated. We identify three regions on the firn area in which FAC 10 where we could fit empirical functions of long-term accumulation and temperature averages (̅ and ̅̅̅ ) to FAC 10 measurements and explain the spatiotemporal evolution of FAC 10  Area (where ̅̅̅ >= -16 o C and ̅ < 600 mm w.eq. yr -1 ) we find an average FAC 10 of 2.9 m 3 m -2 during the 2000-2017 period. FAC 10 observations also indicated that meltwater may percolate deeper than 10 m from the surface making FAC 10 20 insufficient to describe the retention capacity of the firn there. In a similar way, Machguth et al. (2016) showed that under conditions not completely understood, ice formation may prevent meltwater from accessing the entire top 10 m of firn.
Therefore, further investigation of the firn permeability will help to understand how much of the FAC is used for meltwater retention. The firn area delineation and FAC 10 dataset and maps can be used to constrain firn models and monitor the future evolution of the firn and are available for download at www.promice.dk. 25

Acknowledgement and data availability
This work is part of the Retain project funded by the Danish Council for Independent research (Grant no. 4002-00234) and the Programme for Monitoring of the Greenland Ice Sheet (www.PROMICE.dk). The FirnCover field campaigns were funded by the NASA grant NNX15AC62G. Achim Heilig was supported by DFG grant HE 7501/1-1. The source code is available at github.com/BaptisteVandecrux/FAC10_study. We thank Hubertus Fischer (Climate and Environmental Physics, 30 The Cryosphere Discuss., https://doi.org /10.5194/tc-2018-172 Manuscript under review for journal The Cryosphere Discussion started: 28 August 2018 c Author(s) 2018. CC BY 4.0 License.