Elevation change of the Greenland ice sheet due to surface mass balance and ﬁrn processes, 1960–2014

Observed changes in the surface elevation of the Greenland ice sheet are caused by ice dynamics, basal elevation change, basal melt, surface mass balance (SMB) variability, and by compaction of the overlying ﬁrn. The latter two contributions are quantiﬁed here using a ﬁrn model that includes compaction, meltwater percolation, and refreezing. The model is forced with surface mass ﬂuxes and temperature from a regional climate model for the period 1960–2014. The model results agree with observations of surface density, density proﬁles from 62 ﬁrn cores, and altimetric observations from regions where ice-dynamical surface height changes are likely small. In areas with strong surface melt, the ﬁrn model overestimates density. We ﬁnd that the ﬁrn layer in the high interior is generally thickening slowly (1–5 cm yr − 1 ). In the percolation and ablation areas, ﬁrn and SMB processes account for a surface elevation lowering of up to 20–50 cm yr − 1 . Most of this ﬁrn-induced marginal thinning is caused by an increase in melt since the mid-1990s, and partly compensated by an increase in the accumulation of fresh snow around most of the ice sheet. The total ﬁrn and ice volume change between 1980 and 2014 is estimated at − 3295 ± 1030 km 3 due to ﬁrn and SMB changes, corresponding to an ice-sheet average thinning of 1.96 ± 0.61 m. Most of this volume decrease occurred after 1995. The computed changes in surface elevation can be used to partition altimetrically observed volume change into surface mass balance and ice-dynamically related mass changes. percolation and refreezing on the vertical density proﬁles. It is impossible to isolate the initial density of snow that is truly fresh. We use two main sources for ρ 0 -data from the dry-snow zone: observations from North and Central Greenland from Benson (1962), and from the EGIG-line (Expédition Glaciologique Internationale au Groen-land, Morris and Wingham, 2011) in Central Greenland. We ﬁnd the following correlation ( r 2 = 0 . 52 ) between observed ρ 0 and annual mean surface T s (in

agreement would have been if the Ligtenberg (2011) values had been used? If there is no physical interpretation for these coefficients and if there is no reason for introducing them other than to improve the fit to the density profiles then this should be pointed out clearly in the paper. The fact that these coefficients are different in Greenland and Antarctica indicates that they are not representing a physical process whereby accumulation rate influences the rate of compaction directly. Instead, they must be correcting for some other process that is missing from the model, but happens to correlate with accumulation rate (albeit differently for each ice sheet). I think this should be pointed out in the paper.
In the revised manuscript, we acknowledge that the coefficients were altered to obtain a better fit between observed and modelled density profiles. We point out two possible resons for the difference in coefficients over Greenland and Antarctica, but we do not attempt to explain this difference in further detail.
Eqn 6. This equation has the wrong dimensions. It should be divided by the density of ice.
Thank you for pointing this out.

Eqn 6. This equation neglects the presence of liquid water within the firn. This might be a serious limitation under conditions for which a firn aquifer can develop. It would be better to present a more comprehensive treatment that includes liquid water, and then state what assumptions are being made.
We now start the section with the general formula for a mixture of water, ice, and air. We then point out what assumptions we make, and how the equation is simplified as the result of those assumptions.

P. 3550. Line 14. We set v_ice equal to the sum of all other components. Again, should this be the opposite sign from the sum of all other components?
Indeed, and we revised this sentence.
Sections 2.4 and 2.5. These are quite short to be separated as distinct 'methods' sections and could perhaps be combined with the respective discussion sections 3.1 and 3.2.
We acknowledge that sections 2.4 and 2.5 are short compared to 2.1-2.3. Nonetheless, we would like to retain these sections as they are, to keep the separation between models and techniques in section 2, and the results in section 3.
Section 4.4. The present study, with these sensitivity tests included, is just about OK for the time periods under consideration here. However, before doing more long runs using this model I would advise that the use of a reference accumulation rate is replaced by a better representation of the dynamical system representing grain-growth (e.g. Equations B1 and B2 of Appendix B of Arthern et al., 2010).
Within IMAU, we are working on including a dynamical description of the firn system by including grain growth and compaction as a result of overburden pressure rather than by a climatologically mean compaction rate. However, the overburden pressure formulation (in combination with prognostic grain growth) introduces a difficult-to-manage behaviour in the top of the firn column, when the overburden pressure becomes very small. Research on this issue is ongoing.

Figure 1. Needs more tick marks on x-axis. I think it is a log scale, but this is ambiguous.
This is now fixed.

General comments and questions:
I find the manuscript interesting, very well written and easy to follow. The results are convincing and presented in a clear way, and the analysis is thorough. The title states '1960-2013' but the presented results cover the period 1980-2013 because the years 1960-80 are used for the model spin-up. Would it not be more appropriate to change the title to 1980-2013? You are correct that we present results only since 1980. However, the underlying data set, which will be made available, starts in 1960, and the data set can be used to study temporal variability in the period 1960-1979, as long as one keeps in mind that we assume a steady state of the firn layer over that period. We would prefer to keep 1960 in the manuscript title, but if the editor decides otherwise, we will change the title.
You state in the abstract that the model results agree with firn core density data, and in the conclusion you state that you find a very good agreement. While I agree that the model produces convincing results, I think that Figure 3 shows that the model generally predict too high densities -especially considering that the model has been calibrated using the same core densities. I think that this should be pointed out more clearly -already in the abstract.
We have changed the abstract, such that the overestimation of density in the percolation area is emphasized.
I understand the need for the MO correction terms applied to the model results, but it seems that some term (s) are missing in the model since these MO corrections are different from what was found in Antarctica (Ligtenberg et al., 2011). I guess that this could mean that the MO factors could possibly change over time (?) and that this increases the uncertainty when running the model for longer time periods. I think that a more detailed discussion on this should be included in the manuscript.
This point was also raised by Reviewer 1. We have updated the discussion on the choice for MO-parameter values, and we are now explicit about the fact that this is merely done to optimize the fit between the model and observations. Data from 62 firn cores are used as validation, but only 57 cores are used in Figure 3, 59 in Figure 4. The number of cores used for determining the MO shown in Figure 1 is also not 62. The authors should make it clear why not all of the cores are used.
In figure 1, we used only cores with little surface melt, in order to tune the MO parameters for dry firn compaction only. We included one core with somewhat more surface melt (Das2), because it is the core with the highest accumulation rate, thus expanding the fit range. As a result, we used 22 cores, of which 7 extend to z_830. This is now mentioned in the figure caption.
In figure 3, we use 59 cores, since cores H2-1, H3-1, and H4-1 fall within the same grid box as the cores H1-1 and H5-1 shown in the figure. This is now mentioned in the figure caption.
In figure 4, we used the same 59 cores as in Figure 3. This is added to the figure caption.

Maybe I missed it, but have the authors provided information on the horizontal resolution of the firn model?
The horizontal resolution is mentioned in section 2.2 I think that a reference to Simonsen et al., 2013 would be appropriate in this manuscript, as this paper also describes the work of assessing a firn compaction model in Greenland.
We have added a reference to this work in the introduction, acknowledging previous work.
As also stated by the authors, the use of a mean value of the accumulation rate in Eq (4) and (5) has a significant impact on the compaction rate. I think that the authors describe clearly how and why this represents a limitation, but I think that some more discussion on why this is chosen anyways is needed. I reckon that there are valid explanations for not using the temporally varying accumulation rate, but it is not clear (to me at least) as it is now.
Within IMAU, we are working on including a dynamical description of the firn system by including grain growth and compaction as a result of overburden pressure rather than by a climatologically mean compaction rate. However, the overburden pressure formulation (in combination with prognostic grain growth) introduces a difficult-to-manage behaviour in the top of the firn column, when the overburden pressure becomes very small. Research on this issue is ongoing.

More specific comments and questions:
p. 3543, l. 1-3 : Bottom melt could in principle also be responsible for elevation changes. Or is this term included in what you call 'basal elevation change'?
We have included basal melt as a cause for elevation changes.
p. 3543, l. 3: 'and by the compaction of the overlying firn'. What do you mean with overlying here?
By "overlying" we mean the firn that is situated on top of the ice that makes up the Greenland Ice Sheet. While I am not a native speaker, I hope that the term "overlying" is sufficiently clear for the reader. These are the uncertainties of the observed elevation itself. We have changed "elevation change" to "elevation" in this sentence.
The choice for 15% is motivated in lines 1-18 of section 4.4. It is the best estimate of the uncertainties in the RACMO melt and accumulation fields. This is now explained in the figure caption, as indicated above. To us, the strength of using the true resolution of the model data is that the figure shows the vertical resolution of the model. We think that the comparison is needlessly complicated if we downsample the vertical density profiles to the core grid. For example, what to do with the depths at which only model data are availble?    Changed.

Introduction
The mass balance of the Greenland ice sheet has been negative over the last few decades (e.g. Zwally et al., 2011;Shepherd et al., 2012;Hurkmans et al., 2014). A common method to assess ice-sheet imbalance is altimetry, by which elevation changes are monitored by repeated scanning of the ice-sheet surface by active laser or radar instruments onboard airplanes or satellites. A crucial step in translating the observed Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | volume change to a mass change is to determine the density associated with the volume change.
The simplest assumption is that below the equilibrium line altitude (ELA), all mass change is caused either by ice-dynamical thinning or thickening, or by melting of ice. There, the ice density is used to convert from volume to mass change. Above the ELA, it is often assumed that ice-dynamical elevation changes are negligible, and that volume changes are caused solely by changes in the firn layer. Here, a fixed density is adopted to convert from volume to mass, for example a fresh-snow density, or a representative density for the entire firn layer (e.g. Davis et al., 2005;Wingham et al., 2006 for Antarctica, andThomas et al., 2006 for Greenland). This assumption is still commonly made over smaller ice caps and glaciers (e.g. Gardelle et al., 2012;Moholdt et al., 2012;Gardner et al., 2013), where the error associated with this approach is generally assumed to be small. For the Greenland ice sheet, such an approach can compromise the accuracy of the retrieved mass changes, as the current mass loss is divided between surface mass balance (SMB) and ice-dynamical changes. Moreover, the density of the Greenland firn layer is susceptible to change, invalidating the choice for a constant density. In fact, elevations change in the ice-sheet interior is primarily attributed to variability in snow accumulation (McConnell et al., 2000a). More recent altimetry-based mass balance estimates of the Greenland ice sheet (e.g. Sørensen et al., 2011;Khan et al., 2014;Hurkmans et al., 2014) build upon an empirical firn model that takes into account firn compaction, and the formation of ice lenses due to melt and refreezing in annual layers (Reeh et al., 2005;Reeh, 2008). Another firn model for altimetry corrections was developed by Zwally and Li (2002) and Li and Zwally (2011); this model is driven by satellite-observed surface temperature, and a fixed, firn-core derived relation between temperature and accumulation change. Simonsen et al. (2013) developed a firn model similar to ours, and compared it to Ku-band radar observations of annual layering of the Greenland firn.
Here, we present a time series of elevation changes over the Greenland ice sheet due to changes in depth and mass of the firn layer, i.e. variability and change in the surface mass balance and associated firn processes like compaction, percolation, and refreezing. These time series extend from 1960 up to and including 2014, at a horizontal resolution of 11 km× 11 km. The results are obtained using a semi-empirical firn compaction model (Ligtenberg et al., 2011) that is forced with surface mass fluxes and temperature from the polar-adapted regional climate model RACMO2.3 (Noël et al., 2015). Combining the time series of firn depth and mass from this firn model allows one to convert satellite altimetry observations of volume change to mass change, and to partition this mass change into an ice-dynamical and a firn/SMB component. As an added advantage, the firn model explicitly calculates all firn and SMB processes that modify the surface elevation, allowing the analysis of the individual components of the surface elevation change in greater detail.

Model, data, and methods
Elevation change of the Greenland ice sheet surface is simulated with a model of the firn (Sect. 2.1), which is forced at its upper boundary by temperature and surface mass fluxes from the regional climate model RACMO, version 2.3 (see Sect. 2.2, and Noël et al., 2015). In Sect. 3, we evaluate the firn model results against 62 shallow and medium-depth firn cores (details in Sect. 2.4), and against airborne laser altimetry collected in areas of limited ice-dynamical activity (Sect. 2.5).

The firn model
The firn model (IMAU-FDM v1.0, Ligtenberg et al., 2011) describes the temporal evolution of firn compaction, meltwater percolation and refreezing, and temperature in a vertical, 1-D column of firn and ice. The top of the uppermost model layer represents the surface of the ice sheet (h), which moves up and down in time. The vertical velocity of the ice-sheet surface due to firn and SMB processesḣ f is given by: where we neglect vertical displacement of the surface by horizontal velocity divergence or horizontal advection of mass in the snow and firn column. The velocity components v rep-4 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | resent solid precipitation (v acc ), surface sublimation (included in v acc ), snowdrift sublimation (v snd ), snowdrift erosion (v er ), snowmelt (v me ), and firn compaction (v fc ). The solid precipitation v acc represents solid precipitation minus surface sublimation, and thus, it is not the accumulation rate as it is usually understood. Snowdrift sublimation v snd differs from surface sublimation in the sense that it represents sublimation of drifting snow that is whirled up from the surface by surface winds (Lenaerts et al., 2010). Snowdrift erosion represents the horizontal redistribution of surface snow by surface winds, and can be positive (deposition) and negative (erosion).
In a steady-state firn layer, the long-term average vertical mass flux through the lower boundary of the firn column equals the mass flux through the upper boundary. This is represented in the firn model as a constant vertical velocity v ice that equals the mean SMB (v acc + v snd + v er + v me ) but of opposite sign, over a reference period for which steady-state is assumed. The choice of this reference period is discussed in Sect. 2.3. This model setup is similar to earlier process-based models like Zwally and Li (2002).
An important upper boundary condition is the density of fresh snow, ρ 0 . While valuable observations of density in the percolation area of the Greenland ice sheet are available (Braithwaite et al., 1994;Brown et al., 2012), we choose to construct a parameterization of ρ 0 based on observations from the dry-snow zone only. The observations from Braithwaite et al. (1994) are averages of density over the uppermost 1 m of the snow/firn layer, which includes the effect of percolation and refreezing on the vertical density profiles. It is impossible to isolate the initial density of snow that is truly fresh. We use two main sources for ρ 0 -data from the dry-snow zone: observations from North and Central Greenland from Benson (1962), and from the EGIG-line (Expédition Glaciologique Internationale au Groenland, Morris and Wingham, 2011) in Central Greenland. We find the following correlation (r 2 = 0.52) between observed ρ 0 and annual mean surface temperature T s (in degrees Cel-Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | sius) simulated by RACMO2.3 (see Sect. 2.2): (2) After deposition, the fresh snow starts to become denser. The effect of dry-firn compaction is represented by v fc in Eq. (1). Compaction is an increase of firn density ρ in time (t), expressed by Eq. (4) in Arthern et al. (2010) as: Here, C, E c , and E g are constants,ḃ is the mean annual accumulation over a reference period (Sect. 2.3), g is the gravitational acceleration, ρ i the ice density of 917 kg m −3 , and R the gas constant. T is the firn temperature that varies with depth. The rate constant C has a value of 0.07 for ρ ≤ 550 kg m −3 , and a value of 0.03 above that density. This captures the higher densification rate near the surface due to sliding of snow grains relative to each other at low densities (Arthern et al., 2010).
To evaluate Eq.
(3), we compare modelled depths of the 550 and 830 kg m −3 density layers (z 550 and z 830 , respectively) to observations at 62 firn core locations around the Greenland ice sheet (see Sect. 2.4). As we find a systematic departure of the modelled values, we introduce a correction term MO, defined as the ratio of modelled and observed values of z 550 and z * 830 (Ligtenberg et al., 2011), where z * 830 = z 830 − z 550 . Figure 1 shows MO as a function ofḃ. As the MO-values need to represent dry compaction, we selected 22 firn cores with little surface melt. Linear least-squares fitting then yields the following MO-relations for Greenland: The coefficients in these two equations are different from a previous application in Antarctica Ligtenberg et al. (2011). We use these updated coefficients to improve the fit with 6 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | observed density profiles. There is no physical interpretation for these coefficients to be different. The different set of coefficients for Greenland and Antarctica could point to a process not presently captured in the model (that happens to correlate with accumulation rate). Alternatively, the coefficients could be different because the range of accumulation rates on which the Antarctic parameterization (Ligtenberg et al., 2011) is based, is biased to lower accumulation rates than those found in Greenland.
Equation (3) is multiplied with the correction factors MO in Eqs. (4) and (5), that are not allowed to be smaller than 0.25 (Ligtenberg et al., 2011). The accurate performance of these densification expressions is further demonstrated by fully independent comparisons against in-situ (Larson et al., 2015) and remotely-sensed  densification rates. However, the use of a mean value ofḃ in these equations has an important limitation: in reality, compaction is determined by the overburden pressure of the overlying snow. By using a constant meanḃ rather than the instantaneous overburden pressure, the firn compaction variability is dampened significantly. Following an accumulation event, the model only takes into account the increase in compactable material, not the effect of increased overburden pressure on the underlying firn.
Rain is added to the surface snowmelt flux. This liquid water is allowed to percolate into the firn. Each layer has a maximum irreducible water content W c , depending on the density (Coléou and Lesaffre, 1998). Meltwater will refreeze as soon as it encounters a layer that can accommodate both the space of the refreezing water and the latent heat that is released upon refreezing. For details see Ligtenberg et al. (2011) andKuipers Munneke et al. (2015).
Suppose a firn column with total depth z i and an observed or modelled density profile ρ(z), consisting of a mixture of air (with density ρ a ), water (with density ρ w ), and ice (density ρ i ). Using W c from Coléou and Lesaffre (1998) (which can vary with depth z), it can be 7 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | shown that the firn air content F (in m) is given as Here, we will approximate this general expression by assuming that ρ a ρ i , and that W c (z) = 0, yielding As an extreme example, a uniform value of W c = 0.10 (i.e.,10% of pore space filled with water throughout the firn column) yields a 10% overestimation of F .
The firn model is also used to simulate surface elevation changesḣ f in the ablation area. Here, the prescribed RACMO2.3 mass fluxes determine the ice ablation rate assuming the ice density ρ i . Section 2.3 describes the method to derive the SMB-induced surface elevation changeḣ f in the ablation area.

Model forcing from RACMO2.3
At the upper boundary of the 1-D firn column, the firn model is forced with surface temperature and mass fluxes from RACMO2.3 (Noël et al., 2015), a regional climate model that is adapted to simulate climatic conditions over ice sheets. The horizontal spatial resolution of RACMO2.3 is 11 km. Forcing data are available for the period 1 January 1960-31 December 2014, in time steps of 6 h.
RACMO2.3 supersedes RACMO2.1 van Angelen et al., 2012). In the new version, the cloud microphysics, surface and boundary layer turbulence, and radiation transport have been updated (Van Wessem et al., 2014). The most pronounced effect of these updates on the SMB is an increase in summer snowfall events, decreasing the amount of snow and ice melt in the percolation and ablation area (Noël et al., 2015). The 8 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | agreement between RACMO2.3 SMB and mass balance stakes in these areas is improved. The ELA is lower and in better agreement with observations. This is expected to improve the description of firn processes in the percolation area of the ice sheet.

