Development of physically based liquid water schemes for Greenland firn-densification models

As surface melt is increasing on the Greenland Ice Sheet (GrIS), quantifying the retention capacity of the firn layer is critical to linking meltwater production to meltwater runoff. Firn-densification models have so far relied on empirical approaches to account for the percolation–refreezing process, and more physically based representations of liquid water flow might bring improvements to model performance. Here we implement three types of water percolation schemes into the Community Firn Model: the bucket approach, the Richards equation in a single domain and the Richards equation in a dual domain, which accounts for partitioning between matrix and fast preferential flow. We investigate their impact on firn densification at four locations on the GrIS and compare model results with observations. We find that for all of the flow schemes, significant discrepancies remain with respect to observed firn density, particularly the density variability in depth, and that inter-model differences are large (porosity of the upper 15 m firn varies by up to 47 %). The simple bucket scheme is as efficient in replicating observed density profiles as the single-domain Richards equation, and the most physically detailed dual-domain scheme does not necessarily reach best agreement with observed data. However, we find that the implementation of preferential flow simulates ice-layer formation more reliably and allows for deeper percolation. We also find that the firn model is more sensitive to the choice of densification scheme than to the choice of water percolation scheme. The disagreements with observations and the spread in model results demonstrate that progress towards an accurate description of water flow in firn is necessary. The numerous uncertainties about firn structure (e.g. grain size and shape, presence of ice layers) and about its hydraulic properties, as well as the one-dimensionality of firn models, render the implementation of physically based percolation schemes difficult. Additionally, the performance of firn models is still affected by the various effects affecting the densification process such as microstructural effects, wet snow metamorphism and temperature sensitivity when meltwater is present.

1 Introduction

Estimating the properties of the firn layer – and how it evolves under a warming climate – is a critical step in measuring the ice sheets' contribution to sea level rise, yet it remains one of the key sources of uncertainty in present assessments (McMillan et al., 2016). Accurate estimates of firn thickness and density are required for the conversion of space-borne measurements of volume change into mass change (e.g. McMillan et al., 2016; Shepherd et al., 2018). Also, assessments of the Greenland and Antarctic ice sheets' contribution to sea level require information on firn density and spatial distribution in order to calculate meltwater retention potential and the capacity of firn to buffer the flow of meltwater to the ocean (Harper et al., 2012; Machguth et al., 2016; van den Broeke et al., 2016). Surface melting has become more widespread and intense on the Greenland Ice Sheet (GrIS), with annual total melt rates rising by 11.4 Gt yr⁻² between 1991 and 2015 (van Angelen et al., 2014; van den Broeke et al., 2016). This meltwater percolates into the firn layer, where it can refreeze, run off or remain liquid in temperate firn. Refreezing of liquid water in firn, known as internal accumulation, has an impact both on ice sheet mass balance and on heat fluxes from the surface to the ice sheet (van Pelt et al., 2012). As such, understanding physical processes in firn, including in particular the transport of liquid water, is becoming increasingly important in order to accurately constrain and predict the mass balance of the GrIS (van den Broeke et al., 2016).

Densification of dry firn is typically modelled as a function of near-surface air temperature and accumulation (Herron and Langway, 1980; Arthern et al., 2010; Li and Zwally, 2011; Simonsen et al., 2013; Morris and Wingham, 2014; Kuipers Munneke et al., 2015a, b). If applied to wet firn, these models are often modified to include a simplified representation of liquid water percolation, the bucket scheme, which assumes flow and refreezing through the firn column occur in a single time step (Simonsen et al., 2013; Kuipers Munneke et al., 2015b; Steger et al., 2017a). Observations have shown that, in reality, liquid water transport in firn is characterised by flow patterns that are heterogeneous in space and time (Pfeffer and Humphrey, 1996; Humphrey et al., 2012). Incorporation of liquid water schemes representing such flow patterns would enable models to better represent the transport of mass and heat through the firn; these schemes might also improve modelled densification in wet firn conditions (Kuipers Munneke et al., 2014; van As et al., 2016a; Meyer and Hewitt, 2017). Liquid water flow, however, is a complex function of several properties and processes that are difficult to constrain by observations and, as a corollary, are difficult to represent in these models (e.g. presence of impermeable ice layers, snow hydraulic properties, grain size, lateral runoff). The infiltration of water through firn can be partitioned between the progressive advance of a uniform wetting front through the pores, called matrix flow, and fast, localised, preferential flow (Waldner et al., 2004; Katsushima et al., 2013). This dual nature of water flow has been reported in observations of the firn layer of the GrIS, where preferential flow pathways come in the form of discrete vertical conduits and are crucial to effectively transport surface meltwater in deep subfreezing firn (Pfeffer and Humphrey, 1996; Parry et al., 2007; Humphrey et al., 2012; Cox et al., 2015). The detection of perennial firn aquifers (PFAs), in which large amounts (140 Gt) of liquid water are stored year-round in deep firn, further emphasises the importance of firn hydrology on the GrIS mass balance (Forster et al., 2014; Koenig et al., 2014). Snow modellers have developed liquid water schemes based on the Richards equation (RE) to simulate matrix flow (Hirashima et al., 2010; Wever et al., 2014; D'Amboise et al., 2017). The RE is a continuity equation describing water flow in unsaturated porous media and is widely used in hydrological models. Recently, a preferential flow scheme has been included in the SNOWPACK model to account for heterogeneous percolation (Wever et al., 2016). Until now, however, such developments have not been implemented in firn-densification models.

In this study, we describe and compare liquid water schemes of different levels of physical complexity from snow models, and we apply these in combination with firn-densification models in order to evaluate the impact of the treatment of liquid water flow on modelled firn densification and temperature. We use the Community Firn Model (CFM) as the modelling framework for our study; the CFM is able to simulate numerous physical processes in firn and includes a large choice of governing formulations for densification. We use the common bucket approach and also develop schemes for liquid water flow in firn following physically based advances in snow models (Wever et al., 2014, 2016; D'Amboise et al., 2017).

We simulate liquid water flow and firn densification starting from 1980 at four sites on the GrIS: DYE-2, NASA-SE, KAN-U and a PFA site (Fig. 1). These sites were chosen because they are collocated with recently drilled firn cores which allow a direct comparison of model results with observations. By comparing simulated firn densification to observations at these sites, we investigate the sensitivity of the system to the choice of liquid water flow scheme and the sensitivity of the flow schemes to various parameterisations of firn structural properties. Finally, we perform simulations with a range of firn-densification formulae and assess the relative importance of the choice of liquid water flow scheme to the choice of the underlying densification equation.

Figure 1 Study site locations and mean annual melt rates (1958–2017) from RACMO2.3p2.

2 Firn model and data

In this study we use and further develop the CFM, an open-source firn-densification modelling framework. We refer the reader to Stevens (2018) for details and briefly summarise the main characteristics here. The CFM is one-dimensional and works in a Lagrangian framework; it is forced at its upper boundary by observed or modelled values for accumulation, surface temperature, surface density, rain and snowmelt. The CFM includes many of the commonly used dry firn-densification schemes (e.g. Herron and Langway, 1980; Helsen et al., 2008; Arthern et al., 2010; Li and Zwally, 2011; Morris and Wingham, 2014; Kuipers Munneke et al., 2015b). We refer the reader to the original publications for details on the different densification schemes and briefly outline the expressions used in our simulations in this section.

2.1 Dry firn-densification model

As a base case, we use the firn-densification formulation implemented in the snow model CROCUS (Vionnet et al., 2012), Eq. (1). It has previously been used in model studies of firn densification on the GrIS and also on polar ice caps (Gascon et al., 2014; Langen et al., 2017). This model is formulated so that densification is based on the overburden stress:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}t}}=\mathit{\rho }{\displaystyle \frac{\mathit{\sigma }}{\mathit{\eta }}},\end{array}$$

where ρ is the density of the firn (kg m⁻³), σ is the stress due to weight of the upper layers (kg m⁻¹ s⁻²) and η is the snow viscosity (kg m⁻¹ s⁻¹) following the parameterisation:

$$\begin{array}{}\text{(2)}& \mathit{\eta }=f_{\mathrm{1}}{f}_{\mathrm{2}}{\mathit{\eta }}_{\mathrm{0}}{\displaystyle \frac{\mathit{\rho }}{c_{\mathit{\eta }}}}\exp \left[a_{\mathit{\eta }}\left(T_{\mathrm{0}}-T\right)+b_{\mathit{\eta }}\mathit{\rho }\right],\end{array}$$

where η₀=7.62237 kg s⁻¹ m⁻¹, a_η=0.1 K⁻¹, b_η=0.023 m³ kg⁻¹, T is the firn temperature (K) and T₀=273.15 K. The parameter c_η is set to 358 kg m⁻³ as suggested by van Kampenhout et al. (2017) when using Eq. (1) for polar firn. There are two additional correction factors, f₁ and f₂, depending on firn microstructural properties. The factor f₁ accounts for the presence of liquid water:

$$\begin{array}{}\text{(3)}& f_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathrm{60}\mathit{\theta }}},\end{array}$$

where θ is the volumetric water content (m³ m⁻³). In this study, we neglect the change in snow viscosity for grain sizes smaller than 0.3 mm by keeping the constant value f₂=4 (after Langen et al., 2017; van Kampenhout et al., 2017).

Several firn-densification equations have been derived and calibrated for the GrIS specifically. We favoured the use of Eq. (1) as our base case because (i) most of these calibrated schemes were developed for dry firn densification, whereas the CROCUS formulation accounts for the presence of liquid water explicitly; (ii) applying a percolation scheme in a stress-based densification model rather than in an accumulation-rate-based model ensures that the redistribution of mass associated with percolation will affect the densification appropriately; and (iii) the CROCUS densification scheme is currently used by the regional climate model (RCM) MAR and by the Community Earth System Model (CESM) to quantify firn densification on the GrIS (Fettweis et al., 2017; van Kampenhout et al., 2017).

2.2 Climatic forcing

To force the model at its upper boundary we use three-hourly skin temperature, melt, snowfall, rain and sublimation fields simulated by the latest version of the RACMO2 regional climate model (RACMO2.3p2, Noël et al., 2018). This model has a 5.5 km horizontal resolution grid and has been explicitly adapted for use over the polar ice sheets. Discrepancies between the climatic forcing and the real climatic history can bias the firn models' results. For the areas of our study sites (Sect. 2.4), we provide detailed statistical comparisons between RACMO2.3p2 output and available observations (see Sect. S1 in the Supplement). We further refer to Noël et al. (2018) for more discussion about the performance of RACMO2.3p2 on the GrIS scale and related uncertainties. Additionally, Ligtenberg et al. (2018) have demonstrated the impact of recent developments in RACMO2.3p2 on firn modelling, mostly yielding improvements in modelled densification.

