Firn data compilation reveals widespread decrease of firn air content in western Greenland

Abstract. A porous layer of multi-year snow known as firn covers the
Greenland-ice-sheet interior. The firn layer buffers the ice-sheet
contribution to sea-level rise by retaining a fraction of summer melt as
liquid water and refrozen ice. In this study we quantify the Greenland
ice-sheet firn air content (FAC), an indicator of meltwater retention
capacity, based on 360 point observations. We quantify FAC in both the
uppermost 10 m and the entire firn column before interpolating FAC over the
entire ice-sheet firn area as an empirical function of long-term mean air
temperature (Ta‾) and net snow accumulation
(c˙‾). We estimate a total ice-sheet-wide FAC of 26 800±1840 km3, of which 6500±450 km3 resides within the
uppermost 10 m of firn, for the 2010–2017 period. In the dry snow area
(Ta‾≤-19 ∘C), FAC has not changed
significantly since 1953. In the low-accumulation percolation area
(Ta‾>-19 ∘C and c˙‾≤600 mm w.e. yr−1), FAC has decreased by 23±16 % between
1998–2008 and 2010–2017. This
reflects a loss of firn retention capacity of between 150±100 Gt and
540±440 Gt, respectively, from the top 10 m and entire firn column.
The top 10 m FACs simulated by three regional climate models (HIRHAM5,
RACMO2.3p2, and MARv3.9) agree within 12 % with observations. However,
model biases in the total FAC and marked regional differences highlight the
need for caution when using models to quantify the current and future FAC and
firn retention capacity.


temporary sea-level rise (Bindoff et al., 2013;Nerem et al., 2018). Over half of this GrIS mass loss stems from summer surface melt and subsequent meltwater runoff into the ocean (van den Broeke et al., 2016). While most meltwater runoff originates from the low-elevation ablation area, the surface melt area is now expanding into the high-elevation firn-covered interior of the GrIS (Mote et al., 2007;Nghiem et al., 2012). Rather than flowing horizontally, most of the meltwater produced at the surface of the firn area percolates vertically into the underlying firn where it refreezes and thereby does not contribute to sea-level rise . Hence, the meltwater retention capacity of Greenland's firn is a non-trivial parameter in the sea-level budget.
Assessing meltwater retention capacity of the firn in Greenland requires knowledge of both the extent of the firn area, as well as the spatial distribution of depth-integrated firn porosity or firn air content (FAC). The extent of the firn area can be tracked using the firn line, which Benson (1962) described as "the highest elevation to which the annual snow cover recedes during the melt season". Recently, Fausto et al. (2018a) updated the methods from Fausto et al. (2007) and presented maps of remotely sensed end-of-summer snow lines over the 2000-2017 period. These maps effectively provide an annual delineation of Greenland's firn area. FAC is the integrated volume of air contained within the firn from the surface to a certain depth per unit area Ligtenberg et al., 2018). FAC quantifies the maximum pore volume available per unit area to retain percolating meltwater, either in liquid or refrozen form . Previously, ice-sheet-wide firn retention capacity has been estimated using simplifying assumptions (Pfeffer et al., 1991) or unconstrained regional climate model (RCM) simulations . Harper et al. (2012) provided the first empirical estimate of the firn's meltwater retention capacity in the GrIS percolation area using 3 years of observations (2007)(2008)(2009) at 15 sites in western Greenland. While pioneering, their approach did not acknowledge the GrIS's diverse firn regimes Machguth et al., 2016). Ligtenberg et al. (2018) provided an RCM simulation of FAC that generally compares well against observations in 62 firn cores but substantially underestimated FAC in the western percolation area.
The depth to which meltwater may percolate, and therefore the depth range over which FAC must be integrated to constrain meltwater retention capacity, varies with melt intensity and firn permeability (Pfeffer et al., 1991). This makes the maximum depth of meltwater percolation both temporally and spatially variable, as highlighted by the following studies. Braithwaite et al. (1994) and Heilig et al. (2018) reported meltwater refreezing within the top 4 m of firn in western Greenland, respectively, at ∼ 1500 m a.s.l. during the summer of 1991 and at 2120 m a.s.l. during the 2016 melt season. Both studies indicate that, at specific sites and years, meltwater is stored in near-surface firn. However, firn temperature measurements in 2007-2009 at 1555 m a.s.l. in western Greenland  as well as the presence of firn aquifer at a depth greater than 10 m in southeastern Greenland (Miège et al., 2016) both show that meltwater can percolate below 10 m depth in the firn. This deep percolation implies that, for certain firn conditions and given sufficient meltwater, the FAC of the total firn column, from the surface to the firn-ice transition, may be used for meltwater retention. Finally, Machguth et al. (2016) show that percolation depth may not increase linearly with meltwater production, and instead low-permeability ice layers can limit even abundant meltwater from percolating into the entire firn column. Given the complexity of meltwater percolation and the paucity of percolation observations, reasonable upper and lower bounds of the meltwater retention capacity of firn can be estimated by determining FAC through the total firn column (FAC tot ) and within the uppermost 10 m of firn column (FAC 10 ), respectively . FAC tot is also valuable information to convert remotely sensed surface height changes into mass changes (Sørensen et al., 2011;Simonsen et al., 2013;Kuipers Munneke et al., 2015a).
In this study, we first compile a dataset of 360 firn density profiles, collected between 1953 and 2017, and quantify the observed FAC. We then extrapolate these point-scale observations across the entire GrIS firn area as empirical functions of long-term mean air temperature and mean snow accumulation. The point observations are thereby used to resolve the spatial distribution of FAC but also, where possible, its temporal evolution. We use a simple extrapolation to estimate FAC tot from FAC 10 where firn cores do not extend to the firn-ice transition. Spatial integration of FAC 10 and FAC tot over the firn area permits the estimation of the lower and upper bounds, respectively, of the GrIS firn meltwater retention capacity. Finally, we evaluate the FAC simulated by three RCMs that are commonly used to evaluate ice-sheetwide firn meltwater retention capacity but that have never been compared to such an extensive firn dataset.
2 Data and methods

