Improved representation of the contemporary Greenland ice sheet ﬁrn layer by IMAU-FDM v1.2G

. The ﬁrn layer that covers 90% of the Greenland ice sheet (GrIS) plays an important role in determining the response of the ice sheet to climate change. Meltwater can percolate into the ﬁrn layer and refreeze at greater depths, thereby temporarily preventing mass loss. However, as global warming leads to increasing surface melt, more surface melt may refreeze in the ﬁrn layer, thereby reducing the capacity to buffer subsequent episodes of melt. This can lead to a tipping point in meltwater 5 runoff. It is therefore important to study the evolution of the Greenland ﬁrn layer in the past, present and future. In this study, we present the latest version of our ﬁrn model, IMAU-FDM (Firn Densiﬁcation Model) v1.2G, with an application to the GrIS. We improved the density of freshly fallen snow, the dry-snow densiﬁcation rate and the ﬁrn’s thermal conductivity using recently published parameterizations and by calibrating to an extended set of observations of ﬁrn density, temperature and liquid water content at the GrIS. Overall, the updated model settings lead to higher ﬁrn air content and higher 10 m ﬁrn temperatures, 10 owing to a lower density near the surface. The effect of the new model settings on the surface elevation change is investigated through three case studies located at Summit, KAN-U and FA-13. Most notably, the updated model shows greater inter-and intra-annual variability in elevation and an increased sensitivity to climate forcing.


Observations
IMAU-FDM output is evaluated using previously available and newly obtained profiles of firn density, temperature and liquid water content from the GrIS accumulation zone. The observations are from 128 different locations to ensure that the various ice facies and climate zones of the GrIS are well represented (Fig. 1). Vertical profiles of observed firn density from ice cores 65 vary in depth from 9.6 to 150.8 m and have been drilled between 1952 and 2018 in the framework of the Program for Arctic Regional Climate Assessment (PARCA; McConnell et al. (2000); Mosley-Thompson et al. (2001); Hanna et al. (2006); Banta and McConnell (2007)), the Arctic Circle Traverses (ACT, Box et al. (2013)) and the EGIG line (Harper et al. (2012b)), Das 1 and Das 2 (e.g. from Hanna et al. (2006)) and several other cores were retrieved from the SUMup database (SUrface Mass balance and snow depth on sea ice working grouP), (Koenig et al. (2013); Koenig and Montgomery (2019)). 70 Temperature observations include profiles ranging in depth between 4 and 14 m obtained by Harper et al. (2012b) along a transect in the western GrIS and at the NEEM deep ice core drilling site ). Additional firn temperature observations are from Summit, Dye-2 (Vandecrux et al. (2019a), KAN-U (Charalampidis et al. (2015)) and FA-13 ). An additional 14 observations of 10 m firn temperatures are from Polashenski et al. (2014).
For observations of liquid water in firn, we use observations from Dye-2 (Heilig et al. (2018)), obtained using an upward-75 looking ground-penetrating radar (upGPR), which was installed and operated in the summer of 2016. The upGPR was buried ∼ 4.5 m under the snow, and was capable of measuring the liquid water percolation depth, content as well as the changing distance between the instrument and the snow surface. Although the data do not cover a wide spatial (single location) or temporal range (1 May-16 October 2016), they are unique and moreover have high temporal and vertical resolution, making them very valuable for firn model evaluation ), but also to evaluate melt intensity and timing in the 80 forcing time series.

