The modelled liquid water balance of the Greenland Ice Sheet

Recent studies indicate that the surface mass balance will dominate the Greenland Ice Sheet’s (GrIS) contribution to 21st century sea level rise. Consequently, it is crucial to understand the liquid water balance (LWB) of the ice sheet and its response to increasing surface melt. We therefore analyse a firn simulation conducted with SNOWPACK for the GrIS and over the period 1960–2014 with a special focus on the LWB and refreezing. An indirect evaluation of the simulated refreezing climate with GRACE and firn temperature observations indicate a good model performance. Results of the LWB analysis reveal 5 a spatially uniform increase in surface melt during 1990–2014. As a response, refreezing and runoff also indicate positive trends for this period, where refreezing increases with only half the rate of runoff, which implies that the majority of the additional liquid input runs off the ice sheet. However, this pattern is spatially variable as e.g. in the southeastern part of the GrIS, most of the additional liquid input is buffered in the firn layer due to relatively high snowfall rates. The increase in modelled refreezing leads to a general decrease in firn air content and to a substantial increase in near-surface firn temperature in some regions. 10 On the western side of the ice sheet, modelled firn temperature increases are highest in the lower accumulation zone and are primarily caused by the exceptional melt season of 2012. On the eastern side, simulated firn temperature increases more gradually and with an associated upward migration of firn aquifers.


Introduction
The mass balance (MB) of the Greenland Ice Sheet (GrIS) has been negative since the early 1990s ( Van den Broeke et al., 2016). In addition to increased ice discharge through the acceleration of marine-terminating outlet glaciers, the ice sheet is losing mass through increased surface melt and associated meltwater run-off. The latter process has recently become the dominant contributor to mass loss from the ice sheet (Enderlin et al., 2014). The increase in meltwater runoff and associated decrease in the surface mass balance (SMB) is attributed to processes on various spatial and temporal scales, such as the polar amplification (Bekryaev et al., 2010) and the darkening of the GrIS . Some of these processes are further promoted by the hypsometry of the ice sheet Van As et al., 2017). An accurate quantification of the liquid water balance (LWB) of the ice sheet is important, as it determines how much of the liquid input at the surface ultimately reaches the ocean and contributes to sea level rise. A key parameter of the LWB is meltwater storage in the firn (Rennermalm et al., 2013a) by refreezing and liquid water retention. Previous studies suggest that modelled refreezing strongly depends on the model formulation (Reijmer et al., 2012;Steger et al., 2017) and that it exhibits the largest inter-model variation of all SMB components (Vernon et al., 2013). Besides the instantaneous effect of retaining liquid water, refreezing also co-determines the future potential of firn to absorb melt, as it reduces the porosity of the firn  and releases large amounts of latent heat Cox et al., 2015), which decreases the firn's cold content.
The hydrology of the GrIS is a complex system, which involves various ill-constrained processes ( Fig. 1). At the surface, liquid input is determined by rainfall, evaporation/condensation and melt. In areas where the ice sheet is covered by snow and/or firn, liquid water is able to percolate vertically. These snow/firn layers may act as a buffer for run-off if liquid water either refreezes  or remains in its liquid state in perennial firn aquifers . Such aquifers typically form at lo- cations with relatively high amounts of snow accumulation (Kuipers Munneke et al., 2014) and are thus particularly abundant along the south-eastern and north-western margins of the ice sheet Miège et al., 2016). A recent study (Poinar et al., 2017) revealed that some aquifers likely drain into crevasses. To what degree the water refreezes there or reaches the bed of the ice sheet remains largely unknown. Along the south-western and north-eastern margins of the ice sheet, firn aquifers are less abundant. In these areas, percolating water typically refreezes in the firn or runs off over the ice surface. A study by  suggests that horizontal ice layers could inhibit vertical percolation and render underlying pore space inaccessible for liquid water. The water would hence be forced to flow laterally above such obstacles -either as surface run-off or within the firn.
In the bare ice zone, hydrological processes are better understood: liquid water flows along surface rivers and may accumulate in supraglacial lakes (Arnold et al., 2014) or enter the subglacial system via moulins or crevasses. The amount of water stored in supraglacial lakes, which may be buried with snow in winter and hence retain liquid water perennially, is thereby rather small compared to the magnitude of supraglacial river fluxes (Smith et al., 2015;Koenig et al., 2015). Liquid water flowing into moulins or crevasses enters the en-and subglacial (Lewis and Smith, 2009;Lindbäck et al., 2015) hydrological system of the ice sheet. Here, water may refreeze, accumulate in subglacial lakes or flow along channels to the margins of the ice sheet. The relevance of en-and subglacial water storage is currently rather uncertain. Rennermalm et al. (2013b) suggests that for a watershed in south-western Greenland, up to 54 % of meltwater may be retained during one season. It is however possible that this residual is partly caused by uncertainties in, for example, watershed delineation (Rennermalm et al., 2013b) and inter-basin piracy (Lindbäck et al., 2015). A more recent study for a similar catchment yielded little evidence for meltwater storage in en-and subglacial environments . In short, the hydrology of the GrIS represents a complex system of pathways that transport meltwater from the surface of the ice sheet to the ocean (Chu, 2014).
In this study, we quantify the components of the LWB from the GrIS surface to the firn-ice transition, using a state-of-the-art snow/firn model. The upper boundary conditions for the model are provided by the regional atmospheric climate model RACMO2.3 . Potential englacial (below the firn-ice transition) and subglacial liquid water retention is not considered as we only model the upper part of the ice sheet (∼ 40-200 m). The primary goal is to quantify the spatial magnitude of the LWB components and assess how these mass fluxes evolved over the last decades. We evaluate the spatial and seasonal occurrence of refreezing and the impact of this process on firn density and temperature. Finally, we analyse how the horizontal extent of firn aquifers, which act as perennial storage for liquid water, evolves with time. The following section provides a brief description of the model and the observational data used in this study. Subsequently, we discuss the comparison of model output with remote sensing data (GRACE) and in situ measurements (firn temperatures). Section 4 contains the results of the LWB evaluation and a more detailed analysis of refreezing, run-off and changes in different firn properties.
2 Definitions, model and data