If the solid input rate (snowfall – sublimation) is negative over a time step, the CFM treats it as a corresponding mass loss in the surface layer; liquid water is evaporated before solid mass gets sublimated. The temperature of a newly accumulated snow layer is defined as the skin temperature at that time step. Deep firn temperatures in the model are thus mostly determined by the mean surface temperature applied during the spin-up process (Sect. 2.5) together with latent heat release through refreezing. We use a Neumann boundary condition for the temperature at the bottom of the domain and use a 250 m deep column to account for the large thermal mass of the ice sheet during the transient run.

In addition to latent heat release due to refreezing, the CFM accounts for heat conduction through the different layers to determine the temperature profile. In accordance with previous firn modelling studies (Kuipers Munneke et al., 2015a; Steger et al., 2017a), we make the firn conductivity, k_s, a function of density following Anderson (1976):

$$\begin{array}{}\text{(4)}& k_{\mathrm{s}}=\mathrm{0.021}+\mathrm{2.5}{\left({\displaystyle \frac{\mathit{\rho }}{\mathrm{1000}}}\right)}^{\mathrm{2}}.\end{array}$$

Another boundary condition is the density of every fresh snow layer deposited at the surface. To reduce the sources of possible uncertainties, we simply use a constant and site-specific surface density according to the surface value of the drilled firn cores instead of a parameterised formulation (see Sect. 2.4).

2.3 Grain size

The temporal evolution of grain size in firn is poorly understood and observational constraints are scarce. However, the grain size is a key variable for the RE, and the flow schemes used in this study thus require an initial grain size and a grain growth rate. For the former, we use the empirical formulation of Linow et al. (2012) derived from observations of snow samples from Antarctica and Greenland:

$$\begin{array}{}\text{(5)}& r_{\mathrm{0}}=(b_{\mathrm{0}}+b_{\mathrm{1}}\left(T_{\mathrm{av}}-T_{\mathrm{0}}\right)+b_{\mathrm{2}}\dot{b}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}}),\end{array}$$

where ρ_i is the ice density (917 kg m⁻³); ρ_w is that of liquid water (1000 kg m⁻³); $\dot{b}$ is the mean annual accumulation rate (m w.e. yr⁻¹); T_av the mean annual surface temperature (K); and b₀, b₁ and b₂ are calibration parameters taking the values 0.781 m, 0.0085 m K⁻¹ and −0.279 yr (m w.e.)⁻¹ respectively.

For grain growth rate, the relationship proposed by Katsushima et al. (2009) is applied:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}r}{\mathrm{d}t}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{8}{r}^{\mathrm{2}}{\mathrm{10}}^{\mathrm{9}}}}\min \left[{\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}\left(\mathrm{1.28}{\mathrm{10}}^{-\mathrm{8}}+\mathrm{4.22}{\mathrm{10}}^{-\mathrm{10}}{\left({\mathit{\theta }}_{\mathrm{weight},\mathit{\%}}\right)}^{\mathrm{3}}\right),\right.\\ & \left.{\mathrm{6.9410}}^{-\mathrm{8}}\right],\end{array}\end{array}$$

where r is the grain radius (m) and θ_weight,% is the mass liquid water content expressed in percent and is thus related to θ (see Eq. 3):

$$\begin{array}{}\text{(7)}& {\mathit{\theta }}_{\mathrm{weight},\mathit{\%}}=\mathrm{100}{\displaystyle \frac{\mathit{\theta }{\mathit{\rho }}_{\mathrm{w}}({\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho })}{{\mathit{\rho }}_{\mathrm{i}}\mathit{\rho }}}.\end{array}$$

Equation (6) combines a wet snow metamorphism formula and a higher limit of growth rate of ice particles, both derived from laboratory measurements.

To study the sensitivity of the model to the grain-size implementation, we also use an alternative option based on the approach for West Antarctic firn of Arthern et al. (2010); the grain radius in newly deposited layers (r₀) has the constant value of 0.1 mm and the grain growth rate is formulated as

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}r}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}r}}{k}_{\mathrm{g}}\exp \left({\displaystyle \frac{-E_{\mathrm{g}}}{\mathrm{RT}}}\right),\end{array}$$

where E_g is the activation energy for grain growth (42.4 kJ mol⁻¹), R is the gas constant (8.314 J mol⁻¹ K⁻¹) and k_g a parameter that takes the value $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹. Note that Eq. (8) does not take the impact of liquid water presence on firn metamorphism into account.

2.4 Study sites

We perform simulations at four study sites in the percolation zone of the GrIS where the availability of well-documented firn cores allows for model–observation comparisons: NASA-SE, DYE-2, KAN-U and FA13 (perennial firn aquifer) (Fig. 1). NASA-SE (66.48^∘ N, 42.50^∘ W; 2372 m a.s.l.) is located in the upper part of the percolation zone with a mean annual temperature of −20^∘ and relatively low melt rates (50 mm yr⁻¹). DYE-2 (66.48^∘ N, 46.28^∘ W; 2126 m a.s.l.) is a slightly warmer site ($T_{\mathrm{av}}=-\mathrm{18}$^∘), and melt is about 3 times greater than at NASA-SE (150 mm yr⁻¹). KAN-U (67.00^∘ N, 47.03^∘ W; 1838 m a.s.l.) is near the equilibrium line altitude and has warmer temperatures ($T_{\mathrm{av}}=-\mathrm{8}$^∘) and significant melting (280 mm yr⁻¹). FA13 (66.18^∘ N, 39.04^∘ W; 1563 m a.s.l.; $T_{\mathrm{av}}=-\mathrm{13}$^∘) is a location known to contain a firn aquifer (Forster et al., 2014). The persistence of deep saturated layers year-round is due to the coupling of high melt rates (587 mm w.e. yr⁻¹) with high accumulation rates (1002 mm w.e. yr⁻¹) (Kuipers Munneke et al., 2014). Multiple firn core density data exist, a large part of which are available in the SUMup dataset (Montgomery et al., 2018). The selection of these four particular sites is motivated by their variety in climatic and glaciological conditions: a cold site with low melt rates, a cold site with high melt rates, a site close to the equilibrium line with substantial refreezing and a site with the presence of a firn aquifer. We perform transient firn-model simulations for each site until the date that a core was drilled. The cores at NASA-SE and DYE-2 were drilled in spring of the years 2016 and 2017 respectively, as part of the FirnCover project. The cores at KAN-U (Machguth et al., 2016) and at FA13 (Koenig et al., 2014) were drilled in spring 2013. The firn-temperature measurements were given in the sources mentioned for the density data, except at KAN-U for which it comes from the collocated automatic weather station overseen by the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) (van As et al., 2016b). The fixed surface densities for DYE-2, NASA-SE, KAN-U and FA13 are 325, 240, 325 and 365 kg m⁻³ respectively and were taken in accordance with the surface density of the drilled cores.

2.5 Spin-up and domain definition

In order to simulate the evolution of the firn layer in time, we start the transient simulation from an initial state in equilibrium with a reference climate. In accordance with previous GrIS firn studies (Kuipers Munneke et al., 2015b; Steger et al., 2017a), we take the 1960–1979 climate as reference climate because it predates the onset of the general warming of Greenland and the subsequent increase in surface melt. We iterate over the reference climate until 70 m w.e. of snow has been accumulated, which ensures the entire firn column is refreshed. The number of iterations over the reference climate is thus site-specific. This spin-up process starts from an analytical solution for the density profile (Herron and Langway, 1980) with temperatures corrected to account for latent heat release by refreezing (Reeh, 2008). During the spin-up process we use the simple bucket approach, and the more advanced flow schemes, detailed in the next section, are turned on only at the end of the spin-up for the transient simulation. This is because using the advanced schemes over long periods is computationally expensive. The domain on which the flow calculations are applied is a subset of the entire CFM domain; this sub-domain is defined each time the flow routine is called in the transient run. The bottom of the sub-domain is defined as the depth below which all layers have density higher than the pore close-off value (830 kg m⁻³), because infiltration of liquid water becomes negligible at this point. The thickness of the layers deposited in every three-hourly time step determines the vertical resolution, and we apply a merging process only to individual layers less than 2 cm thick (see Sect. S2.8).

3 Liquid water schemes

The water flow schemes are added to the dry-densification model detailed in Sect. 2.1 and are thus also effectively one-dimensional, representing no lateral exchange of heat and mass although lateral runoff is used as a mass sink. In this section, we present the three different flow schemes that we implement in the CFM: (1) the bucket method (BK), (2) a single-domain Richards equation scheme (R1M) and (3) a dual-permeability Richards equation scheme (DPM). Because of its robustness and ease of implementation, BK is the current state of the art in firn-densification models that are interactively coupled to regional climate models. R1M is used in several stand-alone snow models to describe water flow (Hirashima et al., 2010; Wever et al., 2014; D'Amboise et al., 2017), and DPM is entirely based on the scheme implemented in the snow model SNOWPACK (Wever et al., 2016), where dual-permeability means that separate domains for matrix flow and preferential flow coexist with liquid water exchanged between these domains.

3.1 Bucket model

The bucket percolation scheme is commonly used to account for the vertical transport of meltwater in firn models, though the precise form of its implementation is variable. Each layer in the model can refreeze meltwater according to its “cold content”, i.e. the energy required to raise the temperature of the layer to the melting point. Starting from the surface, the meltwater may percolate through successive layers, thus allowing for refreezing at depth. Meltwater is progressively depleted due to refreezing and retention according to each layers' water-holding capacity, which is the part of the water that is stored in some of the available pore space and not subject to vertical transfer. The water-holding capacity acts as an approximation of the effect of capillary forces on water retention. Percolation proceeds until all the meltwater is stored (refrozen or retained) or until it reaches a layer with a density exceeding the impermeability threshold (780–830 kg m⁻³), at which point all the water in excess is instantly treated as lateral runoff. The BK thus requires two parameters: the water-holding capacity and the impermeability threshold. We test two possibilities for the former and three for the latter. The water-holding capacity can be prescribed by the calculations of Coléou and Lesaffre (1998) for the mass proportion of water in a firn layer, W_w:

$$\begin{array}{}\text{(9)}& W_{\mathrm{w}}=\mathrm{0.057}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho }}{\mathit{\rho }}}.\end{array}$$

This mass proportion is then converted to the water-holding capacity, θ_h:

$$\begin{array}{}\text{(10)}& {\mathit{\theta }}_{\mathrm{h}}={\displaystyle \frac{W_{\mathrm{w}}}{\left(\mathrm{1}-W_{\mathrm{w}}\right)}}{\displaystyle \frac{\mathit{\rho }{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}\left({\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho }\right)}}.\end{array}$$

Using constant values of the water-holding capacity is also common practice (Reijmer et al., 2012; Steger et al., 2017a). Our base case scenario uses a fixed θ_h at 0.02, or 2 % of the pore space available for liquid water retention. This low value assumes effective downward percolation and is meant to account for vertical preferential flow (Reijmer et al., 2012). For that reason, we consider this as a good basis for comparison with the DPM that explicitly accounts for such flow.

We test three values for the impermeability threshold; these were selected in accordance with Gregory et al. (2014), who tested firn permeability of Antarctic samples in a lab and reported that impermeability can occur over density values ranging from 780 to 840 kg m⁻³. We thus take our three test values to be 780, 810 and 830 kg m⁻³, respectively the lower bound and middle of this range and a commonly used value of pore close-off density.

3.2 Richards equation

Vertical movement of water in a variably saturated porous medium can be described by the one-dimensional version of the RE:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left[K\left(\mathit{\theta }\right)\left({\displaystyle \frac{\partial h}{\partial z}}+\mathrm{1}\right)\right]=\mathrm{0},\end{array}$$

where K is the hydraulic conductivity (m s⁻¹), h is the pressure head (m), z is the vertical coordinate (m, taken positive downwards) and θ is as defined in Eq. (3). The +1 term accounts for the effect of gravity. The RE is an equation expressing the mass conservation law and Darcy's law, and it includes the “suction head”, i.e. the suction force exerted at the surface of individual grains.

A water-retention curve describes the relationship between θ and h required by Eq. (11). We use the van Genuchten (1980) model, which is typically applied in studies of liquid water flow through snow (Jordan, 1995; Hirashima et al., 2014; Wever et al., 2014; D'Amboise et al., 2017):

$$\begin{array}{}\text{(12)}& \mathit{\theta }={\mathit{\theta }}_{\mathrm{r}}+({\mathit{\theta }}_{\mathrm{sat}}-{\mathit{\theta }}_{\mathrm{r}}){\displaystyle \frac{(\mathrm{1}+{\left(\mathit{\alpha }{\left|h\right|}^n\right)}^{-m}}{\mathrm{Sc}}},\end{array}$$

where θ_r is the residual water content (m³ m⁻³) and θ_sat is the volumetric liquid water content at saturation (m³ m⁻³). Sc is a correction coefficient following Wever et al. (2014). The parameters α, n and m are tuning coefficients, with α being related to the maximum pore size and n and m being related to the pore size distribution. These three parameters, referred to as the van Genuchten parameters, are specific to the modelled porous medium and for snow; a common approach is to use the parameterisation developed by Yamaguchi et al. (2012) in a laboratory study:

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathit{\alpha }=\mathrm{4.4}{\mathrm{10}}^{\mathrm{6}}{\left({\displaystyle \frac{\mathit{\rho }}{\mathrm{2}r}}\right)}^{-\mathrm{0.98}},\text{(14)}& {\displaystyle }n=\mathrm{1}+\mathrm{2.7}{\mathrm{10}}^{-\mathrm{3}}{\left({\displaystyle \frac{\mathit{\rho }}{\mathrm{2}r}}\right)}^{\mathrm{0.61}},\text{(15)}& {\displaystyle }m=\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{n}}.\end{array}$$

Yamaguchi et al. (2012) measured the water-retention curve for a range of grain radii (0.025 to 2.9 mm) and densities (361 to 636 kg m⁻³) in different snow samples by using a gravity drainage column method.

The porosity is the part of the volume not occupied by the solid matrix and, in the case of firn, is defined as

$$\begin{array}{}\text{(16)}& P=\mathrm{1}-{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{i}}}}.\end{array}$$

The volumetric liquid water content at saturation is proportional to the porosity (Wever et al., 2014):

$$\begin{array}{}\text{(17)}& {\mathit{\theta }}_{\mathrm{sat}}=P{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}}.\end{array}$$

Note that water is not assumed to fill the entire pore space in saturated conditions and the correction factor $\frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}$ included in Eq. (17) accounts for the required space to allow the liquid water to freeze. It is reasonable to use this correction factor since Yamaguchi et al. (2010) found that trapped air still occupies 10 % of the porosity in saturated snow.

The parameter θ_sat thus represents the pore space available for liquid water, and from there we can define the effective saturation as

$$\begin{array}{}\text{(18)}& \mathrm{Se}={\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{sat}}-{\mathit{\theta }}_{\mathrm{r}}}},\end{array}$$

and Se must be bounded between 0 and 1. In completely dry layers, a zero effective saturation would lead to infinite values in the head pressure calculation, and, thus, we use a numerical adjustment to avoid this happening (see Sect. 2.3). The residual water content θ_r is defined as the amount of liquid water that cannot be removed by gravity as it is held by capillary tension at the surface of the solid grains. Following Yamaguchi et al. (2010), a constant value of θ_r=0.02 can be taken, but, in case of refreezing, θ can approach zero and θ_r must be adjusted accordingly. We take θ_r following a piecewise function:

$$\begin{array}{}\text{(19)}& {\mathit{\theta }}_{\mathrm{r}}=\min \left[\mathrm{0.02},\mathrm{0.9}\mathit{\theta }\right].\end{array}$$

The numerical requirement of an effective saturation value strictly greater than zero causes the persistence of very low flow rates, even for liquid water contents close to the residual water content. Over long time periods, layers cannot hold any residual water content and eventually dry out under the effect of gravity. By taking the coefficient 0.9 in Eq. (19) instead of 0.75 used in snow models (Wever et al., 2014; D'Amboise et al., 2017), we partially reduce this effect because this lowers the effective saturation value Se for any value of volumetric water content θ approaching zero.

The hydraulic conductivity (K(θ)) is the ability of the fluid to flow through the porous medium under a certain hydraulic gradient dependent on pressure head and gravity. Thus, K(θ) depends on the effective saturation and on the properties of both the porous medium and the fluid; fluid flow is enhanced in highly saturated layers. The hydraulic conductivity is described by the van Genuchten–Mualem model (Mualem, 1976; van Genuchten, 1980):

$$\begin{array}{}\text{(20)}& K\left(\mathit{\theta }\right)=\mathrm{Ksat}{\mathrm{Se}}^{\mathrm{1}/\mathrm{2}}{\left[\mathrm{1}-{\left(\mathrm{1}-{\mathrm{Se}}^{\frac{\mathrm{1}}{m}}\right)}^m\right]}^{\mathrm{2}},\end{array}$$

where Ksat is the hydraulic conductivity in saturated conditions (Se=1). For the case of water flow through snow, it has been inferred using three-dimensional images of the microstructure by Calonne et al. (2012) as

$$\begin{array}{}\text{(21)}& \mathrm{Ksat}=\mathrm{3.0}{r}^{\mathrm{2}}\exp (-\mathrm{0.013}\mathit{\rho })\left({\displaystyle \frac{g{\mathit{\rho }}_{\mathrm{w}}}{\mathit{\mu }}}\right),\end{array}$$

where g is the gravitational acceleration (9.8 m s⁻²) and μ is 0.001792 kg m⁻¹ s⁻¹, the dynamic viscosity of liquid water at 273.15 K. Equation (21) shows that simulated water flow is faster in layers with coarser grains and lower densities. These conditions correspond to cases where the connectivity between the pore spaces is high. With respect to the hydraulic conductivity parameterisation, we additionally modify the permeability of ice layers. The hydraulic conductivity of any layer with density exceeding the impermeability threshold is set to zero, rendering it impermeable to incoming flow and leading to the ponding of water on top of that layer and subsequent enhancement of runoff rates (Sect. 3.4.2). This RE implementation completely describes R1M and provides the basis of DPM, further detailed in the next section.

Details of the numerical implementations that are required to maintain stability and to improve computational efficiency for the RE calculations are discussed in the Supplement.

3.3 Dual-permeability model

Physical models of preferential flow in snow are still scarce (Hirashima et al., 2014; Wever et al., 2016). In this section, we explain how the SNOWPACK dual-permeability model (Wever et al., 2016) is implemented in the CFM. The firn column is separated into two domains and water flow in both is governed by the RE (Sect. 3.2). We define F as the pore space allocated to the preferential flow domain and accordingly 1−F as the pore space for the matrix flow domain. Wever et al. (2016) used a grain-size dependence for F, but their regression was performed on only four data points measured in idealised snow laboratory conditions (Katsushima et al., 2013). The experimental grain sizes ranged from 0.1 to 0.8 mm and the water input from 480 to 550 mm per day, which is not representative of firn conditions in Greenland (Figs. 1 and 2). Moreover, due to the typical grain-size ranges in firn (Gow et al., 2004; Lyapustin et al., 2009), the model would regularly be forced to use for F the minimal value for numerical stability implemented in SNOWPACK. To deal with this uncertain parameter but still retain fidelity with respect to the SNOWPACK implementation, we favour the use of a constant value based on observations in natural snow. Marsh and Woo (1984) and Williams et al. (2010) reported that rapid flow paths occupy respectively 22 % and 5 % to 30 % of the area, and we thus fix the value F=0.2. The extension of the preferential flow area within the snowpack is very likely to be a function of grain size and meltwater influx, and these dependencies are still uncertain (Avanzi et al., 2016). The value of F thus determines the value of the saturated liquid water content θ_sat in both domains, and, instead of Eq. (17), we write

$$\begin{array}{}\text{(22)}& \left\{\begin{array}{l}{\mathit{\theta }}_{\mathrm{sat},\mathrm{m}}=(\mathrm{1}-F)P{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}}\\ {\mathit{\theta }}_{\mathrm{sat},\mathrm{p}}=FP{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}}\end{array}\right.,\end{array}$$

where from hereon, the subscripts m and p stand for matrix and preferential flow domain respectively. Equation (22) shows that the volumetric water content in the preferential flow domain is smaller than that in the matrix flow domain. All the input of meltwater is added to the matrix flow domain. For the regulation of the exchange of water between domains, we also closely follow the transfer processes of SNOWPACK (Wever et al., 2016) which are executed at the same 15 min time step. We briefly summarise the transfer processes below.

Figure 2 Annual surface mass fluxes from RACMO2.3p2 at the study sites (1980–drilling date).