IMAU-FDM
For this work we use the offline IMAU-FDM, a semi-empirical firn densification model that simulates the time evolution of firn density, temperature, liquid water content and changes in surface elevation owing to variability of firn depth. The model has been compared extensively to, and calibrated with observations of firn density and temperature from the ice sheets of 85 Greenland and Antarctica (Ligtenberg et al. (2011);Kuipers Munneke et al. (2015a); Ligtenberg et al. (2018)). IMAU-FDM is forced by three-hourly output of the polar version of the Regional Atmospheric Climate Model (RACMO2.3p2) ). Over glaciated grid cells, the RACMO2 subsurface model uses approximately the same expressions as IMAU-FDM, but with a lower vertical resolution (max. 150 vs. 3000 layers) and less comprehensive initialization to save computation time.
In the current version of IMAU-FDM we do not consider the subsurface penetration of shortwave radiation (Van Dalum et al. 90 (2020)). For both the ice sheets of Greenland and Antarctica, the performance of IMAU-FDM has been comparable to the more physically-based SNOWPACK model (Steger et al. (2017b); Van Wessem et al. (2021); Keenan et al. (2021)). In the following An important boundary condition for the model is the density of freshly fallen snow, ρ 0 . When determined from field observations, fresh snow density is often assumed equal to the near-surface density, loosely defined as the average density of the top 0.5 m of dry snow. As density is highly variable near the surface, the exact chosen depth is critical for the outcome, which hampers a robust comparison between datasets (Fausto et al. (2018)). In firn models, fresh snow density is commonly parameterized as a function of meteorological variables such as temperature and wind speed at the time of deposition, or, when these are not 100 available, using annual average values instead (Keenan et al. (2021)). Several studies have addressed the parametrization of ρ 0 on the GrIS (Kuipers Munneke et al. (2015a);Fausto et al. (2018)). Assuming a linear dependence of the density on mean annual surface temperature T s , this parametrization takes on the following form: In the updated model, a new parameterization for fresh snow density (Fausto et al. (2018)) was adopted. In contrast to previous studies, which typically use the first ∼ 0.5 − 1.0 m of snow, Fausto et al. (2018) used only the upper 0.1 m of snow to define surface density at 200 locations and found: with T a the annual mean near-surface (usually 2 m) air temperature in • C.
Previously, the climatological mean 2 m air temperature has been used in IMAU-FDM (Kuipers Munneke et al. (2015b)), or an instantaneous value (Ligtenberg et al. (2018)). Using a climatological mean value suppresses the year-to-year variability in snow density. This is undesirable, especially because the model will also be used future scenarios, in which long term trends 115 in temperature may have an effect. On the other hand, using instantaneous temperature values may introduce an excessive variability which, in reality, is smoothened by the effects of the snow being subjected to settling by wind and metamorphosis through numerous, daily warming and cooling cycles. As a trade-off, T a in the updated model is calculated as the average 2 m air temperature of the year preceding the snowfall.

Dry snow densification rate
IMAU-FDM is a 1D, vertical Lagrangian model. When new snow accumulates at the surface (model top), the model layers 125 are buried deeper and tracked during their downward motion. At every time step, each layer is compacted under the influence of the pressure exerted by the mass of snow/firn above it. However, in IMAU-FDM the densification rate dρ dt is not directly related to the overburden pressure, but rather follows a semi-empirical, temperature-dependent equation based on Arthern et al. (2010): whereḃ is the annual average accumulation rate (mm w.e. per year) over the spinup-period (1960-1979), ρ i = 917 kg m −3 is the density of glacial ice, g, E c , E g and R are constants, and T is the instantaneous layer temperature in Kelvin. The average annual accumulation rateḃ is provided by RACMO2 as the amount of total precipitation minus sublimation and drifting snow erosion. Different values of C above and below ρ = 550 kg m −3 represent a shift in the dominant densification mechanism (Cuffey and Paterson (2011)). For ρ < 550 kg m −3 , the densification of the firn is dominated by the settling and sliding of 135 grains. For ρ ≥ 550 kg m −3 recrystallisation, deformation and sublimation become dominant and the densification rate is lower, which is reflected in a lower value for C.
Compared to observations of the depth of the 550 and 830 kg m −3 density levels, a structural bias is found in Eq. 3, that in previous studies turned out to depend on the annual average accumulation rate. In order to calibrate Equation 3 to the 2 × 10 2 3 × 10 2 4 × 10 2 6 × 10 2 Mean accumulation (mm w.e. yr 1 )   (2011)).
Firn densification owing to horizontal compression is neglected, although in fast-flowing regions this can be locally important (Horlings et al. (2021)).
In the model update, we recalibrated the dry densification correction factor MO as a function of mean annual accumulation, by using an updated, high-resolution GrIS accumulation field (Noël et al. (2015)) and optimizing the modelled depths at which Since MO corrects for the dry compaction rate, only dry firn cores (i.e. with little surface melt) are used. A core is considered as "dry" if the mean annual melt is less than 5% of the mean annual accumulation rate.. Least squares fitting yields R 2 values Table 1. Values of the old and new linear regression of Eq. 4, their R 2 as well as the standard error in of the new fitting parameters.  Table 1.
With the update and the use of new firn and accumulation data, the linear relation between MO 550 and ln(ḃ) vanishes, and reduces to an almost constant value of 0.7. Despite the difference with previous formulations in IMAU-FDM, this is similar to findings by Robin (1958) and Herron and Langway (1980), who also found that below 550 kg m −3 the densification rate correlates almost linearly with accumulation, and that this correlation became non-linear at higher densities. The high 160 correlation for MO 830 also implies that the relation between densification rate dρ dt and accumulation is non-linear above ρ = 550 kg m −3 .