Firn core dataset and firn area delineation
We compiled 340 previously published GrIS firn density profiles of at least 5 m in depth (Table 1). To these, we added an additional 20 cores extracted in 2016 and 2017, for which firn density was measured at a 10 cm resolution following the same procedure as Machguth et al. (2016). When nearsurface snow densities were missing, we assigned a density of 315 kg m −3 (Fausto et al., 2018b) to the top centimetre and interpolated over the remaining gaps in density profiles using a logarithmic function of depth fitted to the available densities.
We use the end-of-summer snow lines from Fausto et al. (2018a)  snow is always detected during the 2000-2017 period, is taken to represent the GrIS's current firn area. Moving this firn line 1 km inward or outward, the resolution of the product from Fausto et al. (2018a) suggests an uncertainty of ±17 250 km 2 (∼ 1 %). Additional uncertainty applies to the margin of the firn area where transient firn patches may exist outside of our delineation. Owing to the inherent thinness of firn at the lower elevation boundary of the firn area, we expect these omitted firn patches to play a negligible role in the overall meltwater retention capacity of the firn area.

Calculation of FAC 10
For a discrete density profile composed of N sections and reaching a depth z, the FAC in metres is calculated as fol-lows: where, for each depth interval k, ρ k is the firn density and m k is the firn mass. ρ ice is the density of the ice, assumed to be 917 kg m −3 . With 121 cores shorter than 10 m in our dataset, we extrapolate shallow measurements to a depth of 10 m. We do this by finding the core longer than 10 m that best matches the FAC-versus-depth profile of the core shorter than 10 m with the lowest root-mean-squared difference (RMSD) amongst all available cores. We then append the bottom section of this core longer than 10 m to the FAC profile of the core shorter than 10 m (see Fig. S1 of the Supplement). When testing this methodology on the available cores longer than 10 m, from which we remove the deepest 3 m of the FAC profile, we find a mean difference between extrapolated and real FAC 10 < 1 % and an RMSD of 0.15 m.
We assess the accuracy of the firn density measurements, as well as the effect of spatial heterogeneity, by comparing FAC 10 measurements located within 1 km and collected in the same year (Fig. S2). A standard deviation below 0.15 m is found in the majority of the co-located and contemporaneous FAC 10 observations (20 of 27 groups of comparable observations). We correspondingly assign an uncertainty of ±0.3 m, twice this standard deviation, to FAC 10 measurements.