Definitions
In this study, we investigate the LWB of the upper part of the ice sheet, namely the snow/firn layer. This layer ranges from the surface down to the firn-ice transition, usually defined at the density of 830 kg m −3 . If percolating water reaches the bottom of the model domain, it is considered to leave the ice sheet as run-off. Potential en-and subglacial storage of liquid water at greater depth is neglected. The LWB of the firn layer is defined as where M ret is the retained liquid mass, RA, EV and ME are surface mass fluxes of rainfall, evaporation and meltwater respectively, RF is internal refreezing and RU is run-off at the bottom of the model domain. In this study, the term evaporation refers to phase changes of water from liquid to gaseous (evaporation) and vice versa (condensation). The SMB used in this study equals the climatic mass balance (Cogley et al., 2011); i.e. it includes subsurface processes of liquid water retention and refreezing, and is defined as where SN is snowfall, SU sublimation (and resublimation) and SD deposition or erosion by snow drift. The SMB is linked to the LWB through the components rainfall, evaporation and run-off. The Greenland MB derived to validate the modelled SMB with GRACE data is defined as where SMB GrIS and SMB PIC are the SMB of the glaciated area (GrIS and peripheral ice caps/glaciers), D is ice discharge across the grounding line from marine-terminating glaciers and M ts is the tundra snow mass.

Model data
Snow/firn on the GrIS and the peripheral ice caps/glaciers is modelled with the SNOWPACK model (version 3.30). SNOWPACK has been applied in several studies (Groot Zwaaftink et al., 2013;Van Tricht et al., 2016;Steger et al., 2017) to simulate snow and firn in polar regions. The model contains an overburden-dependent densification scheme and simulates the evolution of different microstructural snow properties, which are linked to thermal and mechanical snow quantities (Bartelt and Lehning, 2002;Lehning et al., 2002a, b). We run SNOWPACK on an 11 km horizontal grid and with the same ice mask (Fig. 2) as used in the regional atmospheric climate model RACMO2.3 . At the snow-atmosphere interface, SNOWPACK is forced with mass fluxes (precipitation, evaporation/sublimation, snow drift and surface melt) and with skin temperature from RACMO2.3. Skin temperature is the temperature of an infinitesimally thin layer without heat capacity and is representative of surface temperature. The capability of RACMO2.3 to accurately simulate present-day surface climate on the GrIS was illustrated in an extensive www.the-cryosphere.net/11/2507/2017/ The Cryosphere, 11, 2507-2526, 2017 evaluation by Noël et al. (2015). Vertical water percolation is simulated with a bucket scheme (Bartelt and Lehning, 2002;Wever et al., 2014) and the irreducible water content follows the formulation of Coléou and Lesaffre (1998). We do not consider heterogeneous percolation (Wever et al., 2016;Marchenko et al., 2017) in our simulation due to an insufficient spatial coverage of observational data to calibrate such routines for the entire ice sheet and/or the too-expensive computational demand. Neglecting heterogeneous percolation causes refreezing to occur mostly in the upper snowpack, where temperature and porosity are determined by the recent climate. Lateral flow of run-off is also not considered in our simulation. Fresh snow density is prescribed with an empirical parameterisation that depends on mean annual surface temperature . The enhanced near-surface snow compaction due to strong winds, which is implemented in SNOWPACK for Antarctic simulations (Groot Zwaaftink et al., 2013), is switched off, because the applied fresh snow density parameterisation already accounts for this effect. A more detailed description of the model set-up and the applied spin-up procedure is stated in Steger et al. (2017), where the same SNOWPACK run was used.

Observational data
To derive a MB for Greenland, we use ice discharge data from Enderlin et al. (2014) and a GRACE gravity field solution for Greenland (Groh and Horwath, 2016). The ice discharge data contain annual estimates of ice discharge from 178 marine-terminating glaciers wider than 1 km and are available for the period 2002-2012. Following Van den Broeke et al. (2016), we neglect seasonal variations in ice discharge and assume that all intra-annual variation in the MB is induced by components of the SMB or by tundra snow. The GRACE data we apply are based on the monthly GRACE solution ITSG-Grace2016 (Mayer-Gürr et al., 2016) and is available between mid-2002 and mid-2016. We computed the MB for the overlapping period 2003-2012 where data are available from all sources throughout the year. We use firn temperatures that were recorded along a 2700 km transect in north-western Greenland (Fig. 2), referred to as the NW GrIS transect, to evaluate our simulation. Shallow borehole temperature measurements were conducted at 14 sites between 1952 and 1955 (Benson, 1962) and repeated in 2013 (Polashenski et al., 2014). The former measurements were taken at a range of 3 to 16.75 m depth (predominantly at 8 m) and were corrected for seasonal influences to obtain an intercomparable, mean annual 10 m temperature. The measurements in 2013 were recorded at a depth between 5 and 12 m (mainly at 8.5-12 m) and were corrected with the same methodology (Polashenski et al., 2014).