Modelling strategy
While RACMO2.3 itself contains a multi-layer snowpack with the same compaction and meltwater routines as the firn model, the rationale for using an off-line firn model is the ability to spin up the firn layer with a reference climate, until it is in equilibrium with that reference climate. This circumvents the difficult task of assuming an initial condition of the firn layer at the start of the RACMO2.3 simulation, that is sufficiently accurate for correctly determiningḣ f . Furthermore, because of its computational demands, RACMO2.3 cannot be used for sensitivity tests, in contrast to an offline firn model. Finally, the vertical resolution is higher in the offline model.
Still, the spin-up procedure requires that we define a reference climate, i.e. a period of time in which neither the properties of the firn nor the reference climate forcing exhibit significant trends. Recent modelling and observations reveal that the Greenland ice sheet SMB has decreased since the beginning of the 1990s (e.g. Shepherd et al., 2012;Enderlyn et al., 2014). As a result, thinning has increased sharply since the mid-1990s (Csatho et al., 2014) along the margins of the ice sheet. Clearly, the reference period should not include this period, and should end preferably some years before its onset. Therefore, our modelling strategy is that we choose the first 20 years of RACMO2.3 forcing data (01 January 1960-31 December 1979) and spin up the firn column at each location with a loop over this 20 year period, until the properties of the firn layer have converged to an equilibrium. By doing so, we assume that the pre-1960 climate (i.e. the reference climate) can be represented by a sequence of 20 year periods. In practice, equilibrium is reached when all the mass in the firn layer is refreshed once since the start of the spin-up. The duration of the spin-up is therefore computed as the thickness of the firn layer (from the surface to the depth at which the ice density is reached) divided by the mean annual accumulation rate. A major 9 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | uncertainty in the calculated firn depth changes in this study stems from this assumption of reference climate. We will quantify this uncertainty in Sect. 4.4.
As a second important assumption, we set the surface elevation changeḣ f over the reference period to zero. After all, the modelled firn layer at the end of the reference period (31 December 1979) is the result of the spinup procedure that uses multiple loops over the reference period to reach an equilibrium state, plus 20 years of model integration using the same data as in the spinup procedure. For the accumulation area, it means that the amount of mass leaving the bottom of the firn layer (with a velocity v ice ) is assumed equal to the total mass added to and retained in the firn column by snowfall and refreezing, minus runoff.
For the ablation area, the assumption ofḣ f = 0 during the reference period implies that the downward velocity from a negative SMB (ablation) is balanced by an upward and equal flux of emergent ice. The emergent ice flux is represented by the term v ice . In Eq. (1), we set v ice equal to the opposite sign of the sum of all other velocity components for the reference period. Thereafter, v ice retains the same value but the other parameters are free to evolve.
In this framework,ḣ f (like presented in Fig. 7) in the ablation area represents the surface elevation change due to the anomaly of surface melt with respect to the reference period. To clarify, the change in surface elevation itself does not have to be zero over the reference period: it can change due to ice-dynamical thinning or thickening. Only the ablation-driven surface elevation changeḣ f is assumed to be zero.
Note that our choice of the period 1960-1979 for a representative reference climate implies that modelled surface elevation after 1980 is allowed to evolve freely due to SMB and firn processes. It is not bounded by the requirement thatḣ f should be zero over the entire simulation period, like in studies addressing the Antarctic ice sheet (Ligtenberg et al., 2011;Pritchard et al., 2012). As the firn layer starts to evolve freely from 1980 onwards, we present time series starting in 1980, although the complete time series of elevation change start in 1960.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Firn cores
To evaluate the firn model, we collected vertical profiles of firn density from 62 shallow cores from widely varying locations across the Greenland Ice Sheet (see the map in Fig. 2), drilled between 1995 and 2012. Among the cores are those drilled for PARCA (Program for Arctic Regional Climate Assessment, McConnell et al., 2000b;Mosley-Thompson et al., 2001;Hanna et al., 2006;Banta and McConnell, 2007); cores from the Arctic Circle Traverses (ACT, Box et al., 2013); cores from the lower part of the EGIG line ; and Das 1 and Das 2 (e.g. Hanna et al., 2006).
Vertical profiles of density from firn cores are usually based on the mean density of 0.5-2 m long sections. Some researchers log the midpoint of each section as the depth of the section, others use the top or the bottom of the section. Here, all 62 profiles have been interpolated to give the mean density at the midpoints of each core section. For each core, the collection date is known, and vertical profiles of modelled firn density are extracted from the model at the time closest to the collection date, and from the nearest model grid cell.