Zonation of firn air content
The FAC 10 is calculated from firn density, which depends, among other parameters, on the local near-surface air temperature and snowfall rate (Shumskii, 1964). Air temperature is a proxy for summer melt and subsequent refreezing within the firn, as well as firn temperature and compaction rates. Through these processes, increasing air temperature acts to decrease FAC (Kuipers Munneke et al., 2015b). On the other hand, snow accumulation introduces low-density fresh snow at the surface. Increasing snowfall thus acts to increase FAC. To put our FAC 10 measurements in their climatic context, we extract the long-term (1979-2014) average annual net snow accumulationċ (snowfall − sublimation) and air temperature T a for each FAC 10 measurement location from the nearest 5 km 2 cell of the Modèle Atmosphérique Régional (MARv3.5.2; Fettweis et al., 2017).
Following the terminology of Benson (1962), we define three regions where FAC 10 shows distinct regimes: (1) the dry snow area (DSA, yellow area in Fig. 1a), (2) the lowaccumulation percolation area (LAPA, red area in Fig. 1a) and (3) the high-accumulation percolation area (HAPA, green area in Fig. 1a). The DSA encompasses low temperature regions of high altitude and/or latitude where melt is uncommon and where FAC 10 can be related by a linear function of T a (yellow markers in Fig. 1c). Two distinct firn regimes emerge towards higher T a , meaning lower altitude and/or latwww.the-cryosphere.net/13/845/2019/ The Cryosphere, 13, 845-859, 2019 itude. Firstly, towards lowerċ, in the LAPA, more scatter appears in FAC 10 and the slope of the FAC 10 temperature dependency changes. Secondly, towards higherċ, in the HAPA, the few available FAC 10 observations describe a similar temperature dependency as in the DSA, even though they are in relatively warm regions where melt occurs. FAC 10 observations in the HAPA are up to 5 times higher than at locations with similar T a in the LAPA (Fig. 1c). The boundary that delineates the cold (DSA) and warm regions (LAPA and HAPA) can be defined as the temperature where an inflection occurs in the linear dependency of FAC 10 on T a (Fig. 1c). We interpret the slope break in the temperature dependence of FAC 10 as the upper limit of frequent meltwater percolation and refreezing within the firn which Benson et al. (1962) defined as the dry snow line. While the transition between cold and warm areas is gradual in practice, for our analysis we set this boundary to T a = −19 • C. Our LAPA and HAPA here stretch from the dry snow line to the firn line and therefore also include the so-called wet snow facies defined by Benson et al. (1962). The snowfall boundary that delineates the low-and highaccumulation percolation areas is more difficult to characterise. There are insufficient firn observations available along the transition from LAPA to HAPA. The snowfall boundary could be anywhere between 543 mm w.e. yr −1 (the highestaccumulation LAPA core, Fig. 1b) and 647 mm w.e. yr −1 (the lowest-accumulation HAPA core, Fig. 1b). Acknowledging this uncertainty, we chose the round value ofċ = 600 mm w.e. yr −1 to separate LAPA and HAPA. The spatial delineations of the DSA, LAPA and HAPA are illustrated in Fig. 1a.

FAC 10 interpolation
To interpolate point-scale observations of FAC 10 over the entire GrIS firn area, we describe FAC 10 observations using empirical functions of long-term mean air temperature and net snowfall. The derivation of these empirical functions is described in the following sections and an overview of their general form as well as the data used to constrain them is presented in Table 2.

Dry snow area
In the DSA, the 259 FAC 10 observations obtained between 1953 and 2017 can be approximated by a linear function of their local T a (Fig. 1c). This dependency is the same for the 19 FAC 10 observations from the upper HAPA available between 1981 and 2014. We consequently include these observations so that the linear relationship remains valid in the upper HAPA (Sect. 2.4.2). These 278 FAC 10 observations are then binned into four equal T a ranges to avoid the overrepresentation of clustered data (Fig. 2a). Eventually, a linear function of T a is fitted to the bins' average FAC 10 using a least-squares method to estimate the FAC 10 in the DSA: FAC 10 T a = −0.08 × T a + 3.27. (2) We assign to any FAC 10 estimated in the DSA using Eq.
(2) an uncertainty equal to twice the regression's RMSD: 0.4 m. Although FAC 10 is also dependent onċ, the residuals from Eq.
(2) do not present any correlation with their respectiveċ values. It indicates that because of the intrinsic co-variability ofċ and T a , most of the variations in observed FAC 10 can be explained using eitherċ or T a . Insufficient data are available to separate the role ofċ and T a in FAC 10 variations in the DSA. We therefore choose to use only T a in Eq. (2).