Model evaluations
Although modelled refreezing cannot directly be evaluated with observations, Steger et al. (2017) made a comprehensive assessment of modelled firn density, which is the combined result of dry compaction and refreezing. Results show a reasonable performance of SNOWPACK but a general overestimation of densities in the percolation zone. This bias is likely the result of overestimated near-surface refreezing caused by neglecting heterogeneous water percolation, an overestimation of fresh snow density and/or errors in the atmospheric forcing (Steger et al., 2017). In this study, we use additional observations to evaluate the ability of SNOWPACK forced by RACMO2.3 to simulate the LWB and particularly the refreezing climate of the GrIS. Due to the lack of direct refreezing observations, we assess the model's performance in terms of refreezing indirectly by comparing the modelled spatially integrated SMB and local snow/firn temperatures to the observations.

Model evaluation using GRACE
Due to the large footprint of GRACE, the signal also contains mass variations from Greenland's peripheral ice caps and glaciers and from tundra hydrology, primarily from seasonal snow cover. These signals are thus included in Greenland's MB as explained in Sect. 2.1 Tundra snow cover is not simulated by SNOWPACK but the signal is taken from RACMO2.3 output. In RACMO2.3, seasonal snow is simulated with a single-layer model that does not allow for refreezing and liquid water retention in the snow (Van den . All surface melt is hence immediately transferred to run-off. A comparison between the derived cumulative MB and GRACE is provided in Fig. 3a. The MB is computed by taking the simulated SMB over the glaciated area either from RACMO2.3 or SNOWPACK. Both cumulative MBs indicate an excellent agreement with GRACE (R 2 > 0.99). In terms of linear trends, SNOWPACK agrees better with GRACE due to higher modelled refreezing fractions and thus lower amounts of run-off from the ice sheet. The detrended mean seasonal cycles (Fig. 3b) indicate a good agreement in winter and spring, when changes in cumulative SMB are mainly caused by accumulation of solid precipitation on the glaciated area and the tundra. From May onwards, the derived MBs show an earlier and steeper decrease compared to the GRACE signal. The minima in the MBs occur both earlier and with higher magnitudes than in GRACE, where SNOWPACK performs slightly better due to smaller amounts of modelled run-off. These findings are consistent with earlier studies Alexander et al., 2016), in which the average seasonal cycle of the MB and GRACE were compared. A likely contributor to this mismatch is the neglect of the time it takes run-off to reach the ocean. Van that the monthly error between detrended modelled SMB and GRACE on a GrIS-wide scale could be minimised by delaying simulated run-off by 18 days. A study for a catchment in south-western Greenland revealed that transit times up to 10 days are required to align the modelled surface run-off and observed river hydrograph optimally (Van As et al., 2017). Another uncertainty arises from modelled tundra snow cover and tundra hydrology. The mean seasonal amplitudes of the detrended modelled MBs originate ∼ 30 % from winter accumulation and summer melting of seasonal snow over the tundra (Fig. 3b). A too-early snow ablation in the tundra could hence also contribute to the bias between MBs and GRACE. This assumption is supported by a comparison of the simulated snow cover fraction (SCF) with MODIS/Terra Snow Cover data (Hall and Riggs, 2016), which revealed a too-early decrease in modelled SCF in most basins (not shown). Potential causes for this bias are the neglect of refreezing and liquid water retention in the relatively simple RACMO2.3 snow model and the poor representation of tundra topography at a horizontal resolution of 11 km. Additionally, heterogeneous snow distribution on a subgrid scale could also contribute to the bias (Aas et al., 2017). Finally, run-off may also be retained in the hydrological system of the tundra by refreezing in soil, ponding on frozen ground (Johansson et al., 2015), accumulating in surface lakes (Mielko and Woo, 2006) and storage in terrestrial aquifers. All these processes are currently not represented in our model framework.