Water from the matrix flow domain can enter the preferential flow domain of the layer below if the pressure head in the layer reaches the water entry suction, h_we, of the underlying layer. The parameter can be expressed as (Katsushima et al., 2013; Hirashima et al., 2014; Wever et al., 2016)

$$\begin{array}{}\text{(23)}& h_{\mathrm{we}}=\mathrm{0.0437}{\left(\mathrm{2}r\right)}^{-\mathrm{1}}+\mathrm{0.01074}.\end{array}$$

The amount of water transferred into the preferential flow domain equals the amount of water in excess of h_we. If after the transfer Se in the matrix flow domain still exceeds Se in the preferential flow domain of the underlying layer, their respective Se's are equalised by transferring the appropriate amount of water from the overlying matrix flow domain to the underlying preferential flow domain. In addition, in every individual firn layer where Se in the matrix flow domain exceeds Se in the preferential flow domain, matrix and preferential Se's are equalised by transferring water from the matrix flow domain to the preferential flow domain. This serves to avoid the presence of horizontal pressure gradients in wet snow.

Water can flow from the preferential flow domain to the matrix domain by two processes. The first process is when the saturation in the preferential flow domain exceeds a threshold value Θ. Wever et al. (2016) determined Θ by tuning its value to best match observations. When this threshold is reached, the amount of water corresponding to the cold content of the layer flows back into the matrix domain. If there is still water in excess of the threshold in the preferential flow domain, saturation in both domains is set equal to one another. The second process simulates the heat flow from the preferential flow domain (at the melting point) to the colder surrounding matrix domain. Instead of transferring sensible heat, this process allows liquid water and its inherent latent heat to be exchanged to account for a theoretical heat flow, Q, and thus approximating Fourier's law:

$$\begin{array}{}\text{(24)}& Q=k_{\mathrm{s}}{\displaystyle \frac{(T-T_{\mathrm{0}})}{\left(\sqrt{\frac{\mathrm{1}+F}{\mathrm{2}\mathit{\pi }}}-\sqrt{\frac{F}{\mathit{\pi }}}\right)}}.\end{array}$$

This formulation assumes a linear horizontal temperature gradient in the matrix and a circular shape of the preferential flow path's perimeter. From Eq. (24), the corresponding water transfer is calculated as

$$\begin{array}{}\text{(25)}& \mathrm{\Delta }{\mathit{\theta }}_{\mathrm{p}\to \mathrm{m}}={\displaystyle \frac{\mathrm{2}N\sqrt{\mathit{\pi }F}Q\mathrm{\Delta }{t}_{\mathrm{15}}}{L_{\mathrm{f}}{\mathit{\rho }}_{\mathrm{w}}}},\end{array}$$

where Δt₁₅ is the 15 min time step (s), L_f is the specific latent heat of fusion (335 500 J kg⁻¹) and N is a tuning parameter representing the number of preferential flow paths per square metre (m⁻²). In their study, Wever et al. (2016) arrived at a best parameter set for Θ and N of 0.1 and 0 m⁻² based on comparisons of ice-layer occurrence and runoff amounts with observations in alpine snowpacks. Note that the use of a null value for N is implausible in our case of firn-column simulations. Indeed, this would imply that liquid water would persist and flow deeper in the preferential flow domain in saturation conditions below the Θ value until the bottom of a subfreezing firn column, which can be up to 70 m thick in some areas of the GrIS. Therefore, we use the smallest non-zero value of N tested by Wever et al. (2016), and the parameters Θ and N are fixed to 0.1 and 0.2 m⁻² respectively.

The hydraulic conductivity of ice layers is not artificially set to zero in the preferential flow domain as it is in the matrix flow domain. Preferential flow thus provides a way for water to flow through an ice layer, reproducing observations that ice layers are not totally impermeable barriers and can lead to localised piping events (Marsh and Woo, 1984; Pfeffer and Humphrey, 1998; Williams et al., 2010; Sommers et al., 2017). An exception for this is the bottom of the domain: as preferential flow is stopped at the last layer, it does not percolate through the solid ice.

3.4 Additional processes in the single- and dual-domain schemes

3.4.1 Refreezing process

In R1M and DPM a cold content is calculated for every firn layer, similarly to BK (Sect. 3.1), and refreezing in both flow schemes is executed at the 15 min time step (equivalent to the time step of the transfer processes of DPM, Sect. 3.3).

When refreezing occurs, every layer freezes the maximum of its liquid water content that its cold content allows. For numerical reasons, refreezing cannot dry out a layer completely; instead, a very low value of liquid water remains in every layer (see Sect. S2.3). The refrozen water densifies the firn layer and modifies its hydraulic properties. The remaining liquid water is still subject to flow and infiltrates deeper into the firn column.

In DPM, refreezing is restricted to the matrix flow domain (see Sect. S2.7). In the preferential flow domain, liquid water can percolate through cold layers, as has been observed in field studies on the GrIS (e.g. Pfeffer and Humphrey, 1996; Humphrey et al., 2012). For this liquid water to refreeze, it first has to be transferred back to the matrix flow domain. Preferential flow thus provides a way for liquid water to bypass cold firn layers and subsequently to infiltrate deeper layers.

3.4.2 Aquifer development and lateral runoff

In R1M and DPM, lateral runoff in the firn column is simulated using the parameterisation of Zuo and Oerlemans (1996):

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{\mathrm{dRu}}{\mathrm{d}t}}={\displaystyle \frac{L_{\mathrm{excess}}}{{\mathit{\tau }}_{\mathrm{Ru}}}},\text{(27)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{Ru}}=c_{\mathrm{1}}+c_{\mathrm{2}}\exp (-c_{\mathrm{3}}S),\end{array}$$

where Ru is the amount of meltwater that runs off (m), L_excess is the excess of liquid water amount with respect to the residual water content (m) and τ_Ru is a characteristic runoff time (s). The constants c₁, c₂ and c₃ are parameters derived by comparison with observations by Zuo and Oerlemans (1996) for the GrIS, and S is the surface slope. The meltwater input is immediately treated as lateral runoff if the surface layer is an impermeable ice layer or if it is saturated.

Equation (26) leads to the complete drainage of a layer with a zero slope in only 26 d, precluding the formation and persistence of perennial firn aquifers. Therefore, we do not apply Eq. (26) in the layers at the bottom of the firn column. All the inflow of water reaching this section is added to the aquifer, allowing the model to progressively fill the pore space of the bottom layers in the firn column with meltwater.

Firn aquifers are known to be affected by drainage mechanisms not represented in the model, for example via crevasses (Poinar et al., 2017), and possibly hydrofracture and rapid drainage events (Koenig et al., 2014). Miller et al. (2018) found discharge rates within the firn aquifer to be $\mathrm{4.3}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹ by borehole dilution tests in the field. We tested this approach in our model by applying this value as a constant discharge rate for aquifers formed in our simulations. We found however that, using this approach, an aquifer was not sustained; suggesting that such discharge rates must be dependent on the total amount of water within the aquifer and are likely temporally variable. To account for drainage processes and yet allow the formation of an aquifer, we therefore limited the amount of water stored in the firn aquifer to 1.65 m w.e., i.e. the water level measured in the field by Koenig et al. (2014). In firn aquifers forming at the bottom of the firn column, the saturation in both domains is equalised and the model does not perform flow calculation in this lowest part of the domain (see Sect. S2.5).

3.5 Investigating model sensitivity

In Sects. 2 and 3, we highlight several factors influencing BK, R1M and DPM. For each of the schemes, we analyse results generated using three possible impermeability thresholds: 780 kg m⁻³ (ip780), 810 kg m⁻³ (ip810) and 830 kg m⁻³ (ip830). This provides a way to compare the sensitivity of the simple BK and of the physically based schemes (R1M and DPM) to a common parameter. For BK, we try two different formulations of the water-holding capacity: constant at 0.02 (wh02) and according to the parameterisation of Coléou and Lesaffre (1998), Eq. (9) (whCL). For R1M and for DPM, we test two different grain-size implementations: the Linow et al. (2012) surface grain-size calculation, Eq. (5), coupled to the Katsushima et al. (2009) grain growth rate, Eq. (6) (grLK); and the grain-size implementation of Arthern et al. (2010), Eq. (8) (grA). It is important to examine model sensitivity to the grain-size variable as almost all the hydraulic parameters of the RE depend on it. The different sensitivity tests are summarised in Table 1.

Table 1 Summary of the sensitivity tests.

4 Results

In this section, we describe and discuss the model performance at each of the four sites tested (DYE-2, NASA-SE, KAN-U and FA13). We begin by comparing BK, R1M and DPM in a base case parameterisation: BK wh02 ip810, R1M grLK ip810 and DPM grLK ip810 respectively. Then, we perform various tests to investigate the sensitivity of the flow schemes to variations in their parameter values. We refer to ice layers as layers with a density value exceeding the impermeability threshold in the model and to liquid water input as the total of meltwater and rain influx. The DPM approach features two tuning parameters, N and Θ. Model results and depth–density profiles were found to be weakly sensitive to the value of N and Θ, and so we omit consideration of these from the remainder of our study. Results of simulations and observations are inter-compared based on the firn air content (FAC; the depth integrated porosity in a firn column) over the top 15 m of firn and the temperature at 10 m depth. Comparing the modelled FAC and 10 m depth temperature values with observed data depicts the ability of the tested models to reproduce the bulk condition of the upper firn column. We also qualitatively assess the degree to which the models to form a realistic ice-layer distribution and depth–density profile. One would not expect simulated values of either to match observations precisely given the high spatial variability of firn structure (Marchenko et al., 2017), but it is indicative of the models' performance in reproducing heterogeneity in firn density.

4.1 DYE-2