Percolation areas
In the LAPA and in the HAPA, FAC 10 observations exhibit a more complex dependency onċ and T a (Fig. 1b and c). Additionally, observations are unevenly distributed in space and time. Thus, to reveal the temporal trends in FAC 10 , the observation dataset is divided into two time slices that each contain enough FAC 10 observations to describe the spatial pattern of FAC 10 and constrain our empirical functions.  (2) 6 selected from the firn line in the HAPA Over the 2010-2017 period, 25 FAC 10 observations were made in the LAPA, stretching from the upper boundary of the LAPA down to the vicinity of the firn line. During that same period, 10 firn cores were collected in the HAPA. Unfortunately, in addition to their small number, the cores are located relatively far into the interior of the ice sheet and do not describe how the FAC 10 decreases in parts of the HAPA closer to the firn line. We consequently complement these firn cores with six sites, selected on the remotely sensed firn line, where FAC 10 is assumed to be null (Fig. S3). FAC 10 in the LAPA and HAPA during 2010-2017 is then described by a smoothed bilinear function of T a andċ fitted through least squares method to the available observations (Fig. 3b). We do not allow that function to exceed the linear function of T a that describes FAC 10 measurements in the DSA and in the upper HAPA (Eq. 2) or to predict FAC 10 below 0 m.
Prior to 2010, insufficient data are available to document the FAC 10 in the HAPA. In the LAPA, however, 35 observations were made between 2006 and 2008 and three cores were collected in 1998. These measurements are used to describe the FAC 10 in LAPA during the 1998-2008 period by a smoothed bilinear function of T a andċ. To ensure that our empirical function has realistic values towards the transition with the HAPA, we also include one core collected in the HAPA in 1998. We also include the previously described six locations from the firn line (Fig. 3a). Although observation locations in 1998-2008 and 2010-2017 can be different, few samples available at the same sites (e.g. Crawford Point, DYE-2) in both time slices confirm that FAC 10 changes are more likely due to a temporal evolution rather than from the different spatial coverage of each period's constraining dataset.
The empirical functions used to estimate the FAC 10 in the LAPA and HAPA (Fig. 3) We investigate the robustness of our empirical functions in the HAPA and LAPA using, for each period separately, the following sensitivity analysis. For 1000 repetitions, we apply four types of perturbations to the FAC 10 observations and then refit our empirical functions. The effect of the availability of measurements in the LAPA is tested by randomly excluding four observations in that region (16 % and 11 % of observations in 2010-2017 and 1998-2008, respectively). The effect of uncertainty in the firn line location in the (T a ,ċ) space is tested by adding a normally distributed noise with a mean of zero and standard deviation of 3 • C to the T a of firnline-derived FAC 10 (illustrated in Fig. S3). The effect of the uncertain FAC 10 value at the firn line is assessed by assigning to firn-line-derived points a random FAC 10 value between 0 www.the-cryosphere.net/13/845/2019/ The Cryosphere, 13, 845-859, 2019 and 1 m. Finally, the effect of the smoothing applied to the bilinear interpolation of FAC 10 measurements is assessed by modifying the amount of smoothing applied. Following 1000 repetitions of the above-mentioned four perturbations to the FAC 10 observations, we then calculate the standard deviation of all empirically estimated FAC 10 values within the (T a ,ċ) parameter space. We then double this standard deviation to approximate the 95 % uncertainty envelope for empirically estimated FAC 10 in the LAPA and HAPA. We set 0.3 m, the uncertainty related to FAC measurements (Sect. 2.2), as the minimum possible uncertainty on any empirically estimated FAC 10 .

Estimation of FAC tot
FAC tot should be integrated from the ice-sheet surface down to the depth where firn reaches the density of ice (Ligtenberg et al., 2018). This depth varies in space and time across the GrIS but is poorly documented. Additionally, the RCM HIRHAM5 (evaluated in Sect. 3.3) does not reach ice density at the bottom of its column in certain locations. We therefore calculate FAC tot as the vertically integrate FAC from the surface to a standard 100 m depth. Only 29 of our 360 firn observations reach depths greater than 100 m. We therefore complement these core observations with 13 ground-penetrating radar observations of FAC tot from Harper et al. (2012).Using the least squares method with an intercept of zero, we fit the following linear regression between FAC 10 and FAC tot ( Fig. 4): This function infers that FAC tot is approximately 410 % of FAC 10 . While we acknowledge this relation is straightforward, we highlight that it is statistically robust. We assign 3.6 m, twice the RMSD of the linear regression, as the typical uncertainty for an estimated FAC tot value that can, in theory, vary between 0 and ∼ 25 m. As a result of deriving FAC tot as a function of FAC 10 (Eq. 3), any change in FAC 10 between two dates implies a The Cryosphere, 13, 845-859, 2019 www.the-cryosphere.net/13/845/2019/ proportional change in FAC tot over the same time period. The use of this covariation neglects the fact that near-surface changes in the firn slowly propagate to greater depth with thermal conduction and downward mass advection (Kuipers Munneke et al., 2015b). We therefore note that, for a decreasing FAC 10 (see Sect. 3.2.1), our estimated change in FAC tot corresponds to the maximum possible change associated with the whole firn column having sufficient time to adapt to the new surface conditions.