Model evaluation with firn temperature measurements
ERA40 reanalysis data, which force SNOWPACK via RACMO2.3, are available from 1958 onwards. SNOWPACK output is therefore not available for the years 1952-1955, when the first firn temperature data set along the NW GrIS transect was collected. We assume, however, little change between the 1952 and 1955 observed temperatures and compare them to modelled firn temperatures from 1960, after model spin-up has moderated. We coincidently compare the 2013 temperature observations to modelled data. Figure 4 shows that SNOWPACK forced by RACMO2.3 slightly overestimates firn temperatures in the higher elevations of the NW GrIS transect for both periods. For the earlier period, this bias may be partly caused by the spin-up procedure of the model, where the model is looped over the reference period (1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979) to generate the initial firn profile (Steger et al., 2017). This means that surface temperature evolutions before this reference period are not considered. The bias for the later period is more difficult to explain in the absence of continuous firn temperature measurements and firn density records. The spatially incoherent firn temperature change (between B 4-225 and B 4-000) in the observations is not reproduced by SNOWPACK, which simulates a uniform temperature increase of ∼ 0.3 • . This incoherency in the observations may be partly explained by uncertainties in the measurements caused by errors in the sensor calibration and uncertainties in the applied correction used to retrieve 10 m firn temperature from shallower measurements (Polashenski et al., 2014). Between locations B 2-175 and B 2-070, there is a ∼ 1.6-2.7 • C warming in the observations between 1952-1955 and 2013, likely caused by latent heat release due to refreezing. This temperature increase is larger than the modelled, spatially rather uniform warming of ∼ 0.5 • C. Possible explanations for this bias are the underestimation of meltwater production at the surface or a too-shallow refreezing depth, which enables the released heat to be conducted upwards to the surface and escape to the atmosphere through emission of longwave radiation. In SNOWPACK, percolating water is not allowed to pass unhindered through layers with refreezing capacity, where in reality, liquid water may move to greater depth through heterogeneous percolation Marchenko et al., 2017). At site B 1-010, SNOWPACK simulates a local maximum in firn warming, in agreement with observations. Here, RACMO2.3 simulates a doubling of the liquid water input between the two periods considered. However, the magnitude of warming in SNOWPACK is somewhat smaller (4.1 vs. 5.7 • C), which may again be linked to the neglect of heterogeneous percolation. Along the entire transect, modelled increases in solid precipitation are spatially rather uniform and small (∼ 0.02 m w.e. a −1 ) and, therefore, likely less relevant for explaining changes in firn temperature.
Other snow/firn temperature records are available from the Greenland Climate Network (GC-Net; Steffen and Box, 2001) and for the western percolation zone Charalampidis et al., 2016). The latter two data sets, which cover the periods 2007-2009 and 2009-2013, also (Steger et al., 2017), forced by RACMO2.3, do not reproduce the strong warming observed at these locations. The reason is the overestimation of the bare ice zone on the western GrIS by the IMAU-FDM and SNOWPACK; i.e. the models are incapable of simulating the subsurface warming due to a deficiency of pore space for refreezing. Compared to IMAU-FDM, the overestimation of this zone is less pronounced in SNOWPACK owing to a different densification scheme, which is more accurate for relatively warm conditions (Steger et al., 2017). Inferring the exact extent of the bare ice zone from remote sensing data, by using the different spectral surface properties of snow and ice, is complicated by the formation of near-surface ice layer (Machguth et al., 2016) above porous firn. At higher elevations in western Greenland, SNOWPACK does simulate a pronounced warming of the firn layer. Unfortunately, the subsurface temperature data recorded at GC-Net stations located in this area (DYE-2, Crawford Point 1 & 2 and GITS) suffer from large data gaps and/or unphysical high-frequency fluctuations caused by sensor deterioration (K. Steffen, personal communication, 2017). The data are thus of insufficient quality to verify these changes.
To address the above-mentioned model bias in overestimating the bare ice zone, we briefly assessed fresh snow density, which is a rather uncertain factor in our simulation. The empirical relation  we use to obtain this quantity was derived with samples from the dry snow zone and is subsequently extrapolated to lower elevation on the ice sheet. Snow/firn density profiles from a transect on the western GrIS  allow a comparison between observed and modelled nearsurface densities: averaging over the upper 50 cm and all samples yields a value of ∼ 345 kg m −3 for April (i.e. before the onset of seasonal surface melt). For these locations, our fresh snow density parameterisation returns a mean density of ∼ 405 kg m −3 . The parameterisation, which accounts for near-surface densification due to wind and vapour fluxes, clearly overestimates fresh snow density for this region. A comparison of our fresh snow density parameterisation with near-surface snow density samples obtained on the northern GrIS and for spring  supports the assumption that the applied parameterisation yields too-high densities for comparably warm climate conditions (not shown). To test SNOWPACK's sensitivity to initial snow densities, an experiment with a lower, spatially uniform fresh snow density of 320 kg m −3 was carried out for the west-ern GrIS transect. The selected initial density is comparable to what the recently published parameterisation of Langen et al. (2017) yields for this transect. With this model setting, the mismatch between the observed and modelled bare ice zone extent (and thus the firm warning) was reduced for this specific region (not shown). This model inaccuracy should be addressed in future by testing available or newly derived fresh snow density parameterisations with SNOWPACK for various climate conditions on the GrIS.

Climatology of the liquid water balance
The evaluations of mass changes and firn temperatures with observations presented in the previous sections inspire suf-ficient confidence to use SNOWPACK firn data for a description of the LWB of the GrIS. First, we discuss the mean fields and temporal evolution of the LWB components during the simulation period . Subsequently, refreezing, one of the key components of the balance, and its dependency and influence on firn density and temperature is discussed. And finally, we analyse the temporal evolution of  Figure 5 shows the temporally averaged (1960-2014) LWB components for the GrIS and the peripheral ice caps and glaciers. Mean fluxes of rainfall and evaporation are typically at least 1 order of magnitude smaller than melt, run-off and refreezing. Changes in the retained liquid mass (dM ret /dt) are even smaller than components on the right-hand side of Eq.