DYE-2 has a typical liquid water input between 0.1 and 0.3 m w.e. yr⁻¹ (Fig. 2), which is moderate in the context of our study sites. The extreme melt year of 2012 (Nghiem et al., 2012) is an exception, with an estimated input of more than 0.7 m w.e. Using BK, almost all of this meltwater refreezes locally and runoff is close to zero (Table 2) until the 2012 summer when ice layers (ρ≥810 kg m⁻³) start forming in the top 2 m (Fig. 3a). Runoff increases in the subsequent years because meltwater reaches these ice layers. In R1M and DPM, small amounts of runoff occur between 1980 and 2011 due to the lateral runoff implementation, Eq. (26). Beginning in summer 2004, some ice layers start to form in R1M (Fig. 3b) due to the refreezing of water held close to the surface by capillary forces. Over the 2012 summer, surface layers are progressively melted, bringing ice layers closer to the surface. The ponding and refreezing of water on the top ice layer allows it to thicken. This then acts as an impermeable barrier to vertical percolation from 2012 onwards, resulting in a more than 6-fold increase in runoff (Table 2). In contrast, runoff remains low in DPM, in which several ice layers form in the upper firn as early as summer 1996 (Fig. 3c). These ice layers generally form deeper than 2 m due to more effective water transfer from the near-surface to lower layers; preferential flow provides a path for ponding meltwater in the matrix flow domain to bypass ice layers and continue to percolate vertically, thus maintaining low runoff amounts. Preferential flow brings part of the 2012 meltwater to depths greater than 12 m. For each flow scheme, the modelled FAC underestimates the observed value by 4 %–16 %. This can partly be attributed to the tendency of the CROCUS scheme to slightly overestimate densification rates in the upper part of polar firn (Gascon et al., 2014). FAC is underestimated more strongly in DPM (16 %) than in BK and R1M (4 %) because in DPM the deeper firn is not isolated from surface meltwater percolation (Table 2).

Table 2 Model outputs at DYE-2 site. A slash indicates no data.

Figure 3 Modelled firn density at DYE-2. (a) BK wh02 ip810, (b) R1M grLK ip810, (c) DPM grLK ip810, (d) BK whCL ip810, (e) R1M grA ip810, and (f) DPM grA ip810; black indicates solid ice.

Modelled density profiles using each flow scheme are compared with observations (Fig. 4a). Mean density is reproduced reasonably well using each of the three flow schemes, but no configuration is able to qualitatively reproduce the strong variability in density observed. For example, numerous high-density layers separated by much-lower-density intervals are clear in the observations. Regardless of the flow scheme, only a few ice layers are formed in the model and these tend to be confined to the upper 6 m, which has been affected by the higher melting rates of the recent years. In older firn deposited under lower-melt conditions, the number of density peaks and their amplitude is underestimated even more strongly. Several ice layers are observed in the 10–20 m depth range; only DPM simulates the presence of ice layers here.

Figure 4 Measured and modelled depth–density profiles at DYE-2 on 11 May 2017. Thick vertical lines show ice layers. The modelled densities are averaged at the vertical resolution of the drilled core. CR: CROCUS, HL: Herron and Langway, KM: Kuipers Munneke.

The three flow schemes lead to significantly different firn thermal conditions. The temperatures at 10 m depth of BK and R1M agree well with observations (+0.2 and −0.4 K). In contrast, 10 m temperature is strongly overestimated in DPM (+2.7 K) because it allows percolation at depth, subsequent refreezing and latent heat release. The summer 2012 percolation raises the 10 m depth temperature to within a few degrees of melting using DPM. Since the DPM method seems to exaggerate deep percolation, we tested a lower impermeability threshold (DPM grLK ip780) which should favour the formation of shallow ice layers, the ponding of water in the matrix flow domain, more lateral runoff and colder temperatures at depth. The ice layers do form slightly earlier in the melt seasons but are not noticeably shallower than in DPM ip810. The partitioning between runoff and refreezing is barely affected and the 10 m temperature bias remains (Table 2).

The BK method gives a density profile closer to R1M than to DPM. In order to mimic the behaviour of DPM we increase the impermeability threshold in BK (BK wh02 ip830) to make it more effective in transporting water vertically; however, model results are only weakly affected by this change (Table 2). We also modify the water-holding capacity in BK according to the parameterisation of Coléou and Lesaffre (1998) (BK whCL ip810) which allows more water to be retained in the low-density layers close to the surface. Ice layers appear earlier in the simulation and at shallower depths (Fig. 3d). This increases the amount of runoff in BK whCL ip810 with respect to BK wh02 ip810 (+4 % of the water input over the entirety of the transient model run); however, in the surface layers, where high amounts of water are retained, refreezing dominates. As a result, much less water percolates to the deeper firn and there is less refreezing and latent heat release. All of this leads to a significantly higher FAC (+4 %) and colder 10 m temperature (−1.5 K) relative to BK wh02 ip810.

For models based on the RE (R1M and DPM), we test sensitivity to grain size by implementing a parameterisation for grain growth based on Arthern et al. (2010) (denoted grA). Using this parameterisation, grain sizes tend to be smaller, and so more water tends to be retained and refrozen close to the surface due to stronger capillary forces. Compared to the R1M grLK ip810 experiment, the R1M grA ip810 causes formation of ice layers earlier in the simulation (beginning in 1996) and shallower in the firn column (Fig. 3e), favouring water ponding and subsequent runoff (+7 % of the water input over the entirety of the transient model run). Stronger capillarity also means that saturation is higher for percolation to occur, which in turn increases the simulated runoff since more water is in excess of the residual water content. The enhanced runoff and shallower percolation lead to a higher FAC (+4 %) and a colder 10 m temperature (−0.6 K). In DPM, the flow and refreezing patterns are also altered by the grain-size formulation: DPM grA ip810 produces ice layers much earlier (beginning in summer 1981), at shallower depths and in larger numbers (Fig. 3f). Runoff is however only slightly increased (+2 %). The FAC remains similar to DPM grLK ip810, but the 10 m temperature is 0.3 K lower and the warm bias is thus reduced (an 11 % decrease) (Table 2).

Finally, we investigate differences in the depth–density profiles simulated at DYE-2 attributed to different firn-densification formulations in contrast to those observed due to the use of different flow schemes. We first choose to apply the DPM grLK ip810 flow scheme with the additional firn-densification formulations of Herron and Langway (1980) (HL) and of Kuipers Munneke et al. (2015) (KM), both calibrated for GrIS firn. The HL and KM models are forced with the same three-hourly climatic forcing as the CROCUS (base case) model. The FAC (−5 %), 10 m temperature (+0.4 K) and mean density profile (Fig. 4b) predicted by the HL densification model agree reasonably well with that predicted by CROCUS, although HL predicts greater density variability due to its stronger dependence on the annual temperature cycle. In contrast, the KM model predicts much higher densification rates and thus greater densities, with several thick ice layers in the 3–8 m depth range, some exceeding a metre thickness. This results in a much lower FAC value compared to the CROCUS model (−24 %), and, in this case, differences between flow schemes are small with respect to the choice of the densification formulation. Since the warm bias of DPM can cause temperature-dependent densification formulations to overestimate densities, we also compare the three densification formulations coupled to R1M grLK ip810 (Fig. 4c). Similar to the results using the DPM flow scheme, the HL profile agrees reasonably well with the CROCUS model (FAC value is −7 %) but predicts that a metre-thick ice layer formed at 5 m depth (Fig. 4c) during the 2012 summer. Discrepancies between CROCUS and KM are only slightly reduced using R1M; for example, the FAC predicted by KM is 20 % less than that predicted by CROCUS. This can be attributed to greater densities at depth (>8 m) and to much higher densities in the depth range 3–5 m. The latter corresponds to the layers affected by meltwater refreezing and considerable latent heat release in the 2012 summer.

4.2 NASA-SE

NASA-SE is a site characterised by high accumulation rates, ranging between 0.5 and 0.8 m w.e. yr⁻¹, and low rates of liquid water input, typically between 0.01 and 0.15 m w.e. yr⁻¹ (Fig. 2). Under these conditions, abundant pore space and cold content are available for prompt refreezing of the summer meltwater, so one would expect a smaller sensitivity of the model to the flow scheme applied. In BK, no runoff is produced over the entire simulation (Table 3) since refreezing of small amounts of melt does not lead to the formation of impermeable ice layers. R1M and DPM have very low runoff amounts with a small spike in the summer of 2012 when there was 0.38 m w.e. of liquid water input. No ice layer forms in the top 15 m of the firn column using any of the liquid water schemes, in agreement with the observed core (Fig. 5a). Changing the impermeability threshold results in identical model results since no layer exceeds the lowest possible value in the depth range where water percolates. The three water-transport schemes predict a similar FAC; they all underestimate the observed value by approximately 3 % (Table 3). This is because the mean firn density is well-captured by the model but somewhat overestimated at depths greater than 8 m (Fig. 5a). R1M simulates a single density peak at 8 m depth (Fig. 5a), corresponding to the 2012 summer meltwater percolation, due to capillary forces effectively retaining the relatively high meltwater volume produced in that year close to the surface and exposing it to delayed refreezing once these layers cool below the freezing point. DPM also produces a density peak (albeit a much smaller one) at a similar depth, and more-effective downward percolation results in a uniform increase in density over the next 3 m. Finally, BK also produces a small density peak; however, this is at a greater depth of 9 m since it assumes water flow to be instantaneous in a time step and the major part of the refreezing occurs as water reaches deeper cold layers. Again, none of the percolation schemes capture the observed variability in density. Also, despite the low melt/accumulation ratio, the three percolation schemes overestimate the 10 m temperature by 1.4–2.2 K (Table 3).

Table 3 Model outputs at NASA-SE site. A slash indicates no data.

Figure 5 Measured and modelled depth–density profiles at NASA-SE on 4 May 2016. The modelled densities are averaged at the vertical resolution of the drilled core. CR: CROCUS, HL: Herron and Langway, KM: Kuipers Munneke.

Increasing the water-holding capacity in BK (BK whCL ip810) leads to a minor increase in the FAC (<1 %) and a 0.9 K cooling of the 10 m temperature, because the surface layers have a relatively low density (surface boundary condition of 240 kg m⁻³ at this site) and thus retain high amounts of water with the whCL parameterisation (Table 3). The R1M and the DPM density profiles are weakly sensitive to a change in the grain-size formulation from grLK to grA (Table 3). This is due to the small meltwater amounts with meltwater refreezing only slightly closer to the surface because of the stronger capillarity retention in the grA models. However, we note that simply changing the grain-size formulation in R1M from grLK to grA leads to a 0.4 K colder 10 m temperature and thus decreases the bias with respect to observations by 28 % (Table 3).

We used the R1M grLK ip810 model with the HL and the KM densification formulations in order to prevent the DPM's warm bias from skewing the modelled densification rates. As expected in this relatively dry site, the modelled profiles are much more sensitive to the dry densification than to the percolation scheme (Table 3 and Fig. 5a and b). The maximal difference in FAC between the three densification formulations tested is 20 % compared to less than 1 % between the three flow schemes and their possible parameterisations. In contrast to the DYE-2 simulations, the CROCUS model predicts the fastest densification and thus the lowest FAC. HL and KM predict 20 % and 11 % greater FAC than CROCUS, respectively, and CROCUS is in closest agreement with the observations, underestimating FAC by just 3 %.

Table 4 Model outputs at KAN-U site. A slash indicates no data.