Thermal conductivity
In IMAU-FDM, the vertical temperature distribution and its evolution is obtained by solving the one-dimensional heat transfer in which c is the specific heat capacity of the firn, G the subsurface heat flux, k the thermal conductivity of the firn and L a heat source representing the release of latent heat upon the refreezing of liquid water in the firn or the subsurface absorption of solar radiation. Subsurface penetration of short-wave radiation is neglected in the current model version, which is deemed a reasonable approximation for fine-grained, polar snow. The firn temperature profile is initialized using a spin-up period, 170 see Sect. 2.2.6. Before the spin-up, the firn column is initialised at a constant temperature equal to the annual mean surface temperature during the spin-up period. The lower boundary condition assumes a constant heat flux across the lowest model grid cell, i.e. the deep temperature is allowed to change along with long-term changes in surface temperature or internal heat release. The upper boundary condition for the temperature calculation is provided by the surface ('skin') temperature in RACMO2, obtained by iteratively solving the surface energy balance (Van den ). Subsequently, for 175 computational efficiency Eq. 5 is solved using an implicit/explicit scheme in the absence/presence of liquid water (Helsen et al. (2008)). Due to the Lagrangian character of the model, vertical heat advection is implicitly considered (Helsen et al. (2008)).
Any heat generated by firn horizontal/vertical deformation is neglected.
The thermal conductivity is assumed to depend on firn density and temperature, and in previous versions of IMAU-FDM followed the expression for seasonal snow due to Anderson (1976), which only depends on density. In the updated model, the The equation consists of two parts: one for snow and low-density firn and one for ice and high-density firn. The transition between the two regimes remains smooth through the weight factor θ(ρ). The definition of θ and the thermal conductivities that are used in Eq. 6 are: Here k a represents the thermal conductivity of air, taken from Reid et al. (1959). Figure 4 compares the old and new 195 expressions for various temperatures. As can be seen, the new expression takes on a slightly lower value than Anderson (1976) at densities below ∼ 475−565 kg m −3 , depending on the temperature, but a higher value at densities above that. This difference becomes larger at lower temperatures.