Spatially integrated FAC and retention capacity
We define, for any ice-sheet region, the spatially integrated FAC as the cumulated volume of air within that region either in the top 10 m of firn or for the total firn column (top 100 m). The uncertainty associated with the empirically estimated FAC 10 and FAC tot at a given location is not independent of other locations because the same functions of T a andċ are applied across the GrIS. Consequently, we consider that the uncertainty of the mean FAC in a specific region is the mean of FAC uncertainty values therein and that the uncertainty of spatially integrated FAC is the sum of the uncertainty values in the considered region. We use the estimated FAC to calculate the meltwater retention capacity of the firn. Harper et al. (2012) defined the firn retention capacity as the amount of water that needs to be added to the firn to bring its density to 843 kg m −3 , the density of firn saturated by refrozen meltwater measured in firn cores.

Comparison with regional climate models
We compare our FAC 10 observations and spatially integrated FAC estimates to the firn products available from three RCMs: HIRHAM5, RACMO2.3p2 and MARv3.9. HIRHAM5 output is available at 5.5 km spatial resolution and is presented in Langen et al. (2017). Two versions of HIRHAM5 are used: linear parameterisation of surface albedo (hereafter referred as HH_LIN) and MODIS-derived albedo (hereafter referred as HH_MOD). RACMO2.3p2, presented by Noël et al. (2018), provides FAC at a 5.5 km resolution. MARv3.9, presented in Fettweis et al. (2017), only simulates FAC 10 because of its shallow subsurface domain and has a spatial resolution of 15 km.

Spatio-temporal distribution of FAC
In the DSA, we consider the absence of a temporal trend in the deviation between measured FAC 10 and FAC 10 estimated using the linear function of T a (Fig. 2b) as evidence of unchanging FAC 10 in that area between 1953 and 2017. This inference of widespread stable FAC in the DSA is confirmed at point scale by firn cores in our dataset taken at the same sites but decades apart, showing the same FAC (e.g. Summit Station, Camp Century). This result is also corroborated by recent firn modelling at weather stations located in the DSA (Vandecrux et al., 2018).
Using the 5 km × 5 km T a andċ grids from Fettweis et al. (2017) and the empirical functions presented in Fig. 3, we map the FAC 10 and its uncertainty across the GrIS firn area (Fig. 5). From these maps we calculate an average FAC 10 of 5.2 ± 0.3 m in the DSA over the 1953-2017 period and of 3.0 ± 0.4 m in the HAPA during the 2010-2017 period. Within the LAPA, we calculate an average FAC 10 of 3.9 ± 0.3 m during the 1998-2008 period, which decreases by 23 % to 3.0 ± 0.3 m in the 2010-2017 period. Spatially, the FAC 10 loss in the LAPA is concentrated in a 60 km wide band above the firn line in western Greenland (Fig. 5b).
We find that during the 2010-2017 period, the entire firn area contained 6500±450 km 3 of air within the top 10 m and up to 26 800 ± 1840 km 3 within the whole firn column (Table 3). About 83±5 % of this air content is found in the DSA, which represents 74 % of the firn area. The HAPA, covering 12 % of the firn area, contains 8 ± 1 % of GrIS FAC, both for the top 10 m and the whole firn column.
The LAPA, which comprises 14 % of the firn area, contained 9 ± 1 % of ice-sheet-wide firn air content in the period 2010-2017. Decreasing FAC 10 between 1998-2008 and 2010-2017 yields a loss of 170 ± 120 km 3 (23 ± 16 %) of air from the top 10 m of firn. The corresponding decrease in FAC tot indicates that, as an upper estimate, 700 ± 490 km 3 of air may have been lost from the total firn column. In this we assume that the FAC 10 decrease propagated to the entire firn column (see Sect. 2.5), which might not be accurate. Insufficient data are available to determine precisely how much FAC was lost below 10 m and we can only give a hypothetical upper bound to the FAC tot decrease.
Recent studies have identified increasing surface melt and meltwater refreezing as major contributors to increasing near-surface firn densities and subsequent loss of FAC (de la Peña et al., 2015;Charalampidis et al., 2015;Machguth et al., 2016;Graeter et al., 2018). However, firn density and FAC are also dependent on annual snowfall, with decreasing snowfall driving increasing firn density and decreasing FAC (e.g. Vandecrux et al., 2018). Nevertheless, the lack of widely distributed observations of snow accumulation for the 1998-2017 period and the contradicting trends in precipitation calculated by RCMs (Lucas-Picher et al., 2012;van den Broeke et al., 2016;Fettweis et al., 2017) complicate the partitioning of the melt and snowfall contributions to changes in GrIS FAC.
To investigate how uncertainties in T a andċ impact our FAC 10 maps, we repeat our procedure using the 1979-2014 T a andċ estimated by Box (2013) and Box et al. (2013) (hereafter referred to as "Box13"). The Box13-derived FAC 10 fits equally well (RMSD < 0.3 m) to the FAC 10 observations, leading to spatially integrated FAC values within the uncertainty of the MAR-derived values. However, due to differwww.the-cryosphere.net/13/845/2019/ The Cryosphere, 13, 845-859, 2019  ing model formulations and atmospheric forcings, the spatial patterns of air temperature and snowfall are different between Box13 and MARv3.5.2 (detailed in Fettweis et al., 2017), especially in the southern and eastern regions of the firn area. This leads to different estimations of FAC 10 in these regions (Fig. S4). Additionally, in these regions no firn observations are available to constrain our FAC 10 estimates. More observations in the sparsely observed southern and eastern regions would improve FAC 10 estimates and help better elucidate which T a andċ source best describes the spatial pattern in FAC 10 .