The liquid water balance
(1), particularly when integrated over basins, and are thus not presented. Rainfall rates are particularly significant along the southern margin of the ice sheet and in the western ablation zone. For the north-eastern part of the GrIS, the contribution of rainfall to the LWB is small and liquid water input at the surface is dominated by melt. The highest melt rates on the GrIS occur along the western ablation zone with a maximum of 129 Gt a −1 for Basin 6. The mean spatial runoff pattern is comparable to the one of melt but attenuated by the buffering effect of refreezing. Run-off also peaks in Basin 6 with a value of 85 Gt a −1 , which accounts for a third of the total GrIS run-off. Averaged over the entire ice sheet, SNOWPACK simulates that almost half (47 %) of the liquid water input at the surface refreezes in snow or firn. This fraction has a high spatial variability and is relatively low for the north-eastern basins and for Basin 6, where precipitation is low and bare ice extent relatively large. As a result, refreezing rates in these regions peak more inland in the lower accumulation zone just above the equilibrium line. Refreezing in the ablation zone is, in terms of absolute liquid water retention, only relevant on intra-annual scales. The highest overall refreezing fractions, up to 75 %, are modelled along the wet south-eastern margin of the ice sheet (Basin 4). Time series of the four most relevant LWB components (melt, run-off, refreezing and rainfall) for the eight basins show no distinctive trends for the first half of the simulation period  but do exhibit large interannual variability, particularly for surface melt (Fig. 6). For the second half , however, there is a statistically significant in-crease in melt in all basins (Table 1). Changes are particularly large for basins 5 and 6, where melt increases by 0.36 and 0.38 m w.e. a −1 , respectively. The dominant cause for these large changes is the comparably high increase in melt in the ablation area of the GrIS, especially in the south-west. Modelled snow melt in the ablation zone is particularly sensitive to temperature increases due to the albedo difference between snow and ice, where bare ice with a lower albedo is more rapidly exposed through accelerated melting of snow. The lowered surface albedo subsequently enhances melting of bare ice. A secondary cause is the relatively flat hypsometry of these basins, where 58 % and 47 % of the area is below 2000 m a.s.l. (compared to 39 % for the GrIS). Rainfall, as a further contributor to liquid input, does not exhibit a significant trend for the majority of the basins. Linear trends are comparably high for Basin 5 (1.22 mm w.e. a −2 ) and Basin 6 (0.43 mm w.e. a −2 ) but statistically insignificant. Remarkably, the north-western Basin 8 is the only region with a significant positive trend in rainfall of 0.56 mm w.e. a −2 . This increase is not caused by a significant change in total precipitation but by a significant increase of the rainfall fraction in this area. For all basins, melt rates peak in 2012 when the GrIS experienced unprecedented surface melt both in spatial extent (Nghiem et al., 2012) and magnitude. The exposure of relatively high-elevated regions with cold and porous firn to surface melt is the main reason that refreezing also peaks in all basins during this year. In response to the positive trends in melt, run-off also exhibits a significant increase in all basins between 1990 and 2014 ( Table 1). The increase in runoff accounts in most basins for more than half of the increase in melt (∼ 55-80 %); i.e. most of the additional melt is not stored in the firn layer but is running off the ice sheet. As for melt, the south-western basins, 5 and 6, show the strongest increase per area. An exception is Basin 4, where run-off increases with only ∼ 30 % the rate of melt. In terms of refreezing and refreezing fraction, the response of the basins to increased surface melt is spatially less uniform: the majority of the basins do not indicate a significant trend in refreezing. This means that, for example, in basins 1 and 2, most of the additional melt is not absorbed in the firn but runs off, similar to what happens to northern ice caps not connected to the main ice sheet . Basin 4, which has the highest overall mean refreezing fraction (75 %), is an exception. Refreezing in this basin shows a distinctive positive trend. This is linked to the high amounts of solid precipitation in this basin, which provide enough pore space to absorb the increase in surface melt. Refreezing is also significantly increasing in the north-western basins 7 and 8 but with a lower trend than run-off. Significant trends in the refreezing frac-tion are only apparent in Basin 1 and particularly in Basin 8, where the fraction decreases by ∼ 16 % in 25 years.
For the entire GrIS, melt, run-off and refreezing indicate significant positive trends between 1990 and 2014. The increase in run-off is roughly twice that of refreezing, which leads to a significant decrease in the GrIS-integrated refreezing fraction of ∼ 9 % over the 25 years (Table 1). The different responses of the eight basins to increasing surface melt are related to refreezing, which in turn is linked to firn porosity and temperature.