Meltwater percolation, retention and refreezing
IMAU-FDM employs a tipping bucket method to treat the percolation, irreducible (capillary) retention and (re)freezing of 200 water, by filling up subsequent deeper layers to maximum capacity in a single model time step (i.e. quasi-instantaneous).
Magnusson and others (Magnusson et al. (2015)) show that, in spite of its simplicity and shortcomings, the tipping bucket method is a robust and useful method to deal with liquid water transport in the snowpack when compared to more sophisticated methods, especially when capturing general firn properties at the larger, multi-kilometre horizontal scale for which IMAU-FDM is designed. In IMAU-FDM, the fraction that is retained in a model layer by capillary forces ('irreducible water content') 205 depends on the available pore space according to the expression by Coléou and Lesaffre (1998): where P is the porosity of the firn layer, defined as P = 1−ρ/ρ i . The maximum amount of water that is stored thus decreases with increasing density of the firn layer. Standing water and lateral runoff over ice-layers are currently ignored; the latter is  considered a fair assumption, because on the spatial scales at which the model is employed (i.e. the RACMO2 grid of 5.5 by 210 5.5 km) it is assumed that within a model grid cell the meltwater can usually find a way to flow around a layer of ice. Varying the irreducible water content by, e.g., multiplying Eq. 7 with a constant factor or using a constant volume or mass fraction, did not improve the result, and it was decided to leave the liquid water scheme unchanged.

Firn thickness and elevation change
IMAU-FDM tracks the total firn thickness, and changes in it. The resulting vertical velocity of the ice-sheet surface due to 215 changes in the firn layer ( dh dt ) is given by: The total vertical surface velocity dh dt can thus be decomposed into separate contributions from accumulation (v snow ), surface sublimation (included in v snow ), sublimation by snowdrift ( snd ), erosion or deposition by snowdrift (v er ), snowmelt (v melt ), and firn compaction (v f c ). The term v ice is defined as the mean surface mass balance (SMB, v ice = v snow +v snd +v er +v melt ) 220 with an opposite sign. It represents the long-term average vertical mass flux through the lower boundary of the firn column, which equals the mass flux through the upper boundary in a steady-state firn layer. In Sect. 4 we show surface elevation change and the individual components for three case study locations on the GrIS.