Firn retention capacity
The decrease in FAC 10 in the LAPA between 1998-2008 and 2010-2017 translates to a loss in meltwater retention capacity of 150 ± 100 Gt in the top 10 m of firn (Table 3). This lost The Cryosphere, 13, 845-859, 2019 www.the-cryosphere.net/13/845/2019/ retention capacity represents 0.4 ± 0.3 mm sea-level equivalent (s.l.e.). For the total firn column, we estimate an associated upper bound loss of 540 ± 440 Gt (1.5 ± 1.2 mm s.l.e.). While these volumes are small compared to the average GrIS mass loss (∼ 0.47±0.23 mm s.l.e. yr −1 for 1991-2015 in van den Broeke, 2016), the impact of reduced retention capacity has an important time-integrated effect in amplifying meltwater runoff each year. This amplification can be non-linear as when, for instance, a succession of anomalously high melt years and reduced firn permeability resulted in an abrupt increase in western Greenland runoff in 2012 (Machguth et al., 2016). Harper et al. (2012), using observations from 2007 to 2009, estimated that 150 000 km 2 of firn residing within the lower percolation area (as delineated in an earlier version of MAR) could potentially store between 322 ± 44 Gt of meltwater in the top 10 m of firn and 1289 +388 −252 Gt within the entire firn column. We note that the Harper et al. (2012) estimate is based solely on observations in the LAPA, while 68 % of the percolation area to which they extrapolate is located in the HAPA. By contrast, we find that the warmest 150 000 km 2 of our firn area in 2010-2017 can retain only 150 ± 66 Gt of meltwater in the top 10 m of the firn. We estimate a total storage capacity of 310±270 Gt within the whole firn column in this part of the firn area. Our relatively low estimate of the retention capacity might reflect the recent decrease of FAC in the LAPA but also, for the values derived from FAC tot , our simplifying assumption that this decrease has propagated through the whole firn column (Sect. 2.5). Yet, beyond these integrated values, our approach allows the quantification of the firn retention capacity and the corresponding uncertainty at any location of the firn area. Our product can therefore be used in combination with, for instance, remotely sensed melt extent to derive which areas of the firn actively retain meltwater and evaluate the retention capacity there.
We use the same infiltration ice density as , 843±36 kg m −3 as determined from firn core segments saturated by refrozen meltwater. However, Machguth et al. (2016) measured an infiltration ice density of 873 ± 25 kg m −3 with a similar technique in western Greenland. Using the latter value increases our estimated firn storage capacity of the top 10 m of firn by 8 % to 13 %, depending on the region, but remains within our uncertainty intervals (Table 3). Additional field measurements are needed to ascertain the spatial and temporal dependence of infiltration ice density on climatic drivers. Our definition of retention capacity assumes that retention occurs through the refreezing of meltwater and neglects potential liquid water retention seen in firn aquifers . Nevertheless, recent work in southeastern Greenland showed that meltwater resides for less than 30 years in the aquifer before it flows into nearby crevasses and eventually leaves the GrIS (Miller et al., 2018). Meltwater refrozen within the firn can be retained for much longer periods, until it is discharged at a marine-terminating outlet glacier or reaches the surface of the ablation area. By neglecting liquid water retention in firn, our study focuses on long-term meltwater retention.
3.3 Regional climate model evaluation