4.3 KAN-U

KAN-U is a high-melt site with an average melt rate over the 1980–2013 period of 0.33 m w.e. yr⁻¹, and, in the last 3 years of our simulation (2010–2013), the RCM calculates annual melt exceeding annual accumulation (Fig. 2). Since surface temperatures are relatively high (annual mean around −8^∘), refreezing of the summer meltwater depletes the cold content over large depth ranges. Beginning in summer 1990 in the BK simulation, some ice layers are present in the depth range 3–8 m (Fig. 6a), allowing part of the meltwater to run off and impeding percolation to greater depths. At the start of 2012, there is a thick ice layer in the upper 4 m and another one forms at the surface during the summer. As a result, refreezing is constrained to the uppermost firn layers and a large part of the water input runs off (Table 4). In R1M, the high water content and the almost-continuous presence of ice layers in the upper 5 m from summer 1986 onwards (Fig. 6b) cause relatively high runoff rates throughout the simulation (28 % of the water input over the entirety of the transient model run). As in the BK simulation, runoff is particularly high in 2012 due to ice layers impeding vertical percolation below 1 m (Table 4). In the DPM simulation, the preferential flow mechanism leads to the formation of multiple ice layers in the depth range 4–10 m from 1987 onwards (Fig. 6c). Runoff rates remain low but there is a notable increase in 2012. This is due to the formation of ice layers close to the surface, which allows ponding of water in the matrix flow domain. The preferential flow domain is unable to accommodate all the ponding water, and part of it is treated as lateral runoff (Eq. 26). While matrix flow typically remains constrained to the upper 5 m (Fig. 7a), the recent (2010 to 2012) high-melt summers cause preferential flow to reach much greater depths (e.g. up to 35 m in the 2012 summer; Fig. 7b). Since preferential flow can transfer water below ice layers, the refreezing process can fill the pore space available at depth, leading to substantial thickening of the ice layers. As a result, the FAC is much smaller in the DPM simulation than in the BK (−39 %) and the R1M (−35 %) simulations, in which runoff limits the amount of meltwater refreezing.

Figure 6 Modelled firn density at KAN-U. (a) BK wh02 ip810, (b) R1M grLK ip810, (c) DPM grLK ip810, (d) BK whCL ip810, (e) R1M grA ip810, and (f) DPM grA ip810; black indicates ice layers.

Figure 7 Volumetric water content at the KAN-U site for DPM grLK ip810 in (a) matrix flow domain and (b) preferential flow domain; note the difference in scales.

The observations reveal a thick, almost-continuous ice slab over the depth range of 1–7 m (Fig. 8a). Below it, the density is more variable but remains generally high, causing a low FAC (Table 4). Both the BK and the R1M simulation significantly overestimate the FAC (+59 % and +50 %). In contrast, the average FAC of the DPM simulation is very close to the observed value (−2 %); however, the DPM density profile shows an almost-continuous ice slab from 3 to 17 m depth (Fig. 8a) and does not reproduce the lower-density intervals observed. This demonstrates an important limitation of the liquid water schemes: since water cannot be retained in layers exceeding the impermeability threshold, these layers can only further densify by the dry-densification mechanism and not by water refreezing. The overestimation of the ice slab thickness in the DPM profile is thus compensated for by the underestimation of its density, which leads to the good agreement with the observed FAC value. BK reproduces the presence of the ice slab at 1 m depth, but it underestimates its thickness and simulated a thick (2 m) low-density region (Fig. 8a). Below the observed ice slab, the agreement with observed average density is reasonable but variability in density is underestimated. Despite also underestimating the thickness of the ice slab, the R1M profile agrees better with the observed density profile: it produces only two thin, low-density layers in the slab and more high-density peaks and ice layers below 7 m, which is in better agreement with the observed density variability.

Figure 8 Measured and modelled depth–density profiles at KAN-U on 28 April 2013. Thick vertical lines show ice layers. The modelled densities are averaged at the vertical resolution of the drilled core. CR: CROCUS, HL: Herron and Langway, KM: Kuipers Munneke.

With respect to the 10 m temperature, the BK method is biased cold but gives results in reasonable agreement with the observations (−1.7 K). This bias is more pronounced in R1M (−2.6 K). In contrast, DPM largely overestimates the 10 m temperature (+4.5 K), as a result of its overestimation of percolation and subsequent refreezing at depth.

Changing the impermeability threshold for DPM (DPM wh02 ip780 and ip830) does not alter the pattern of the modelled depth–density profile, but the corresponding changes to the density of the ice slab have an impact on the FAC (+15 % for ip780 and −9 % for ip830). Other factors further affect the FAC: runoff rates slightly decrease with higher impermeability thresholds (Table 4); and the mass of the ice layers increases the overburden stress on the firn column below, increasing the densification rate. In addition, higher (lower) impermeability thresholds lead to warmer (colder) 10 m temperatures (+1.1 K for ip830 and −1.6 K for ip780), due to enhanced latent heat release. Compared to BK wh02 ip810, decreasing the impermeability threshold (BK wh02 ip780) leads to the formation of ice layers in earlier years and closer to the surface and thus more runoff (+3 % of the water input over the entirety of the transient model run), which in turn increases the FAC (+5 %) and decreases the 10 m temperature (−0.5 K). Increasing the threshold (BK wh02 ip830) has the opposite effect (−8 % for the FAC and +0.8 K for the 10 m temperature compared to BK wh02 ip810). If we instead allow for a greater water-holding capacity (BK whCL ip810), the partitioning between runoff and refreezing remains very similar (Table 4). However, the FAC and the 10 m temperature are changed (+3 % and −1.7 K compared to BK wh02 ip810). The lower temperature is due to latent heat release from refreezing being more concentrated in the surface layers (Fig. 6d). The formation of ice layers earlier in the year and at shallower depths allows part of the underlying firn to remain free of refreezing, which increases the FAC. Furthermore, colder temperatures cause a higher firn viscosity, thus decreasing the densification rates. Since the R1M formulation both overestimates the FAC and underestimates the 10 m temperature, we test an increase in its impermeability threshold (R1M grLK ip830), allowing for deeper percolation. Both the decrease in FAC (−1 %) and increase in 10 m temperature (+0.1 K) compared to R1M grLK ip810 are minor.

With the grA formulation in DPM (DPM grA ip810), water is more efficiently transferred vertically through the preferential flow domain, which causes an increase in the number of ice layers formed during the simulation (Fig. 6f), a slight decrease in FAC (−2 %) and a slight increase in the 10 m temperature (+0.1 K) relative to DPM grLK ip810. The nearly continuous ice slab, which extends to 17 m depth below the final winter accumulation, explains the weak sensitivity of the final FAC and 10 m temperature values of DPM to grain size. In contrast, applying the grA formulation in R1M (R1M grA ip810) leads to a considerable increase in FAC (+11 %) and a decrease in 10 m temperature (−0.7 K) compared to R1M grLK ip810. This is due to higher water content during percolation events and, especially in the most recent years of our simulation, refreezing and ice-layer formation at shallower depths (Fig. 6e). This increases the runoff and isolates the deeper firn from meltwater percolation. As in the cases of DYE-2 and NASA-SE, the change in FAC due to different grain-size formulations in R1M is greater than the change due to switching from BK to R1M (Table 4).

The modelled depth–density profiles also differ according to the densification formulation used (Fig. 8b). We compare the different densification formulations using R1M grLK ip810, thus avoiding the effect of the strong temperature bias of DPM on the densification process. Densification in KM is sensitive to high firn temperatures, and it predicts the highest densities: it produces the highest density values in the ice slab range, the most ice layers below the ice slab and the lowest FAC value (−27 % compared to the CROCUS formulation). HL behaves in a similar way to CROCUS in the upper 5 m, apart from a much-lower-density interval in the 2–2.5 m depth range. In deeper firn, densities simulated using HL tend to lie between those simulated using KM and CROCUS, and its FAC difference with the CROCUS (−9 %) is less than that of KM. The DPM scheme simulates a depth–density profile of an ice slab over a 14 m range, which is in stark contrast with BK and R1M. Apart from this, the choice of the densification formulation has a greater influence on the model than the choice of liquid water scheme and of any of their respective parameterisations presented here, in spite of the high water input at this site.

4.4 FA13

The FA13 site is representative of conditions in the southeast part of the GrIS; it has both high accumulation and high melt rates (mean 1980–2012 rates of 1.09 and 0.64 m w.e. yr⁻¹ respectively, Fig. 2). This favours the insulation of summer percolating meltwater from winter atmospheric temperatures, typically leading to the formation of PFAs (Kuipers Munneke et al., 2014). Here, the initial conditions and the spin-up process cause the deep firn to be close to the melting point at the start of the transient run.

The warm firn, combined with the high water influx, allows liquid water to reach greater depths than at the other sites in all three flow schemes. Additionally, the firn–ice transition depth becomes important in the FA13 simulations. The observed core shows that the 810 kg m⁻³ density is reached and maintained from 24 m depth. The CROCUS densification scheme predicts that this density horizon occurs at 60 m depth. Since CROCUS has been developed for seasonal snow, the densification at high overburden stress is probably not well-captured by the model (Stevens, 2018). Because of this, we base our simulations for FA13 on the HL densification model, which predicts this transition depth to be around 21 m.

Table 5 Model outputs at FA13 site. A slash indicates no data.

The total refreezing rates are similar for the three flow schemes (Table 5). Since the deep firn is close to the melting point, the total refreezing amounts are essentially determined by the cold content provided in winter, and the precise behaviour of the percolation has a minor impact. However, variability of refreezing with depth differs between schemes, which leads to differences in the 15 m FAC values (Table 5) and in the modelled depth–density profiles (Figs. 9a, b, and c and 10a). FAC is consistently underestimated (−23 to −30 %) because firn density is overestimated above 10 m. R1M and DPM overestimate density most strongly with FAC values 9 % and 10 % smaller than BK respectively, and both schemes simulate the presence of a thick ice layer in the upper 10 m of the firn, which is not observed in the core. The BK model produces only a single thin ice layer in the upper 10 m (0.2 m thick at 9 m depth), which is in good agreement with the observations (showing a single thin ice layer at 7.5 m depth). Below 10 m, the modelled densities are generally in better agreement and all the schemes produce several ice layers (Figs. 10a and 9a, b and c).

Figure 9 Firn density at FA13. (a) BK wh02 ip810, (b) R1M grLK ip810, (c) DPM grLK ip810, (d) BK whCL ip810, (e) R1M grA ip810, and (f) DPM grA ip810; black indicates solid ice.