Model initialisation
The latest IMAU-FDM model runs span the period 1 January 1960 -31 December 2020. The initial model density, temperature At the upper boundary of IMAU-FDM, mass accumulation (solid precipitation minus sublimation minus drifting snow erosion), liquid water fluxes (melt plus rainfall minus evaporation) and surface temperature are prescribed from the regional atmospheric climate model RACMO2.3p2, which has been used to simulate the climate and surface mass balance of the GrIS and its immediate surroundings for the period 1958-2020 at a horizontal resolution of 5.5 km. This version of RACMO2 has been extensively evaluated over the GrIS (Noël et al. (2018b)). At the lateral boundaries, using a relaxation zone of 24 240 gridpoints, RACMO2 is forced by European Centre for Medium-Range Weather Forecasts (ECMWF) re-analysis data, i.e. ERA-40 between 1958and 1978, ERA-Interim between 1979and 1990and ERA-5 between 1991. For the forcing of IMAU-FDM the full spatial resolution of 5.5 km is used and a temporal resolution of 3 hours was selected, as an acceptable trade-off between robustly resolving the daily cycle and keeping manageable file sizes. IMAU-FDM typically uses a timestep of 3 min (explicit temperature calculation scheme) to 3 h (implicit temperature calculation scheme), for which we linearly 245 interpolate the forcing between the RACMO2 forcing time steps. compare modelled and observed vertically integrated firn air content (FAC), i.e. the vertical distance over which the firn layer can be compressed until reaching the density of glacier ice across the entire firn column.
Here, n z is the number of layers in that firn profile, ∆z j is the thickness of layer j and ρ j is the density of that layer. Note 255 that here the FAC is calculated over the complete depth at which observations are available, as opposed to FAC 10 shown in  With the newly adopted parametrizations, the simulation of FAC in dry locations has significantly improved (Fig. 5  and without temperature dependency for surface snow density. These results show that including temperature as a predictor for the surface density does improve model performance. Thus, we opt for the temperature dependent formulation, since this also 275 allows capturing the effect of long-term temperature trends on the surface density. more realistic shape and reduced variability. It increases the pore space and brings simulated FAC in better agreement with the 280 observed density profile. One of the main reason for the increased performance is the change from an instantaneous surface density parameterization to one that is based on annual mean values. This leads to greatly reduced "peaks" in the density profile, which is much more in line with observations since the surface density. For FA-13, it also seems that the lower surface density matches the upper 25 m of the density profile better.

Firn temperature 285
Modelled and measured 10 m firn temperatures at 31 locations are compared in Fig. 7. The new settings improve results, especially for the warmer locations with significant melt, which are mostly locations from Harper et al. (2012a) in west Greenland.
Here, the cold bias has been significantly reduced; for locations with T 10 > −20 • C, the mean RMSE decreased from 4.7 to 2.7 • C, respectively (-43%).
The main reason for the improvement is a better representation of the density at those locations, which allows for more 290 realistic refreezing and the associated enhanced latent heat release, increasing the temperature in melt-prone locations. Additionally, the lower conductivity due to the lower density leads to a less efficient cooling of the deeper snow during winter. In spite of the clear improvement, a cold-bias remains for some of these locations, which can be partly attributed to a cold-bias in the RACMO2 forcing. For the low-melt locations (T 10 < −20 • C), a persistent warm model bias remains. is warm and wet, has a cold bias. For both locations the new surface density parameterization has decreased the density in the upper layers. This in turn leads to a lower thermal conductivity since the thermal conductivity increases monotonically with density (Fig. 4). Therefore, heat diffuses slower in the upper layers, increasing the temperature gradient there, which can be seen clearly at Summit. For both locations the depth at which the thermal maximum occurs also decreased slightly.

300
Dye-2 now clearly shows a temperature maximum at the depth at which refreezing occurs in contrast to the old model settings. This is also attributed to a decrease in the thermal conductivity: previously, heat generated by refreezing was able to more rapidly diffuse to greater depths or the atmosphere, but now it remains "trapped" around the depth at which refreezing occurs. Another factor that contributes is that refreezing occurs at a greater depth than before, see Sect. 3.3.

305
The liquid water percolation and retention schemes have not been updated, but the changes made to the parameterizations that impact density and temperature do influence water percolation, and therewith liquid water content (LWC), and these changes are discussed here. Very few in-situ, vertically resolved observations of LWC are available. Here we used data from a recent study that used upward looking ground penetrating radar (upGPR) at Dye-2 in the higher percolation zone of the southwestern GrIS (2120 ma.s.l., see Fig. 1, Heilig et al. (2020)). gradients at Dye-2. The increase in temperature means that the water needs to percolate deeper into the firn pack before it can refreeze, which is reflected in the increased penetration depth. Simultaneously, the decrease of the surface density means that there is more pore space near the surface that can retain irreducable water, explaining the increase in volume fraction.
Overall, the penetration depth now agrees better with the observations, although in IMAU-FDM, the meltwater still refreezes too quickly in more shallow layers than observed. 320 4 Pilot application to firn-induced surface elevation change In this section we compare time series (1958-2020) of firn-induced surface elevation (i.e. firn depth) changes at three key locations: Summit in the cold and dry ice sheet interior, KAN-U in the relatively warm and dry southwestern percolation zone and FA-13 in the wet and relatively mild southeastern firn aquifer region Forster et al. (2014), as indicated by the green circles in Fig. 1). Table 2 provides geographical and climatological information of these locations. The 325 three locations represent three very different climates and are therefore useful for investigating how the new model settings affect firn depth change, which represents the surface elevation change in the absence of contributions from ice dynamics, basal melt and/or bedrock elevation change.

Summit
Summit is located at the centre of the GrIS at a high elevation and therefore it experiences a low amount of snowfall and a neg-330 ligible amount of rain and melt. The evolution of its elevation is therefore closely linked to changes in the temperature (higher temperatures lead to a higher compaction rate) and accumulation (higher accumulation leads to a higher surface elevation).   Table 2. Location and climate climate of the three case study sites. The annual mean accumulation are calculated over the whole simulation period .