Comparison with FAC observations
All models reproduce the FAC 10 observations in the DSA and HAPA with a bias ≤ 0.2 m and RMSD ≤ 0.4 m (Fig. 6,  Table 4). RACMO2.3p2, MARv3.9 and HH_LIN tend to underestimate the FAC 10 in the LAPA, while HH_MOD does not show a pronounced bias there. The RCMs all present a RMSD less than 12 % of the mean FAC 10 for our entire dataset. The RCMs are also evaluated against the 29 directly observed FAC tot (Fig. 6, Table 4). Both versions of HIRHAM5 overestimate FAC tot in the DSA (bias > 3 m), while RACMO2.3p2 performs better in that area (bias = 0.1, RMSD = 1.8). HH_LIN and RACMO2.3p2 compare relatively well with the three FAC tot observations available in the LAPA, while HH_MOD presents a larger positive bias. These three FAC tot observations are located in the upper LAPA and therefore do not include regions where RCMs underestimate FAC 10 . All models overestimate the only FAC tot observation available in the HAPA by more than 3 m. Compared to all FAC tot measurements, RACMO2.3p2 gives a RMSD equivalent to 9 % of the mean observed FAC tot when HIRHAM5's RMSD reaches 20 % with HH_MOD. None of the RCMs therefore simulate both FAC 10 and FAC tot accurately.