Refreezing and latent heat release
Refreezing is a process that strongly depends on local climate, particularly on surface temperature, which is the main driver for melt, and therewith on seasonality and elevation ( Fig. 7). At the beginning of the melt season, modelled refreezing primarily occurs in the lower parts of the ice sheet, where the melt onset is earliest and meltwater percolates into the cold winter snow layer. For basins 3-7, low-level refreezing peaks in mid-to late May while for the northern basins, 1, 2 and 8, the maximum occurs in mid-June. During the course of the melt season, the lower regions are gradually depleted of spore space or cold content and the area of peak refreezing moves upward. For the majority of the basins (e.g. basins 1, 2 and 7), the availability of pore space is the limiting factor for refreezing at lower elevations (Fig. 7). Particularly for basins 4 and 5, however, this is not the case. In these basins, refreezing at lower elevations persists throughout the melt season but with lower rates than in spring due to a gradual decrease in the firn cold content. Even if the entire firn column has become temperate, modelled refreezing persists due to the diurnal temperature cycle, which periodically refreshes the near-surface cold content during night. This underlines the importance of using atmospheric forcing data that resolve variations on subdaily timescales. Peak refreezing rates pass the equilibrium line altitude in June for western basins, 6-8, and in July for northern basins, 1 and 2. For basins 3-5, it is not possible to define a mean equilibrium line altitude because these basins have a very narrow ablation zone and at 11 km resolution, many model grid cells close to sea level have a positive SMB due to high accumulation rates. In July, peak refreezing moves beyond the run-off line in all basins. Basins 4 and 5 reveal a percolation zone that stretches over a relatively large vertical extent, with substantial refreezing as high as 500-750 m above the modelled run-off line. A further notable feature of refreezing is its seasonal asymmetry in comparison to melt (Fig. 7). Melt peaks in July in all basins and the seasonal increase and decrease are rather symmetric around this maximum. Refreezing, on the other hand, peaks at the beginning of the melt season in all basins and at all elevations. As mentioned above, this is mainly caused by the decrease in either pore space or cold content during the melt season. A similar feature was found by Cullather et al. (2016) in the regional climate model MAR, when seasonal run-off is plotted as a function of melt area. Run-off was thereby found to be higher in the second half of the melt season. As Cullather et al. (2016) states, this is an important finding for refreezing or run-off parameterisations that do not take seasonality into account.
Refreezing rates increase in all basins during 1990-2014 ( Fig. 6) but not always at a significant level (Table 1). Generally, increases in refreezing are restricted to elevations above 1000 m a.s.l. in all regions (Fig. 8a). The peak of this increase is around 1500 m a.s.l. for the northern basins 1, 2 and 8 and at higher elevations for the more southerly located regions. This increase is primarily caused by a gradual expansion of www.the-cryosphere.net/11/2507/2017/ The Cryosphere, 11, 2507-2526, 2017 the melt area to higher elevations, which causes melt and refreezing to occur in formerly dry snow/firn. The increase in refreezing induces both a corresponding decrease in the modelled firn air content (Fig. 8b) and an increase in firn temperature (up to 4 • C, Fig. 8c) due to latent heat release. Particularly for basins 4 and 5, firn air content also decreases in areas below 1000 m a.s.l. This reduction is not related to changes in refreezing but rather caused by increases of melt and the subsequent transformation of formerly porous firn to bare ice. In contrast to other basins, Basin 2 reveals a relatively constant firn temperature increase at lower elevations. This increase is not only caused by the rather small increase in refreezing and the associated latent heat release but also by an enhanced vertical heat flux from the surface through an increase in surface temperature (Fig. 8d). Surface temperatures changes show a distinct spatial variability, with the largest increases occurring in the north-eastern part of the ice sheet, where temperature increases by more than 1.5 • C. The exceptional melt season of 2012 has an even stronger influence on firn temperatures according to our model simulation: Fig. 9a shows the corresponding refreezing anomaly for this year and Fig. 9b the resulting increase in firn temperature. Almost the entire ice sheet experienced exceptional refreezing rates above the equilibrium line, particularly in the southern area of the GrIS where refreezing anomalies up to +0.8 m w.e. a −1 are modelled. The increase in firn temperature largely reflects the refreezing anomaly, with the strongest warming (6 • C and higher) being modelled in the south-western percolation zone of the ice sheet. Unfortunately, no observational data are available to confirm the pronounced warming. The closest available record is from the KAN_U automatic weather station of the Greenland Analogue Project (GAP) and the Programme for Monitoring of the Greenland Ice Sheet (PROMICE), which is located in the lower accumulation zone of Basin 6. There, firn temperature increased by approximately 4.7 • C during 2012 . To discuss changes in the vertical firn properties over the simulation period in more detail, we present cross sections of firn density, temperature and volumetric water content along a southern GrIS transect (Fig. 2) for the beginning of the simulation period (April 1960, Fig. 10) and as relative changes for the end (April 2014, Fig. 11). In 1960, the transition from bare ice to porous firn is modelled around station KAN_U on the western side of the ice sheet. On the eastern side, no bare ice zone has formed due to the high accumulation rates in this region. Increasing accumulation rates from west to east induce the downward bending of high-density layers east of the ice sheet divide (Fig. 10a). The larger amount of pore space on the eastern side permits larger refreezing fractions, which leads, through release of latent heat, to temperate firn condition close to the margin of the ice sheet and to the formation of a firn aquifer (Fig. 10c).
During the 55 years of the simulation, the firn layer along this transect experienced some distinctive changes: nearsurface density increased both on the eastern and western side of the ice sheet (Fig. 11a) with a shift of the transition between bare ice and porous firn to higher elevations on the western side. For 2012, SNOWPACK simulates a bare ice profile for KAN_U while a retrieved density profile for this year revealed layers with porous firn . A potential cause for this overestimation in firn density is discussed in Sect. 3.2. Modelled firn temperature increased substantially on both sides of the GrIS but with different patterns (Fig. 11b). In the west, 15 m firn temperature in the percolation zone is relatively stable until 2010 and abruptly increases afterwards by ∼ 10 • C, particularly due to the exceptional melt season of 2012. Refreezing during this year induces a substantial warming of the firn down to a depth of 40 m. On the eastern side, the initial 15 m firn temperature is higher by ∼ 5 • C and the simulated warming of the firn is more gradual. A major reason for the less pronounced and shallower firn temperature increase on the eastern side of the ice sheet is the temporal distribution of liquid input in 2012 (inset panel Fig. 11c). On the eastern side, liquid input at the surface is rather evenly distributed throughout the melt season. On the western side, there are several distinctive peaks with liquid input up to ∼ 40 mm w.e. day −1 . These high fluxes, together with the fact that percolating water is able to bypass layers without pore space in our model, cause the relatively deep maximum in firn warming.