Lon.
Lat. Elevation T2 m Acc. Melt and v f c act on different timescales. v f c is fairly constant in time and changes in tandem with the seasonal changes in temperature. v snow on the other hand is much more variable, as snowfall is very episodic and highly variable at many time scales. As a result, the surface elevation increases more during a snowfall event and decreases during dry episodes, and these ratios have increased in the new model, leading to larger interannual variations. This also implies that the firn model is now more sensitive to changes in the forcing, reacting more strongly to a decrease or increase in accumulation and skin temperature in the future.

KAN-U
Situated in the southwest and at a lower elevation, KAN-U is warmer than Summit, but most importantly, melting occurs every year during the summer, which greatly affects the firn properties at its location (Fig. 11). The average influence of surface melt When comparing the vertical surface velocities between the old and the new model settings, we again see an increased accumulation velocity v snow , due to a lower fresh snow density. This is again compensated for by a more negative compaction velocity, resulting in a very similar net velocity (v tot ). As accumulation decreased after 2005, the net effect on the surface is a slight lowering. Following significant warming and increased melt at this site (Fig. 11), the contribution of v melt to firn 370 depth changes increases and that of v f c decreases, making the former the dominant process leading to surface lowering at KAN-U. The strong increase in melt causes a larger downward velocity of the surface, leading to thinning. v melt is also larger in magnitude with the new settings because the melted snow at the surface is of a lower density.
Just like at Summit, the elevation change is more sensitive to its forcing with the new model settings than previously was the case. This is especially apparent in the years 1983 and 2012. Throughout the time series, interannual variability in firn depth 375 is first dominated by snowfall (v snow ), but towards the end of the time series the contribution to the total variability made by v melt increases rapidly.

FA-13
At the two previous sites the new model settings produce similar results for the overall surface elevation change. However, at the third site of FA-13, a site with a firn aquifer, we found that the new model settings produce a significantly different 380 result compared to the old ones. This location experiences a warmer and wetter climate than KAN-U, which leads to a rapid densification in the upper part of the firn column (Fig. 6).
Here, the signal is dominated by large seasonal oscillations in firn depth of up to ∼ 1 m yr −1 between 1960 and 1985. From 1985 onwards, the firn depth decreases until 2012, but at a higher rate in the updated than in the previous model (∼ 0.35 vs.

Summary and outlook
Temporal and spatial variability in firn thickness is highly relevant for the mass balance of the Greenland ice sheet (GrIS), 400 because it directly impacts its refreezing efficiency. Moreover, firn thickness change is an important component of surface elevation change, and improved knowledge is required to accurately convert remotely sensed GrIS volume to mass changes.
In this paper, we presented improvements in the offline version of the IMAU firn densification model (IMAU-FDM v1.2G), forced by three-hourly output of the regional climate model RACMO2.3p2. Taking advantage of improved climate forcing and newly available observations of surface and subsurface firn density and temperature, the improvements are systematically 405 implemented in the parametrizations of surface density, dry snow densification and thermal conductivity. The treatment of liquid water is not changed, owing to a lack of sufficient observations to justify changes in the current configuration.
The updated model predicts predicts higher firn air content (FAC), which at three selected sites in the interior GrIS and in the southwestern and southeastern percolation zone results in a larger sensitivity of firn thickness to intra-and interannual variations in snowfall, melt and temperature. As an important consequence of a change in fresh snow density parameterization, 410 the inter-and intra-annual variations in elevation have increased, owing to an increased sensitivity to changes in its forcing. In a warmer climate, firn thinning owing to increased surface melt becomes increasingly important at the marginal sites, both in the mean and as a component of interannual variability. Future applications of the improved model include a full GrIS assessment of contemporary and future firn mass and thickness changes, as well as explaining areas where firn aquifers and ice slabs currently occur, and their future changes.
Author contributions. MB, PKM, WJvdB and MRvdB started this project, decided on its scope, which parts of the model required further development and interpreted the results. MB performed the model simulations, implemented the changes to the model, comparisons and led the writing of the manuscript. All authors contributed to discussions on the manuscript.