Laser altimetry
Since 1992, NASA's Airborne Topographic Mapper (ATM) has carried out laser surveys of the Greenland ice sheet surface. If sufficient repeat observations are available, a time series of observed surface elevation change can be constructed, spanning the period 1992-2013.
To do so, we use the Level 2 "Qfit" product which provides the waveform-fitted elevations for the centroid of each laser return. In order to derive elevation changes, we interpolated the Qfit point cloud for each campaign to a reference grid with 30 m spacing. We then selected reference grid points with at least one observation in each of five epochs: 1993-1996, 1997-2000, 2001-2005, 2006-2009, and 2010-2013. Of these points, we selected 15 within the central and northern interior of the ice sheet, where a relatively small contribution of ice dynamics is expected. Based on differences in elevation obtained from crossovers within a few weeks, we estimate 1σ errors in the elevation observationsto be ±10 cm, which Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | includes interpolation error. Results of the comparison between the firn model and the ATM data are found in Table 1 and Figure 5.

Vertical profiles of density
We use vertical profiles of density from 62 firn cores to assess the performance of the firn model. This evaluation is not independent, as we used the depths of the 550 and 830 kg m −3 density levels from these cores to tune the densification parameterization in Eqs. (4) and (5). Still, we can compare the shape of the profiles beyond these two levels, and we can assess the impact of melt and refreezing on the vertical density profile. Figure 3 shows the observed and modelled density profiles for all core locations. Each panel includes the mean accumulation and melt from RACMO2.3 (in mm yr −1 ) for 1960-2014, and the ratio R ma of these melt and accumulation averages. The vertical resolution of the firn core data is typically 0.5-2.0 m, thereby smoothing out the effect of ice lenses and higher-density layers. The model data in Fig. 3 is shown at full resolution, i.e. with layers of 5-10 cm thickness. The high-density layers usually represent thick layers of refrozen meltwater with a density close to that of solid ice.
Up to an R ma value of 0.3-0.4, the agreement between the firn model and the observations is good. But for higher R ma , the firn model starts to overestimate the density throughout the firn column. Figure 4 compares the observed and modelled firn air content F , showing R ma in colour. The model bias clearly increases for higher R ma . It means that there are three possible causes for the misfit, that are not mutually exclusive: (1) RACMO2.3 simulates too much melt in the percolation areas, causing the firn to fill up quickly with too much refrozen meltwater; (2) RACMO2.3 simulates too little accumulation, providing insufficient pore space to store meltwater; and (3) the firn model should allow for more and more rapid downward percolation of meltwater without letting it refreeze. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Regarding a possible overestimation of melt in RACMO2.3, there is limited evidence that the amount of melt observed by an automatic weather station at location S10 (67 • 00 N 47 • 01 W, 1850 m a.s.l.) is indeed about 20 % smaller than simulated by RACMO2.3 (Noël et al., 2015). Further north, Harper et al. (2012) find the equilibrium line around the EGIG line at about 1100-1200 m a.s.l., while the equilibrium line altitude in RACMO2.3 is at ∼ 1650 m a.s.l., about 45 km further inland. For firn cores H1-1 down to H5-1 (Fig. 2), it is clear that under RACMO2.3 forcing (with melt larger than accumulation) a firn layer cannot be sustained, whereas in reality there is a shallow firnpack with infiltration ice layers. There is very limited reliable information about melt fluxes from other parts of the percolation area around the ice sheet, so we cannot conclude whether the overestimation of modelled melt flux is structural.
The other possible source for the misfit is the percolation scheme in the firn model itself. The firn model adopts a so-called "tipping bucket"-approach, where meltwater is allowed to move downward from one discrete layer to the next whenever the first layer is saturated. In practice, the percolation is more complex, and vertical meltwater transport through confined channels (pipes) is known to occur (Marsh and Woo, 1984;Humphrey et al., 2012). Piping of meltwater is a way to evacuate more meltwater towards the bottom of the firn layer, reducing the amount of refreezing in the firn itself. Alternatively, intermediate thick ice layers may serve as an impermeable surface along which the water can run off. Both processes increase the runoff and decrease refreezing and density. At present, we cannot assess the performance of the firn model in more detail, since we cannot easily isolate it from errors in the model forcing from RACMO2.3.
It is unclear what exactly the model bias in the percolation zone implies for the modelled rates of surface elevation change. We speculate that, if too much refreezing is the cause for the density overestimation, then a prescribed increase in surface melt would underestimate the rate of surface lowering.
In the high interior of the Greenland ice sheet, horizontal surface ice velocities are low (generally less than 10 m yr −1 , Joughin et al., 2010), and elevation changes resulting from ice-dynamical effects are expected to be small. Figure 5 shows time series of observed surface elevation change from the ATM lidar, along with the surface elevation change predicted by the firn model. Surface elevation change rates at the 15 test sites range from −6.6 to 5.1 cm yr −1 over the altimetry record (map in Fig. 5, Table 1). The sites in the central east (site 1, 2 and 3) had the highest rates of surface rise, with rates increasing inland. Sites 8, 9 and 10 near the northwest margin uniquely show decreasing elevations. The time series of observed surface elevation change (panels in Fig. 5) show the substantial variability between nearby locations in both time and space.
The firn model provides the change in surface elevation due to only variations in snow accumulation and firn density, assuming constant vertical ice motion. Therefore, the difference between the observed change and the modeled elevation change represents the elevation change due to vertical ice motion (ice dynamics) and error. We assume that in the ice sheet interior, variations in ice dynamics occur over timescales that are long (centuries) relative to the observational record, and can therefore be approximated by a linear trend. Under this assumption, the residual between the observed and modeled surface elevations will decrease or increase at a rate equal to the difference between the reference and actual submergence rates. The trend in residuals is therefore the anomaly in the submergence rate from the reference, which is assumed to approximate steady state, and provides an estimates of the contribution of ice-dynamical change to surface elevation. These trends in residuals are given in Table 1. In most cases, these trends are not statistically significant, indicating a submergence velocity close to the reference state. At site 3, the trend in residuals is nearly 7 cm yr −1 , which accounts for more than the 6 cm yr −1 of observed increases, indicating thickening. At site 6, a negative trend in residuals of 1 cm yr −1 opposes the 2 cm yr −1 of observed surface raising, suggesting opposing contributions from dynamics Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | and accumulation. At sites 8 and 10, the strongly negative trend in residuals is larger than the observed surface lowering, indicating that increased accumulation is partially offsetting relatively rapid dynamic thinning.
If the trend in residuals between the observed and modeled surface elevations provides the linear contribution in ice dynamics plus the error, the error is then assessed by the rootmean-square (RMS) of the detrended residuals (Table 1). This is equivalent to adjusting the firn model time series to the ice dynamic trend (shown as green curves in the panels of There is a remarkable contrast between the firn in the NW and the NE, with the NW having higher F . This can be explained by more snowfall in the NW, and higher sublimation in the NE due to a lower relative humidity. By adding all the velocity components in Eq. (1), we findḣ f , the firn thickness change per unit of time due to all firn and SMB processes. We can accumulate the thickness changes over longer periods to get multi-year maps over the ice sheet. Figure 7 showsḣ f (in cm yr −1 ) for the periods 1980-2014, 1980-1995, and 1995-2014, respectively. This surface elevation change is with respect to the reference period 1960-1979, during which zero surface elevation change (due to firn and SMB processes) is assumed. Again, there is a pronounced pattern of modest thickening of the firn layer in the interior (most notably towards the east), and moderate to strong thinning of the firn layer around the margins. The interior thickening is of the order of 1-5 cm yr −1 , or 34-170 cm over the entire 34 year period. The marginal thinning rates are much larger, and can be up to 20-50 cm yr −1 or 6-18 m over the entire period, with the highest values in the southeast. In contrast to the Summit dome firn layer, that of the southern dome of the Greenland Ice sheet is thinning.
It is clear that the patterns have changed over this period. The surface elevation change map over 1980-1995 (Fig. 7b) shows thinning along the southeast, south, west, and northwest coasts. Thickening is occurring in the interior (mainly east of the divide), and along the north and northeast coasts. Since 1995, thinning has intensified, and spread over the entire coastal margin. The thickening moved to the west of the interior. The aggregate picture for the period 1980-2014then shows moderate thickening up to 5 cm yr −1 in the eastern and northern interior. Thinning occurs all around the margins, with the smallest rates (0-15 cm yr −1 ) in the northern and eastern coastal regions. The largest rates (exceeding 40 cm yr −1 ) occur in the southeast, and in the western ablation area.

Decomposing the trends
The firn model allows for a decomposition of theḣ f -signal into its velocity components (Eq. 1). The upper panels in Fig. 8 show this decomposition for the period 1980-2014. The thickening in the eastern interior (Fig. 7a) can be almost fully ascribed to a positive accumulation anomaly (Fig. 8a), offset by a small increase in firn compaction due to this extra Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | firn (Fig. 8g). In the lower accumulation area, the positive accumulation anomaly is offset by a significant increase in surface melt, giving zero or slightly negativeḣ f in the western percolation area. In the south, increased surface melt dominates the thinning signal, whereas accumulation, firn compaction, and snowdrift anomalies play a minor role. In the southeast, melt has increased and accumulation decreased significantly. As the absolute values of both accumulation and melt are large in this region, we find here the largest rates of firn-driven surface lowering seen in Greenland.
For the period 1980-1995, the accumulation anomaly differed from the 1995-2014period, as shown in panels b and c of Fig. 8. A negative accumulation anomaly (partly offset by a positive firn compaction anomaly, panel h) explains the firn-driven surface lowering in the northwest. In the absence of significant melt anomalies (panel e), the thickening in the eastern and northeastern interior can be almost fully ascribed to a positive accumulation anomaly (panel b). This is mirrored in a small negative firn compaction anomaly (panel h).
Over almost the entire ice sheet, with the exception of the southeast, the period 1995-2014shows a positive accumulation anomaly (panel c). At lower elevations however, the firn thickness change is dominated by the strong melt anomaly over this period (panel f).
The velocity components that always lead to a decrease ofḣ f , melt and firn compaction, are negative by definition. To complement, we can add up the snowdrift and surface sublimation velocities whenever they lead to a surface lowering. The partitioning of the surface lowering into these components of negative velocities is shown in Fig. 9. The lowering is dominated by firn compaction (panel b) in the interior, and more and more by melt around the margins. There, the firn layer is thinner, which reduces the compaction potential. In the dry northeast, there is a relatively large contribution from sublimation (up to 30 %, panel d).
This is caused by a combination of strong winds and a relatively low humidity, promoting snowdrift sublimation. For another part, the relative contribution increases as the firn compaction is small due to lower firn temperatures, and due to the relatively small thickness of the firn layer.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Error estimate
An important source of uncertainty inḣ f is the steady-state assumed for the spin-up of the firn model. As explained in Sect. 2.3, the present model setup assumes that the climate under which the firn was formed can be represented by a loop over the forcing data from 1960 to 1979. A reconstruction based on firn cores and a previous version of RACMO2 found large interdecadal accumulation variability over the past 400 years, and an accumulation increase by 12 % over the period 1600. The mean reconstructed, ice-sheet wide accumulation over this period is 782 Gt yr −1 . For 1960-1979, it is 786 Gt yr −1 , i.e. 0.5 % larger than the long-term average. Other reconstructions show the 1960-1979 solid precipitation flux to be 0-10 % different from the period between ∼ 1850 and the present day (Wake et al., 2009;Hanna et al., 2011). For the error analysis, we therefore assume that the 1960-1979 accumulation fluxes (precipitation minus evaporation/sublimation) can differ by up to 15 % from the long-term accumulation history. We denote this accumulation uncertainty by σ˙b (in mm yr −1 ).
For snowmelt, we also assume a maximum deviation of 15 % of the 1960-1979 melt flux compared to the long-term history. There is limited evidence for this, but according to reconstructions from Hanna et al. (2011), the runoff flux for 1960-1979 differs about 5-10 % from the 1870-2010 mean. This uncertainty is given as σṁ (mm yr −1 ).
For the suite of 62 firn core locations, we perform 4 sensitivity tests, in which we increase and decrease melt or accumulation by 15 % after completion of the spin-up, respectively.
The resulting drift inḣ f for 1960-2014gives an error in the firn thickness change. The error in surface elevation change due to σ˙b is written as σḣ ,ḃ . Analogously, the error in surface elevation change due to σṁ is given as σḣ ,ṁ . By relating σḣ ,ḃ and σḣ ,ṁ empirically to accumulation and melt at the core locations, we can expand the error product over the entire ice Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | sheet. These relations are: σḣ ,ḃ = σ˙b(0.107 + 3.609 × 10 4ḃ ) (8) σḣ ,ṁ = σṁ(0.225 + 1.064 × 10 3ṁ ), withḃ andṁ in mm yr −1 .
The uncertainties σ˙b and σṁ cannot be regarded as independent. The SMB module in RACMO2.3 contains interactions between accumulation and melt. We identify the meltalbedo feedback as the most important interaction. As an example, a negative bias in summer snowfall could lead to an excess of summer melt because albedo will be underestimated. To capture this dependence in the error analysis, we assume the errors in surface elevation change due to uncertainties in the melt and accumulation fluxes as dependent, and add them up linearly (σḣ = σḣ ,ḃ + σḣ ,ṁ ).
The resulting total error and its components are shown in Fig. 10. In the interior, the total error (panel c) is understandably dominated by the accumulation uncertainty (panel a). Towards the ice-sheet edge, the melt uncertainty (panel b) starts to dominate the total uncertainty. The total uncertainty increases from 0.2-1.0 cm yr −1 in the interior, to 10-20 cm yr −1 in the lower percolation area. The largest total error (more than 40 cm yr −1 ) is found in the southeast, where high snowfall rates coincide with large amounts of melt.
For low and medium values ofḃ, the term between brackets in Eq. (8) is smaller than 1. It means that an uncertainty in the accumulation rate leads to a smaller uncertainty in the elevation change. The explanation for this behaviour is simple: the densification rate in Eq.
(3) linearly depends on the mean accumulation rate, so that an elevation increase due to more snowfall is partially offset by a more rapid densification. Forḃ > 2427 mm yr −1 , the term between brackets in Eq. (8) becomes larger than unity. There is no physical explanation for this behaviour: it is caused by the empirical nature of Eq. (8). However, we decided not to cap the uncertainty, for the following reason: the densification rate in Eq. (3) depends on a 20 year mean accumulation rate, whereas true densification at a particular depth in the firn depends on the immediate overburden pressure from overlying firn, which can be more variable than the long-term mean. This was already shown to dampen the modelled Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | variability in densification rate compared to observations . Therefore, the observed, short-term firn thickness change variability could be of the same order of magnitude as the accumulation rate variability for large accumulation rates. Figure 11 shows the cumulative volume change of the Greenland ice sheet as a consequence of changes in firn and SMB processes. Until the late 1990s, the total volume change was small. Since 2000, the total volume has decreased by about 3295 ± 1030 km 3 due to firn and SMB processes alone. Averaged over the ice sheet, this is a mean surface lowering of 1.96 ± 0.61 m. Almost all of this total volume loss took place in the part of the Greenland ice sheet where the surface is under 2000 m a.s.l. Fig. 11 show the partitioning of the volume change for the entire ice sheet and the part elevated above 2000 m a.s.l., respectively. Over 1980-2014, the volume loss through melt was slightly over 4700 km 3 . This loss was partly compensated for by an increase in snow accumulation of about 1500 km 3 . On most of the ice sheet, firn compaction has accelerated (become more negative) due to an increase in accumulation.

Panels b and c of
Above 2000 m, the effect is clearly visible (Fig. 11c). Integrated over the ice sheet however, we see a small slowdown in firn compaction, corresponding to about 500 km 3 (Fig. 11b).
The firn compaction anomaly is dominated by the southeastern part of the ice sheet, where snowfall has decreased strongly in an absolute sense (less so in a relative sense) and firn compaction has slowed down ( Fig. 8h and i).
Up to about 2005, firn volume change above 2000 m a.s.l. has been dominated by accumulation variability (consistent with e.g. McConnell et al., 2000a), and below 2000 m a.s.l., the volume change was mainly melt-driven. This is consistent with the original speculation in early reports of what would happen to the Greenland ice sheet in response to global warming. After 2005however, the total firn volume above 2000 m a.s.l. has started to decrease, mainly because surface melt has migrated inland (e.g. Fettweis et al., 2011), but also because the accumulation increase, clearly visible between 1980 and 2000, stagnated 20 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | in the 2000s. It remains to be seen if the paradigm of interior thickening under atmospheric warming can stand up against the inland migration of the area of surface melt.
The extreme melt season of 2012 (Nghiem et al., 2012;Tedesco et al., 2013) is clearly visible in the results of the firn model. A large drop in total volume of 1386 km 3 is seen in the summer of 2012, of which 1150 km 3 is contributed to melt. Melt in the part of the ice sheet above 2000 m contributed almost one third (371 km 3 ) to this volume anomaly. For perspective, the volume loss above 2000 m a.s.l. in the summer of 2012 is equal to the volume gained by snowfall in the interior in the 16 years between 1980 and 1996.

Altimetry correction and mean density of firn-related mass loss
The present dataset of firn thickness and mass change is primarily aimed at the correction of altimetry products, allowing for the extraction of an ice-dynamical thinning/thickening signal. The procedure for doing so is as follows: suppose that an altimetry sensor measures a surface elevation changeḣ(t 0 , t 1 ) ≡ (h(t 1 ) − h(t 0 ))/(t 1 − t 0 ) between initial time t 0 and time t 1 > t 0 . The firn model computes a surface elevation change due to firn and SMB processesḣ f (t 0 , t 1 ) ≡ h f (t 1 ) − h f (t 0 ). The supposed ice-dynamical contribution (neglecting vertical bed movement, but this could be included; Sørensen et al., 2011) is theṅ h d =ḣ −ḣ f . The associated mass changeṁ d is simply computed aṡ The mass changeṁ d is caused by horizontal convergence or divergence of ice flux, arising from, e.g. ice-flow acceleration propagating from the margin, long-term changes in the icesheet viscosity (Colgan et al., 2015), and transient variations in ice flow due to long-term changes in accumulation. The mass change due to the SMB,ṁ f , is computed directly from the RACMO2.3 forcing, using SMB anomalies with respect to the appropriate reference period (1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979). This is by far the simplest approach and completely consistent with the firn model, that uses the same forcing and the same reference period. The total mass change at the given locationṁ is then computed aṡ m =ṁ d +ṁ f .

Conclusions
In this study, we used a time-dependent, semi-empirical model for firn compaction, meltwater percolation, and refreezing. We forced the model with surface mass fluxes and temperature from the regional climate model RACMO2.3 for the period 1960-2014. By forcing the model with all mass fluxes, including melt, the result is a data series of surface elevation change over the entire ice sheet, not only due to firn processes, but also due to anomalies in the SMB. We defined a reference period from 1960 up to and including 1979, in which we assumed the surface elevation change to be zero due to firn and SMB processes. In the ablation zone, the computed surface elevation change represents the ablation anomaly with respect to the 1960-1979 mean ablation.
The firn model was calibrated against vertical profiles of firn density from more than 60 shallow and deep firn cores collected around Greenland in the past two decades. This ensured a very good agreement between observed and modelled vertical density profiles, especially in regions where the annual surface melt flux is small (less than about 20 %) compared to the mean annual accumulation. In regions with higher melt, the firn model overestimates the density below the surface. Potentially, this underestimates the presented rates of surface lowering in the percolation area.
The computed surface elevation change was compared against lidar observations of surface elevation change at locations where we expect the ice-dynamical changes to be small or gradual in time. At most locations, we find a good fit of the modelled elevation change rates to the observed ones.
Between 1980 and 2014, we see a pronounced pattern of small thickening of the firn layer in the high interior, of 1 to 5 cm yr −1 , caused predominantly by an accumulation increase over this period. Around the margins of the ice sheet, in the percolation and ablation areas, the surface is lowering, at rates of up to 20-50 cm yr −1 . This is mostly caused by an increase in surface melt, augmented in the southeast by a decrease in accumulation of snow. The thinning signal in the margins of the ice sheet has accelerated between 1980 and 2014, in line with observations of increased surface melt.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | During the period 1980-2014, the surface elevation increase in the interior shifted from the east and northeast towards the center of the ice sheet, and stagnated towards the end of the time series. The contribution from surface melt to interior surface lowering has increased markedly in this period, with the largest firn volume decrease due to surface melting in the extreme melt summer of 2012.
The time series of surface elevation change due to SMB and firn processes (ḣ f ) is suitable to isolate ice-dynamical thinning from altimetry-based observations of surface elevation. Combining it with the next generation of altimetry products, e.g. from Cryosat-2, allows for further improved assessment of the current imbalance of the Greenland ice sheet.  Table 1. Coordinates of analyzed elevation change points labelled in Fig. 5 with the number of repeat observations N , observed surface elevation trend from the laser altimetry, the best fit trend in the residuals between the observations and the firn model heights, and the root-mean-square of the residuals between the observations and fitted model (in cm). Trends in bold are significantly different than zero at 95 % confidence.