Comparison with the spatially integrated FAC
Agreement between RCM-simulated and observationderived spatially integrated FAC is model-and regiondependent (Fig. 7). RCMs simulate a spatially integrated FAC 10 within the uncertainty of our observation-derived estimation in the DSA. Models also show lower spatially integrated FAC 10 in the LAPA and higher values in the HAPA compared to our estimate ( Fig. 7b-d). These regional differences cancel out when spatially integrating FAC 10 over the entire firn area (Fig. 7a). Our estimation of spatially integrated FAC tot is subject to more assumptions, as uncertainty is introduced in our conversion of FAC 10 to FAC tot (Sect. 2.5). In the DSA, HH_MOD simulates a spatially integrated FAC tot 20 % higher than our estimation while RACMO2.3p2 simulates spatially integrated FAC tot within our uncertainty range (Fig. 7e). In the LAPA, the decrease in spatially integrated FAC tot is more pronounced in our estimate than in the RCMs. This might indicate that, in the RCMs, the FAC loss is concentrated in the near-surface firn and has not yet propagated through the entire firn column. Our estimate assumes that any change in FAC 10 immediately propagates to the entire firn pack (see Sect. 2.5   ues than our estimate (Fig. 7h), contributing to the higher spatially integrated FAC tot across the entire firn area in the RCMs compared to our estimation (Fig. 7e). This is partly due to the fact that in our estimation, FAC decreases with elevation and is set to zero at the firn line. In the RCMs, modelled FAC remains higher than our estimate in the lower HAPA and in the vicinity of the firn line. No FAC observations are available in the lower HAPA to confirm this. Future measurements will help to quantify FAC in the surroundings of the firn line, allowing for better evaluations of our assumptions and further assess to the RCMs' performance in that area.
The differences between RCM outputs may stem from their respective surface forcings. As an illustration, HH_MOD uses a higher albedo than HH_LIN and thus calculates less surface melt and refreezing and, as a conse-quence, higher FAC 10 in the LAPA. Noël et al. (2018) found that the surface mass balance of RACMO2.3p2 in the accumulation area was on average slightly lower than observations, indicating excessive sublimation or runoff relative to snowfall in the model. This surface bias could explain the model's underestimation of FAC 10 in the LAPA at point scale (Fig. 6, Table 4) and on spatially integrated values (Fig. 7). On the other hand, MARv3.9 has slight positive biases in surface mass balance compared to observations (Fettweis et al., 2017). And although this RCM simulates too much precipitation relative to melt, it also underestimates FAC 10 in the LAPA. Surface forcing is therefore not the only factor influencing the FAC estimates by the RCMs.
Differences in RCM-simulated FAC 10 can also be explained by the way firn densification is treated in the snow model of each RCM. For instance, the overestimation of FAC tot in the DSA by HIRHAM5 potentially arises from the use of a firn compaction law originally developed for seasonal snow (Vionnet et al., 2012). RACMO2.3p2 produces more realistic FAC tot in the DSA, most likely because the densification law it uses has been tuned to match eight deep firn density observations (Kuipers Munneke et al., 2015a). It is nevertheless difficult to disentangle the roles of surface forcing and model formulation in the performance of RCMs.
In agreement with our observation-derived FAC 10 estimates, the RCMs calculate a decreasing FAC 10 in the LAPA (Fig. 7c) initiating in the early 2000s and accelerated during the extreme summers of 2010 and 2012. In the DSA, RCMs show a FAC 10 decrease ranging from −120 km 3 in MARv3.9 to −282 km 3 in RACMO2.3p2 between 1998 and 2017. These decreases contradict our conclusion that FAC has not changed significantly in the DSA over that period (Sect. 3.1). The different FAC 10 dynamics in our dataset and in RCMs could be due to (i) the RCMs not capturing an increase of snowfall in the DSA, which could in theory counterbalance the densification expected from the recent warming in the firn area (McGrath et al., 2014); (ii) an overestimated response of firn compaction rates to increasing temperatures in the models; or (iii) the spatial heterogeneity and uncertainty of FAC observations leading to spurious conclusions from our dataset. Yet, finding identical firn density profiles decades apart at several sites (e.g. Summit Station, Camp Century) adds confidence to our findings.

Conclusions
Using a collection of 360 firn density profiles spanning 65 years we quantified the firn air content (FAC) of the Greenland ice sheet as function of long-term air temperature and net snow accumulation averages (T a andċ). For the 2010-2017 period we calculate that the Greenland firn contained 26 800 ± 1840 km 3 of air, of which 6500 ± 450 km 3 is in its top 10 m. We find that over the 1953-2017 period, FAC remained constant within uncertainty in the dry snow area (DSA, where T a ≤ −19 • C). We note that the vast majority of the ice sheet's FAC (83 ± 5 %) resides within the DSA and represents a potential meltwater storage volume of 12 800 ± 1170 Gt. In the low-accumulation percolation area (LAPA, where T a > −19 • C andċ ≤ 600 mm w.e. yr −1 ), we calculate that the FAC decreased by 23±16 % between 1998-2008 and 2010-2017. This decrease translates into the loss of meltwater retention capacity of 150 ± 100 Gt (0.4 ± 0.3 mm sea-level equivalent) in the top 10 m of the firn and potentially up to 540 ± 440 Gt (1.5 ± 1.2 mm sea-level equivalent) in the entire vertical extent of the firn layer. This decreased FAC and meltwater retention capacity is focused in the lower-accumulation area of central western Greenland. Thus, in contrast to the relative stability of the DSA, the LAPA is the focal area of the firn's response to recent climate change. The firn in the high-accumulation percolation area (LAPA, where T a > −19 • C andċ > 600 mm w.e. yr −1 ) has the capacity to store 370 ± 70 Gt in its top 10 m and up to 960 ± 290 Gt in its entire vertical extent. Yet, this area is covered by fewer observations and would highly benefit from future field surveys.
The outputs from three regional climate models (HIRHAM5, RACMO2.3p2 and MARv3.9) indicate that our calculated decrease in FAC may have been initiated in the early 2000s and accelerated after 2010. The RCMs also provide estimates of FAC in regions where no measurements are available. But the mismatch between RCMs and www.the-cryosphere.net/13/845/2019/ The Cryosphere, 13, 845-859, 2019 our firn core dataset illustrates that RCMs should be used with caution when assessing meltwater retention capacity, or when converting ice-sheet volume changes into mass changes in the firn area. Finally, our study highlights the importance of assimilating in situ firn density measurements to document the climate response of ice-sheet firn as a nontrivial component of the sea-level budget. More broadly, this work illustrates how new insight can be obtained from the synthesis of historical data sources and thus emphasises the tremendous value of open-access data within the scientific community.
Author contributions. BV, MM, HM, WTC, DvA, AH, CMS, CC, EMM, EMT, LK, LNM, CM and SBS collected firn density profiles either in the field or from previously published sources. BV, RSF, JEB and TIN designed the study. BV wrote the manuscript, to which all co-authors contributed.
Competing interests. The authors declare that they have no conflict of interest.