Figure 10 Measured and modelled depth–density profiles at FA13 on 10 April 2013. Thick vertical lines show ice layers. The modelled densities are averaged at the vertical resolution of the drilled core. CR: CROCUS, HL: Herron and Langway, KM: Kuipers Munneke.

In the absence of any shallow ice layer throughout most of the simulation (Fig. 9a, b and c), meltwater is free to percolate through the winter accumulation layers and to deplete their cold content. The flow schemes have different abilities to store liquid water, which leads to small variations in runoff and refreezing rates. In BK, water is retained according to the water-holding capacity (Fig. 11a) and refreezes during subsequent winters. In contrast, DPM allows percolation down to the firn–ice-sheet transition where it ponds to form an aquifer (Fig. 11d and e). This leads to a significant reduction of runoff amounts during the aquifer build-up (−6 % of the water input over the entirety of the transient model run compared to BK), and the water remaining in the firn column is essentially constrained by the maximal amount of water we allow in the aquifer (1.65 m). In theory, the same mechanism could be simulated by R1M, but the percolating water is depleted before it reaches the bottom of the firn column (Fig. 11b). This is due to refreezing, to the lateral runoff parameterisation and to the presence of ice layers in the upper 10 m. No water persists through the winter seasons, which illustrates the model artefact that the effective saturation must be strictly positive for the stability of the RE (Sect. 3.2). Thus, the refreezing rates are slightly lower than in BK since no residual water is stored and later exposed to winter refreezing (Table 5).

Figure 11 Volumetric water content at FA13. (a) BK wh02 ip810, (b) R1M grLK ip810, (c) BK whCL ip810, (d) matrix flow domain of DPM grLK ip810, and (e) preferential flow domain of DPM grA ip810; note the difference in scale for panel (e).

The build-up of the aquifer starts very early (in the summer of 1981) when DPM is turned on in the transient run due to the low refreezing capacity of the deep firn. The depth of this aquifer is constrained by the impermeability threshold applied, which determines where the model places the firn–ice transition. This depth is at 33 m in 1981 and 21 m in 2013, with the decrease being caused by enhanced densification. The aquifer is fed only by preferential flow (Fig. 11e) since matrix flow cannot reach the water table due to runoff, refreezing and the presence of ice layers in the firn column.

From 1994 and onwards the total simulated water content in summer is only regulated by the maximum allowed in the model (1.65 m). Since the water table is at a shallow depth towards the end of the simulation (7.5 m), the propagation from the surface of the cold winter temperatures can refreeze part of the saturated layers. This leads to the formation and progressive thickening of the shallow, thick ice layer. Also, the shallowness of the aquifer causes 23 % of the porosity in the top 15 m to be filled with liquid water and the 10 m temperature to be at the melting point.

The higher impermeability threshold in DPM grLK ip830 increases the depth of the calculated firn–ice transition, producing a deeper aquifer that extends between 12 and 29 m depth at the end of the simulation, similar to the 12–37 m depth range observed by Koenig et al. (2014). Compared to DPM grLK ip810, the increased depth leads to less refreezing in the shallowest layers of the aquifer and thus a higher FAC value (+3 %) and a 10 m temperature below the melting point.

The grain-size formulation following Arthern et al. (2010) (DPM grA ip810) reduces the ability of preferential flow to transport water down to the firn–ice transition but instead favours formation of discrete ice layers in the firn column (Fig. 9f). In this case the aquifer does not start to form until summer 1988, but the final aquifer structure (also between 7.5 and 21 m), the FAC value (−3 % for grA), and the partitioning between refreezing and runoff are similar to those simulated using grLK (Table 5). In R1M, the sensitivity to grain size is noticeable in the firn-structure evolution with differences in ice-layer formation between R1M grLK ip810 and R1M grA ip810 (Fig. 9b and e). The final FAC value (+4 % for grA) and the meltwater partitioning remain similar (Table 5) between R1M grLK and R1M grA, as for the case of DPM grLK and DPM grA. This can be explained by the total refreezing's stronger dependence on the firn thermal structure than on the percolation pattern at this site.

Increasing the water-holding capacity in BK (BK whCL ip810) leads to a significantly lower FAC value (−12 %): more water refreezes in the near-surface layers, which reduces the runoff and enhances densification in the entire underlying firn column. Also, more water remains stored at depth throughout the different winter seasons (Fig. 11c), and some is still present at the end of the simulation (0.09 m) between 16 and 23 m depth. However, this small amount retained by the water-holding capacity is much less than is stored in the saturated layers of the aquifer simulated in DPM.

We compare the three different densification models (CROCUS, HL, KM) using the R1M grLK ip810 flow scheme and these show important differences in the final modelled depth–density profiles (Fig. 10b). CROCUS agrees reasonably well with HL in the top 6 m, but, as mentioned above, it has a strong low-density bias at greater depths. Since CROCUS simulates lower densification rates, its underestimation of the FAC value in the upper 15 m (−21 %) is smaller than in HL (−30 %), but it is clearly not representative of the density conditions below 15 m. KM predicts a firn column below the last winter's accumulation entirely at the ice density. The model thus identifies a firn–ice-sheet transition at shallow depth (∼2 m), which the water can reach before being depleted by the lateral runoff parameterisation and saturated layers can thus build up. This further amplifies the densification since the saturated layers at the transition depth are exposed to refreezing. Hence, in 2012, runoff combined with refreezing exceeds the liquid water input (Table 5). This occurs because some layers wherein water had been stored in previous years reach the 810 kg m⁻³ density, causing the stored water to be considered as runoff by the model. Whereas the FAC values are generally close for the different flow schemes and their parameterisations (maximal difference of 12 %), CROCUS and KM reach values 14 % higher and 34 % lower than HL respectively.

5 Discussion

The three liquid water schemes show consistent behaviour between sites. R1M generally predicts slower downward percolation of water than the other schemes, which leads to more near-surface refreezing, the formation of near-surface ice layers, more lateral runoff, and thus lower densities and lower temperatures in deeper firn. As a result, when compared to observations R1M tends to reach higher FAC values and to underestimate 10 m temperatures. The BK formulation with the Coléou and Lesaffre (1998) parameterisation for the water-holding capacity leads to the same effects, but they are amplified. The underestimation of the 10 m temperature is stronger, suggesting that BK whCL does not allow for deep enough percolation. BK with the lower water-holding capacity (BK wh02) leads to a partitioning of the water input between refreezing and runoff similar to the more complex R1M at the four sites. As a result, the FAC values predicted by BK wh02 and R1M generally agree (maximum difference less than 10 %), as do the temperatures at 10 m depth (maximum difference less than 1 K). The FAC values and 10 m temperatures of R1M at the end of the model runs always lie in the range of the ones obtained with different parameterisations of BK. This suggests that BK can produce results similar to R1M, provided it is parameterised appropriately.

DPM exhibits a different behaviour: it effectively brings water to greater depths, depleting the deep-firn pore space and cold content. Even in the presence of shallow ice layers hindering matrix flow, the preferential flow implementation still ensures efficient vertical water transport, and runoff amounts remain low. This suggests that transfer mechanisms to the preferential flow domain implemented in DPM are more effective in draining ponding water than the lateral runoff parameterisation. Due to large FAC underestimation and 10 m temperature overestimation, the data–model mismatch of DPM with respect to these variables is significantly greater than that of R1M and BK. DPM is better at producing density variability in depth, which is underestimated in all schemes at all sites. Also, in contrast to the two other schemes, DPM can form ice layers even in summers of average melt, and it is able to simulate the persistence of deep saturated firn layers at the FA13 site. In this respect our findings support those of Wever et al. (2016), who highlight the tendency of DPM to produce ice layers at various depths in alpine snowpack simulations and thus to reproduce depth–density variability. It is important to bear in mind that we only use the dual-permeability water scheme of SNOWPACK in DPM and not the other physics of this model; the results produced by the full SNOWPACK model would be different because it has its own formulations for snow mechanical and thermal properties. In particular, DPM relies heavily on the grain size, and it would thus benefit from better representations of the firn's structural properties. Moreover, the primary purpose of the DPM implementation in SNOWPACK is to reproduce the occurrence of ice layers in a seasonal alpine snowpack (Wever et al., 2016), whereas in this study we evaluate its ability to simulate representative firn depth–density profiles over the course of numerous decades.

BK with low water-holding capacity is usually used to mimic preferential flow (Reijmer et al., 2012). However, our findings suggest that, in fact, this more closely represents matrix flow as modelled using the Richards equation. We suggest that, in order to use BK in this way, percolation of some meltwater in the presence of ice layers should be considered.

The lack of variability in density in the modelled profiles cannot only be attributed to inaccuracies in the percolation–refreezing process. This is demonstrated in the example of NASA-SE: the layers of the density peak observed around 1 m depth (Fig. 4a and b) were deposited during the final winter of the simulation (2015–2016). As such, these have only been influenced by the percolation and refreezing of negligible amounts of liquid water. The consistent underestimation of density variability across all schemes indicates that one or several other factors that are not or poorly represented by firn models likely play a crucial role in firn evolution. These factors may include horizontal water flow, prolonged ponding of water in soaked firn close to the surface, variability in fresh snow density, the effects of firn microstructure and impurity content on densification, wind packing and short-term weather fluctuations. Moreover, the validity of the firn model relies on the accuracy of the climatic forcing.

At all sites, interchanging the HL, KM and CROCUS formulations for firn densification generally leads to more variability in the results than using different water flow schemes. This highlights the need to understand the causes of disparity between densification models under various climatic conditions and thus to proceed to further model intercomparison experiments (Lundin et al., 2017). The existing firn-densification formulations are likely not suited for representing densification in conditions of high water contents and high refreezing rates. Our study indicates that firn-densification models could be improved by accounting for the latent heat source as well as the effects of liquid water and of refreezing cycles on firn viscosity and densification rates. For example, the KM and HL densification equations were established for dry firn (Herron and Langway, 1980; Kuipers Munneke et al., 2015a). In the CROCUS scheme, firn viscosity is adjusted according to the water content, but our results show that the modification in the parameterisation is insufficient to reproduce the observed densities at KAN-U. A simple example of the densification schemes in HL and KM not representing reality becomes apparent when applying the percolation–refreezing schemes: in reality, densification is dependent on the overburden stress, but these models use accumulation rate as a proxy for stress. Consequently, in these models the redistribution of mass due to runoff and percolation does not affect the densification rates, despite the effect it has on the firn-column mass. The absence of a preferential flow scheme is often presented as a possible explanation for firn-density overestimation close to the surface (e.g. Gascon et al., 2014; Kuipers Munneke et al., 2014; Steger et al., 2017a). However, our results suggest that simply adding a one-dimensional preferential flow scheme, although physically detailed, to firn-densification models does not solve this issue. The water that is transported quickly from the shallow layers must flow back into the matrix domain at some point. If this occurs in the shallow layers, the density-overestimation issue remains; if this occurs in deep layers it can lead to unrealistic temperature signatures. There are several other possible factors for densification errors at high-melt sites, including exaggerated sensitivity of the model to temperature because the densification model is not suited for conditions with substantial latent heat release.