Firn aquifer
In accordance with Steger et al. (2017), we classify any firn with a vertically integrated liquid water content of more than 200 kg m −2 in April as an aquifer, irrespective of water saturation. This is necessary because our model is not able to simulate saturated conditions in the used configuration due to the neglect of impermeable layers. Introducing saturated conditions in SNOWPACK would require a definition of the pore space fraction available for liquid water storage. This quantity is rather uncertain and is assumed to be in the range of 40 % (Jansson et al., 2003) to 100 % . The above-mentioned threshold for firn aquifer delineation is based on a sensitivity estimation of the NASA Operation IceBridge accumulation radar to detect liquid water in firn (Miège et al., 2016).
The eastern part of the southern GrIS transect crosses the region where perennial firn aquifers were discovered by in situ observations in 2011  and mapped from radar measurements for 2010-2014 (Miège et al., 2016). The grey-shaded area in Fig. 11 indicates the horizontal extent of these mapped aquifers. The combination of RACMO2.3 and SNOWPACK underestimates the the upper limit of the firn aquifer's horizontal extent by approximately 50 m in elevation if one assumes only small changes in aquifer extent between 2010 and 2014. A brief sensitivity test of SNOWPACK with a lower fresh snow density, as described in Sect. 3.2, yields a firn aquifer that reaches higher elevations and thus reduces the mismatch (not shown). The reason for this improvement is that the lower near-surface firn density reduces the conductive heat loss of the aquifer to the atmosphere in winter. At lower elevations, SNOWPACK simulates a larger horizontal extent of the aquifer than inferred from radar data (Figs. 10 and 11). Apart from model uncertainties, this disagreement may be caused by different criteria used to delineate firn aquifers, where radar-derived mapping relies on the detection of a water table. The realisation of such a table may be prevented in this area by drainage of the aquifer into crevasses, where water either refreezes or enters the subglacial drainage system (Poinar et al., 2017).
Observations of the vertical extent of firn aquifers in this region return an average depth of 16.2 m for the water table and 27.7 m for the aquifer base (Montgomery et al., 2017). Comparing the depth of the modelled firn-aquifer top to observations is difficult, because observations derived from radar measurements return the depth of the water table and not the transition from dry to wet (but unsaturated) firn. The depth of the water table may roughly be in line with our simulations (Figs. 10 and 11) if an unsaturated wet layer between the aquifer top and the dry firn is assumed. The depth of the modelled firn aquifer base is clearly overestimated (> 40 m). This mismatch is likely related to a positive temperature bias of the deeper firn in our simulation (Steger et al., 2017). Observations indicate that refreezing conditions typically prevail at the base of aquifers (Montgomery et al., 2017)  the bottom of the model domain (Steger et al., 2017). Due to the present setting of SNOWPACK, which does not allow for saturated condition, the observed mean liquid water content of 16 % (Montgomery et al., 2017) is higher than modelled values. Figure 11 shows an expansion of the firn aquifer to higher elevations during the simulation period . This trend is in line with observations (2010)(2011)(2012)(2013)(2014)(2015)(2016), which indicate an inland expansion of aquifers in this area (Miège et al., 2016;Montgomery et al., 2017). Aquifer expansion is also apparent for other regions, where significant firn aquifer areas are modelled by SNOWPACK (Fig. 12). The highest fractions of firn aquifer area are simulated in the southeastern basins 4 and 5. In both basins, firn aquifers considerably expanded inland with time, particularly in Basin 4. This expansion to higher areas is partially compensated by a decrease in aquifer area at lower elevations, where porous firn is transformed to bare ice by increasing melt amounts (Fig. 8b). In basins 4 and 5, the mean surface elevation at which firn aquifers are modelled rises by ∼ 200 m during the simulation period . Upward migration of firn aquifers is also apparent in other basins, where basins 3 and 6 reveal smaller changes of 125 and 90 m, respectively, and Basin 8 a larger elevation increase of 215 m. The modelled aquifer extent of ∼ 59 000 km −2 for the entire GrIS, average over 2010-2014, is substantially larger than an estimate based on remote sensing data of 21 900 km −2 for the same period (Miège et al., 2016). The observational derived estimate provides a minimum extent of the aquifer due to inconsistent flight patterns of airborne data. Further potential causes for the deviation in estimates are discussed in Steger et al. (2017).
The formation of firn aquifers requires comparably high melt rates during the summer season and high annual accumulation rates (Kuipers Munneke et al., 2014). The dependence of firn aquifers on these parameters is also apparent in our simulation when modelled GrIS grid cells are plotted as a function of snowfall and liquid input (Fig. 13). The occurrence of firn aquifers is thereby restricted to a rather well separated space, which supports the hypothesis that snowfall and liquid input are the principal predictors for aquifer formation. The period of 1960-1979 has been selected for this analysis because it is identical to the spin-up period of our simulation. The relation is thus computed for steady-state climate without any long-term trends in the forcing. To assess the influence of a transient climate, firn aquifer occurrence as a function of snowfall and liquid input has also been computed for the period 2010-2014, which is identical to firn aquifer ob- servations by remote sensing (Miège et al., 2016). The zone of modelled aquifers shifts to a region with a higher ratio of liquid input to snowfall. This shift is likely caused by the changing climate conditions, where the spatial firn aquifer extent has not yet equilibrated to the new forcing.

Run-off partitioning
Run-off from the ice sheet can either originate from the melting of bare ice in the ablation zone, in which case run-off is assumed instantaneous, or from melting of snow/firn in the ablation or accumulation zone, in which case meltwater can be retained or refrozen. Partitioning run-off in these two classes yields insights in basin characteristics and in-dicates shifts in the accumulation and ablation area extent. Basins with high fractions of snow/firn run-off exhibit likely a higher uncertainty in run-off estimates due to potential storage of liquid water at the source location or along the routing path, the latter of which is not explicitly modelled. To distinguish run-off from ice and snow/firn, we apply a threshold for firn air content of 0.02 m. Basins 1 and 2 show relatively similar characteristics in terms of run-off partitioning (Fig. 14). In these dry northern regions with relatively wide ablation zones, run-off from snow/firn melt is at least an order of magnitude smaller than run-off from ice melt. Both basins reveal a strong positive trend in run-off from ice between 1990 and 2014 with an increase of 30.3 and 15.6 Gt a −1 over the 25 years. The eastern Basin 3 has a higher run-off fraction from snow/firn than the northern basins. Run-off from snow/firn is increasing in the later period (2.4 Gt a −1 ), although not on a statistically significant level. A very different picture emerges from Basin 4, which has a very narrow ablation zone. In this basin, approximately 87 % of run-off originates from snow/firn. Still, there is a considerable increase in run-off from ice during the second half of the simulation period (2.5 Gt a −1 ), which is caused by the decrease in pore space (Fig. 8b) and the gradual increase in the ablation zone. As a result, the snow/firn run-off fraction in this basin exhibits a significant negative trend (−9 %). Basin 5 reveals a similar pattern, but here run-off from ice and snow/firn is comparable in magnitude, particularly towards the end of the simulated period. This basin also reveals the highest interannual variability in snow/firn run-off fraction, which is related to the high interannual variance of winter (October-March) snowfall in this region (σ = 0.13 m w.e.). Variance in winter snowfall is also high in Basin 4 (σ = 0.14 m w.e.), but the sensitivity of the snow/firn run-off fraction on winter snowfall is lower due to a smaller ratio of ablation to accumulation area. The westerly basins 6-8 exhibit comparable runoff partitionings: all three basins are dominated by run-off from bare ice in the ablation zone and all these fluxes reveal a statistically significant positive trend in the second half of the simulated period. These trends are particularly strong in basins 6 and 8, where run-off from ice increases by 54.6 and 31.5 Gt a −1 over 1990-2014. In the more northerly basins, 7 and 8, there is a small but still significant increase in run-off from snow/firn. For the entire ice sheet, both run-off originating from ice and snow/firn increase at significant rates over the period 1990-2014 by 171.7 and 25.6 Gt a −1 , respectively. The run-off fraction from snow/firn decreased over this time by 6 %.