We compare the results reached at KAN-U with previous firn modelling studies at this site. Langen et al. (2017) also used the CROCUS densification scheme and a water flow scheme conceptually comparable to R1M but with a simplified solving process of the RE. Their results show a density profile entirely at ice density from 4 m depth and an overestimation of the temperature at 10 m depth (+3 to +5 K). In our study, the R1M results show lower deep densities and a 10 m temperature underestimation. These discrepancies can be attributed to differences in (1) model implementation, (2) climatic forcing and (3) details of the water flow scheme. Model resolution is likely of importance here: Langen et al. (2017) used a coarser vertical resolution, disfavouring the formation of thin impermeable layers, allowing liquid water to flow more readily to greater depths and refreeze, thus causing greater density and higher temperature values. Steger et al. (2017b) used the SNOWPACK densification model (recall: we use different densification physics to SNOWPACK in this study) with a bucket scheme, similar to BK, and the same climatic forcing as in this study. Their model output shows a firn column fully compacted to ice at KAN-U. They attribute this to the overestimation of densification rates by the densification model used and also argued that some other effects could be at play, such as the surface density applied. The 10 m temperature is underestimated in their result, likely due to the absence of percolation of water through ice and thus no latent heat release at depth.

Both DPM and R1M exhibit significant sensitivity to the choice of the grain-size formulation (grLK or grA). Modifying this formulation in R1M affects the model results more than changing to the use of BK at all sites apart from FA13 where the magnitudes of change are comparable. For example, the change in grain size at FA13 in R1M leads to a difference of 4 % in FAC and 1.2 K in 10 m temperature, and, in DPM, the patterns of ice layer formation are significantly different at all sites where such ice layers form. This highlights another significant difficulty for percolation schemes: the dependence of water flow on the firn's structural properties. Field evidence demonstrates the crucial role of structural transitions, even at the scale of centimetres, on the behaviour of water flow in firn (Marsh and Woo, 1984; Pfeffer and Humphrey, 1996, 1998; Williams et al., 2010). With respect to this, the advanced flow schemes applied in this study have some limitations. Firstly, the structural properties of grains in the firn layer are poorly constrained by observations. Secondly, the parameterisations linking the structural and hydraulic properties on which R1M and DPM rely were derived from a limited number of laboratory experiments. These are typically performed at a very small scale (e.g. shallow snow columns with diameter of 5 cm in the experiments of Yamaguchi et al., 2012) and are mostly based on homogeneous snow in terms of grain size and temperature (Yamaguchi et al., 2012; Katsushima et al., 2013). The much larger scale of the GrIS firn layer, the spatial and temporal heterogeneity of its structural properties, and its climatic and glaciological settings render the validity of these idealised parameterisations questionable. Finally, the density dependence of the parameters makes them sensitive to errors in the densification process. Thus, a better knowledge of firn structural properties would only be profitable to water flow schemes if we have a clear understanding of the link between snow structure and its hydraulic properties and vice versa.

Another major limitation of the implementation of physically detailed liquid water flow schemes in one-dimensional firn models is the fact that water flow is in reality three-dimensional. Water can flow horizontally on top of buried ice lenses or on thin, near-surface ice crusts caused by daily refreezing (Marsh and Woo, 1984; Pfeffer and Humphrey, 1996). Even at depths greater than 10 m, large masses of liquid water can persist through the winter and move laterally over considerable distances (Humphrey et al., 2012). In one-dimensional models, the key to solving this issue is to accurately partition between vertical percolation and lateral flow; this likely requires a better approach than the lateral runoff parameterisation we implement here. As an example, at FA13 all three water-transport schemes overestimate the density in the upper 10 m except in the last winter's accumulated layers (0–2 m depth), where there is a good agreement with observations (Fig. 10a). This suggests a consistent overestimation of the summer meltwater refreezing and underestimation of lateral runoff. Also, the need to use a limit for the PFA water content demonstrates that some processes not represented in the model must regulate its water volume; these are likely lateral movement driven by hydraulic pressure gradients and connections with englacial and subglacial hydrological systems. A one-dimensional preferential flow scheme aims to correctly partition the water input between matrix flow and fast preferential flow; there are several other difficulties with this approach. These include accurately determining how deep the water can be transported by preferential flow, how much water refreezes and how much is stored as liquid water, and the amount of lateral flow at different depths in the firn column. Similarly, the same considerations apply to liquid water flowing in from upstream grid cells. The representation of preferential flow physics thus requires improvements to address some deficiencies. We highlight (1) that too much water is transported through preferential flow or at exaggerated depths as demonstrated by the consistent overestimation of 10 m temperature; (2) the need to consider lateral flow as the tendency to underestimate FAC shows; and (3) the large sensitivity of the ice layer formation process to grain size, which should be overcome or should be addressed by further observational studies of grain metamorphism in firn.

Our results also suggest that more observations of firn-temperature variability in time and depth would likely be useful for the evaluation of existing flow schemes and the development of new ones. Modelled temperature profiles show both the depth and volume of refreezing due to the release of latent heat. Moreover, deep meltwater refreezing causes marked and long-lasting temperature increases due to insulation from the overlying firn, and temperature measurements in depth can be powerful indicators of the occurrence of deep percolation and refreezing events. On the other hand, comparisons between modelled and observed density profiles are strongly affected by the choice and accuracy of the densification formulation, the variability of surface density, several other factors influencing model outputs mentioned above and possible uncertainties in field measurements. Such uncertainties are related to the strong spatial variability of firn structure (Marchenko et al., 2017), which can be observed by comparing density profiles of cores drilled at nearby locations. However, the modelled temperature profile also depends on the accuracy of the climatic forcing, of the heat-transport scheme and of the thermal-conductivity parameterisation. The latter is a function of density; thus erroneous depth–density profiles inevitably lead to an inaccurate heat-transport process. As an example of the influence of various sources of errors, the warm 10 m temperature bias of all the schemes at NASA-SE (Table 3) is unlikely to be only due to the percolation–refreezing process. RACMO2.3p2 has a slight negative precipitation bias at this site (−0.21 m w.e. yr⁻¹; see Sect. S1), which likely leads to underestimated advection of cold snow down into the firn column and could partly explain the firn models' temperature bias. Similarly, air temperature is overestimated here by RACMO (+1.5 K), which could enhance the problem. It is possible that the consistent over- and underestimation of FAC at FA13 and at KAN-U could be due to snowfall over- and underestimation respectively by RACMO; however, no observations of precipitation at these sites are available for comparison.

6 Conclusions

We implemented three liquid water schemes of different levels of physical complexity in a firn model using a fine vertical resolution: a bucket scheme, Richards equation in a matrix flow scheme and Richards equation in a preferential flow scheme. To our knowledge, this is the first study to apply the Richards equation as well as a preferential flow scheme in firn-densification simulations on the GrIS.

Our three liquid water flow schemes predict significantly different vertical patterns of refreezing and consequently modelled densities, firn air content values and 10 m temperatures. The preferential flow scheme effectively evacuates meltwater from the surface layers and leads to underestimation of firn air content and overestimation of 10 m temperatures. Compared to the preferential flow scheme, the single-domain Richards equation scheme generally showed biases of the opposite signs and of much lower magnitudes in both FAC and 10 m depth temperature, suggesting it slightly underestimates percolation depths. The simpler bucket scheme predicted refreezing rates, firn air contents and 10 m temperatures similar to those obtained by the single-domain Richards equation; by adjusting its water-holding capacity and impermeable density parameters, it could produce the same results. Using the Coléou and Lesaffre (1998) parameterisation for the water-holding capacity in the bucket scheme led to underestimation of percolation depths. The bucket scheme with lower water-holding capacity and the single-domain Richards equation scheme predicted firn air contents and 10 m temperatures in closest agreement with observations. However, the preferential flow scheme was found to perform better than the simpler flow schemes in reproducing the density variability with depth and the water-saturated conditions at the bottom of the firn column at a site of a perennial firn aquifer.

We identified the multidimensionality of liquid water flow as the prominent challenge for water percolation schemes. Because firn models are currently one-dimensional, an accurate partitioning between horizontal and vertical flow is likely to be at least as difficult and as important as the separation between slow matrix and rapid preferential flow. Other difficulties related to water flow representation include the uncertainties in firn hydraulic properties and in firn micro- and macrostructure on the GrIS. This is further demonstrated by our results showing the sensitivity of the Richards-equation-based schemes on the grain-size formulation. However, the absence of any large-scale field observations of water flow in firn makes it difficult to constrain its implementation and to validate model behaviour. Here, we aimed to assess the applicability of flow schemes developed for seasonal snow models to the Greenland firn. Whilst we did apply some modifications to account for the differences between snow and firn, we suggest that more modifications are likely required since the spatial scales and the structural characteristics of seasonal snowpacks and firn are different.

There is no large-scale detailed observation available of liquid water content and percolation pattern during melting events in firn. This renders the validation of a particular flow scheme difficult, and validation relies on temperature and density profiles. However, there are a number of effects that influence firn density and temperature, all potentially contributing to mismatch between modelled and observed values. As an example, the density variability in depth was largely underestimated regardless of the flow scheme. This suggests that there are uncaptured complexities in the percolation and refreezing mechanisms that need to be better taken into account, such as water ponding and the effect of grain shape on water flow. Furthermore, this shows that firn-model development must focus on including complex processes currently poorly or not represented, such as surface-density variability, wet firn densification and firn structural effects on densification. We note also that our assessment considers only a single profile from each site, captured at one time point, and further work is required to assess the impact of liquid water percolation on transient firn evolution. A comprehensive exploration of the various firn models and their parameter spaces could help identify priorities for further model developments based on minimising data–model mismatch and overall uncertainty. In line with this, we showed that output from three common firn-densification models shows greater variability than the output from a single densification model using the different flow schemes. In order to capture the multiple impacts of liquid water on firn densification, future models require an improved liquid water flow scheme, accurate boundary conditions and formulations developed explicitly to simulate densification of wet firn.