Conclusions
In this study, we analysed a SNOWPACK simulation carried out for the glaciated area of Greenland and for the period 1960-2014 with a focus on the liquid water balance (LWB) of the firn layer. The model was forced by output from the regional atmospheric climate model RACMO2.3 at the upper boundary. A comparison of the cumulative MB, derived with modelled SMB values and ice discharge data from observations, indicates an excellent agreement (R 2 > 0.99) with GRACE. The linear trend in cumulative MB improves when the SMB of Greenland's glaciated area is simulated by SNOWPACK instead of RACMO2.3 due to higher refreezing rates in SNOWPACK and thus reduced run-off from the ice sheet. However, the detrended mean seasonal cycles of these signals reveal significant discrepancies during the melt season. This mismatch can likely be attributed to neglecting run-off transit times and inaccuracies in the modelled tundra (snow) hydrology. The model also agrees well with observed changes in firn temperature along a 2700 km transect in north-western Greenland and with firn aquifer occurrence in the south-east. A direct comparison with temperature records from the western percolation zone of the ice sheet is not possible due to an overestimated bare ice zone extent in the model. Among other potential causes, such as climate biases in RACMO2.3, this mismatch is at least partly related to a bias in the fresh snow density parameterisation.
Temporally averaged LWB components over the simulation period  reveal that the balance is dominated by melt, run-off and refreezing in all basins. Modelled changes in retained liquid mass, evaporation and rainfall are typically at least 1 order of magnitude smaller, even for the more southerly basins. SNOWPACK simulates a mean refreezing fraction of 47 % averaged over the entire ice sheet. This quantity reveals a high spatial variability and is smallest for the northern GrIS (30 %) and largest in the south-east (75 %), where snowfall rates are highest. During the first half of the simulation period , there are no distinctive trends in the components of modelled LWB, but this changes for the second half , when surface melt fluxes significantly increase in all basins. These increases are reflected in run-off, particularly in the south-western area of the ice sheet where run-off increases by 0.31 m w.e. a −1 . Simulated trends in run-off generally exceed those in refreezing, which implies that the majority of the additional liquid water input runs off and thus contributes to sea level rise. The only exception is Basin 4 in the south-east, where most of the additional liquid input (∼ 70 %) is buffered in the firn. The simulated increase in refreezing, which is linked to the gradual expansion of the melt area in all basins, impacts firn properties by decreasing firn air content and increasing firn temperature. The exceptional melt in 2012 particularly causes a substantial warming of the firn, with a peak in the western percolation zone where modelled firn temperatures averaged over 2-10 m depth locally increase by more than 6 • C. SNOWPACK also simulates a migration of the firn aquifer area to higher elevations, which is, at least for an area in south-western Greenland, in line with the observations. Partitioning run-off according to its source (melting ice or snow/firn) shows that run-off from ice dominates on the ice sheet scale (78 %), with the highest run-off fractions (87 %) from snow/firn modelled in the south-east of the GrIS. Thus, this basin likely exhibits the highest uncertainty in runoff estimates due to possible retention of run-off in snow/firn at the place of origin or along the routing path.
The evaluation of our SNOWPACK simulation with various in situ and remote sensing observational data revealed several model inaccuracies, which are discussed in Steger et al. (2017). The current study emphasises the uncertainties in the applied fresh snow density parameterisation and the thermodynamic conditions beneath firn aquifers. The positive bias in the applied fresh snow density parameterisation for comparably warm climate conditions may be addressed by a comprehensive sensitivity test of SNOWPACK with different fresh snow density parameterisation for various climatic conditions. A particular focus should be placed on the disentanglement of processes influencing near-surface density (e.g. decrease in snow particle size during wind drift, vapour fluxes) and the statistical or physical representation of these processes in the parameterisation or the snow/firn model. The uncertainties in the simulated thermodynamic conditions beneath firn aquifers may be constrained with the increasing availability of in situ observations. Finally, our study reveals lateral routing of run-off as an additional relevant process that is not considered in our model. This is likely a less relevant issue for horizontal nearsurface redistribution of mass and energy, as surface melt typically reaches higher-elevated areas later in the season, which means that lower areas are already depleted of pore space and/or cold content and thus do not provide any more storage volume for upstream run-off. However, neglecting this process complicates comparisons of modelled SMB and GRACE on seasonal timescales. This shortcoming could be addressed by coupling SNOWPACK to an offline routing scheme for the GrIS. Nonetheless, the modelled MB and firn temperatures presented in this study compare favourably with remote sensing and in situ data and allow for a detailed evaluation of the GrIS's LWB between 1960 and 2014.
Data availability. The SNOWPACK data set presented in this study is available from the authors without conditions. Competing interests. The authors declare that they have no conflict of interest.
Special issue statement. This article is part of the special issue "Mass balance of the Greenland Ice Sheet". It does not belong to a conference.