Elevation change of the Greenland Ice Sheet due to surface mass balance and firn processes, 1960–2014

Abstract

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 firn. The last two contributions are quantified here using a firn model that includes compaction, meltwater percolation, and refreezing. The model is forced with surface mass fluxes and temperature from a regional climate model for the period 1960–2014. The model results agree with observations of surface density, density profiles from 62 firn cores, and altimetric observations from regions where ice-dynamical surface height changes are likely small. In areas with strong surface melt, the firn model overestimates density. We find that the firn layer in the high interior is generally thickening slowly (1–5 cm yr${}^{-\mathrm{1}}$). In the percolation and ablation areas, firn and SMB processes account for a surface elevation lowering of up to 20–50 cm yr${}^{-\mathrm{1}}$. Most of this firn-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 firn and ice volume change between 1980 and 2014 is estimated at$-$3295 $\pm$ 1030 km${}^{\mathrm{3}}$ due to firn and SMB changes, corresponding to an ice-sheet average thinning of 1.96 $\pm$ 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.

Introduction

The mass balance of the Greenland Ice Sheet has been negative over the last few decades e.g.. 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 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 was assumed in earlier studies that ice-dynamical elevation changes are negligible and that volume changes are caused solely by changes in the firn layer. A fixed density was adopted to convert from volume to mass, for example a fresh-snow density, or a representative density for the entire firn layer (e.g.  and for Antarctica, and for Greenland). This assumption is still commonly made over smaller ice caps and glaciers e.g., 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 . Therefore, more recent altimetry-based mass balance estimates of (parts of) the Greenland Ice Sheet e.g. make use of firn models that take into account firn compaction and accumulation variability. More specifically, the empirical firn models used by , and take into account firn compaction, as well as the formation of ice lenses due to melt and refreezing in annual layers . The firn model used by was developed by and . It is driven by satellite-observed surface temperature and a fixed, firn-core derived relation between temperature and accumulation change. Not directly applied to altimetry-based assessments of mass loss but nonetheless very similar to our firn model, a firn model was compared by 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\mathrm{km}\times 11\mathrm{km}$. The results are obtained using a semi-empirical firn compaction model that is forced with surface mass fluxes and temperature from the polar-adapted regional climate model RACMO2.3 . 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. ), which is forced at its upper boundary by temperature and surface mass fluxes from the regional climate model RACMO, version 2.3 see Sect. and. In Sect. , we evaluate the firn model results against 62 shallow and medium-depth firn cores (details in Sect. ), and against airborne laser altimetry collected in areas of limited ice-dynamical activity (Sect. ).

The firn model

The firn model IMAU-FDM v1.0; 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 ${\dot{h}}_{\mathrm{f}}$ is given by $${\dot{h}}_{\mathrm{f}}=v_{\text{acc}}+v_{\text{snd}}+v_{\text{er}}+v_{\text{me}}+v_{\text{ice}}+v_{\text{fc}},$$ 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$ represent solid precipitation ($v_{\text{acc}}$), surface sublimation (included in $v_{\text{acc}}$), snowdrift sublimation ($v_{\text{snd}}$), snowdrift erosion ($v_{\text{er}}$), snowmelt ($v_{\text{me}}$), and firn compaction ($v_{\text{fc}}$). The solid precipitation $v_{\text{acc}}$ represents solid precipitation minus surface sublimation, and thus it is not the accumulation rate as it is usually understood. Snowdrift sublimation $v_{\text{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 . Snowdrift erosion represents the horizontal redistribution of surface snow by surface winds, and it can be positive (deposition) or 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_{\text{ice}}$ that equals the mean SMB ($v_{\text{acc}}+v_{\text{snd}}+v_{\text{er}}+v_{\text{me}}$) but is of opposite sign, over a reference period for which steady-state is assumed. The choice of this reference period is discussed in Sect. . This model setup is similar to earlier process-based models like .

An important upper-boundary condition is the density of fresh snow, ${\mathit{\rho }}_{\mathrm{0}}$. While valuable observations of density in the percolation area of the Greenland Ice Sheet are available , we choose to construct a parameterization of ${\mathit{\rho }}_{\mathrm{0}}$ based on observations from the dry-snow zone only. The observations from are averages of density over the uppermost 1 $\mathrm{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 ${\mathit{\rho }}_{\mathrm{0}}$ data from the dry-snow zone: observations from north and central Greenland from , and from the EGIG line Expédition Glaciologique Internationale au Groenland; in central Greenland. We find the following correlation ($r^{\mathrm{2}}=0.52$) between observed ${\mathit{\rho }}_{\mathrm{0}}$ and annual mean surface temperature $T_{\mathrm{s}}$ (in degrees Celsius) simulated by RACMO2.3 (see Sect. ): $${\mathit{\rho }}_{\mathrm{0}}=481.0+4.834T_{\mathrm{s}}.$$

After deposition, the fresh snow starts to become denser. The effect of dry-firn compaction is represented by $v_{\text{fc}}$ in Eq. (). Compaction is an increase of firn density $\mathit{\rho }$ in time ($t$), expressed by Eq. (4) in as $$\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}t}=C\dot{b}g({\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho })\exp \left(\frac{-E_c}{RT},+,\frac{E_g}{RT}\right).$$ Here, $C$, $E_c$, and $E_g$ are constants; $\dot{b}$ is the mean annual accumulation over a reference period (Sect. ); $g$ is the gravitational acceleration; ${\mathit{\rho }}_{\mathrm{i}}$ is the ice density of 917 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$; and $R$ is the gas constant. $T$ is the firn temperature, which varies with depth. The rate constant $C$ has a value of 0.07 for $\mathit{\rho }\le 550$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{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 .

Ratio of modelled vs. observed $z_{550}$ and $z_{830}^{*}$ (${\text{MO}}_{550}$ (red circles and fit line) and ${\text{MO}}_{830}$ (blue squares and fit line), respectively) for 24 dry firn cores (8 of which are deep enough to include $z_{830}^{*}$), as a function of mean annual accumulation (in $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$).

[Figure omitted. See PDF]

To evaluate Eq. (), we compare modelled depths of the 550 and 830 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ density layers ($z_{550}$ and $z_{830}$, respectively) to observations at 62 firn core locations around the Greenland Ice Sheet (see Sect. ). 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}^{*}$ , where $z_{830}^{*}=z_{830}-z_{550}$. Figure  shows MO as a function of $\dot{b}$. 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: $$\begin{array}{}& & {\text{MO}}_{550}=1.042-0.0916\ln \left(\dot{b}\right)\text{for}\mathit{\rho }\le 550{\text{kg m}}^{-\mathrm{3}},& & {\text{MO}}_{830}=1.734-0.2039\ln \left(\dot{b}\right)\text{for}\mathit{\rho }>550{\text{kg m}}^{-\mathrm{3}}.\end{array}$$ The coefficients in these two equations are different from a previous application in Antarctica . We use these updated coefficients to improve the fit with 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 is based is biased to lower accumulation rates than those found in Greenland.

Equation () is multiplied with the correction factors MO in Eqs. () and (), which are not allowed to be smaller than 0.25 . The accurate performance of these densification expressions is further demonstrated by fully independent comparisons against in situ and remotely sensed densification rates. However, the use of a mean value of $\dot{b}$ 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 $\dot{b}$ 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_{\mathrm{c}}$, depending on the density . 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 and .

Suppose a firn column with total depth $z_{\mathrm{i}}$ and an observed or modelled density profile $\mathit{\rho }\left(z\right)$, consisting of a mixture of air (with density ${\mathit{\rho }}_{\mathrm{a}}$), water (with density ${\mathit{\rho }}_{\mathrm{w}}$), and ice (density ${\mathit{\rho }}_{\mathrm{i}}$). Using $W_{\mathrm{c}}$ from (which can vary with depth $z$), it can be shown that the firn air content $F$ (in m) is given as $$F=\underset{\mathrm{0}}{\overset{z_{\mathrm{i}}}{\int }}\frac{{\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho }\left(z\right)}{{\mathit{\rho }}_{\mathrm{i}}/(\mathrm{1}-W_{\mathrm{c}}(z\left)\right)-{\mathit{\rho }}_{\mathrm{w}}{W}_{\mathrm{c}}\left(z\right)/(\mathrm{1}-W_{\mathrm{c}}(z\left)\right)-{\mathit{\rho }}_{\mathrm{a}}}\mathrm{d}z.$$ Here, we will approximate this general expression by assuming that ${\mathit{\rho }}_{\mathrm{a}}\ll {\mathit{\rho }}_{\mathrm{i}}$ and that $W_{\mathrm{c}}\left(z\right)=\mathrm{0}$, yielding $$F=\underset{\mathrm{0}}{\overset{z_{\mathrm{i}}}{\int }}\frac{{\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho }\left(z\right)}{{\mathit{\rho }}_{\mathrm{i}}}\mathrm{d}z.$$ As an extreme example, a uniform value of $W_{\mathrm{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 ${\dot{h}}_{\mathrm{f}}$ in the ablation area. Here, the prescribed RACMO2.3 mass fluxes determine the ice ablation rate assuming the ice density ${\mathit{\rho }}_{\mathrm{i}}$. Section  describes the method to derive the SMB-induced surface elevation change ${\dot{h}}_{\mathrm{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 , a regional climate model that is adapted to simulate climatic conditions over ice sheets. The horizontal spatial resolution of RACMO2.3 is 11 $\mathrm{km}$. Forcing data are available for the period 1 January 1960–31 December 2014, in time steps of 6 $\mathrm{h}$.

RACMO2.3 supersedes RACMO2.1 . In the new version, the cloud microphysics, surface and boundary layer turbulence, and radiation transport have been updated . 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 . The 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 offline 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 ${\dot{h}}_{\mathrm{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 the properties of neither 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.. As a result, thinning has increased sharply since the mid-1990s 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 (1 January 1960–31 December 1979) and spin up the firn column at each location with a loop over this 20-$\mathrm{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-$\mathrm{year}$ periods. In practice, equilibrium is reached when all the mass in the firn layer is refreshed once following 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 uncertainty in the calculated firn depth changes in this study stems from this assumption of reference climate. We will quantify this uncertainty in Sect. .

As a second important assumption, we set the surface elevation change ${\dot{h}}_{\mathrm{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 spin-up procedure that uses multiple loops over the reference period to reach an equilibrium state, plus 20 $\mathrm{years}$ of model integration using the same data as in the spin-up procedure. For the accumulation area, it means that the amount of mass leaving the bottom of the firn layer (with a velocity $v_{\text{ice}}$) is assumed equal to the total mass added to and retained in the firn column by snowfall and refreezing, minus run-off.

For the ablation area, the assumption of ${\dot{h}}_{\mathrm{f}}=\mathrm{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_{\text{ice}}$. In Eq. (), we set $v_{\text{ice}}$ equal to the opposite sign of the sum of all other velocity components for the reference period. Thereafter, $v_{\text{ice}}$ retains the same value, but the other parameters are free to evolve. In this framework, ${\dot{h}}_{\mathrm{f}}$ (as presented in Fig. ) 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 ${\dot{h}}_{\mathrm{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 ${\dot{h}}_{\mathrm{f}}$ should be zero over the entire simulation period, like in studies addressing the Antarctic Ice Sheet . 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.

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. ), drilled between 1995 and 2012. Among the cores are those drilled for PARCA Program for Arctic Regional Climate Assessment;; cores from the Arctic Circle Traverses ACT;; cores from the lower part of the EGIG line ; and Das 1 and Das 2 e.g..

Map of Greenland showing the 62 firn cores used for validation of the firn model. Solid circles represent firn cores longer than 30 $\mathrm{m}$; open circles are firn cores shorter than 30 $\mathrm{m}$. The inset in the lower right shows cores from in the west Greenland percolation area around 69${}^{\circ }$42${}^{\prime }$ N. Dashed contours are surface elevation isolines with a spacing of 500 $\mathrm{m}$, the uppermost being 3000 $\mathrm{m}$.

[Figure omitted. See PDF]

Vertical profiles of density from firn cores are usually based on the mean density of 0.5–2 $\mathrm{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 $\mathrm{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$\mathit{\sigma }$ errors in the elevation observations to be $\pm 10$ $\mathrm{cm}$, which includes interpolation error. Results of the comparison between the firn model and the ATM data are found in Table  and Fig. .

Coordinates of analysed elevation change points labelled in Fig.  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 $\mathrm{\%}$ confidence.

No.	Lat	Lon	Elev.	$N$	Obs.	Resid.	rms
					trend	trend	fit
	${}^{\circ }$ N	${}^{\circ }$ W	m		$\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$	$\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$	cm
1	69.90	32.4	2749	6	2.6	$-0.2$	16.7
2	71.22	35.7	3139	8	3.7	1.4	8.7
3	70.81	41.3	2935	7	5.1	6.9	38.5
4	74.36	35.5	2933	6	4.3	6.3	49.4
5	75.23	40.4	2910	6	3.5	1.0	7.7
6	75.91	45.2	2827	8	2.0	$\mathbf{-}$1.0	6.1
7	76.53	51.1	2458	8	1.6	$-0.8$	17.3
8	75.90	54.2	1971	8	$\mathbf{-}$3.1	$\mathbf{-}$5.3	13.5
9	76.32	55.3	1923	10	$\mathbf{-}$2.0	$-0.5$	22.9
10	76.92	56.2	2019	8	$\mathbf{-}$5.8	$\mathbf{-}$6.6	18.5
11	77.89	56.5	2191	7	0.4	$-$1.4	11.1
12	78.42	52.2	2239	7	1.6	0.5	9.1
13	78.04	43.7	2564	6	0.5	$-$0.5	13.7
14	79.07	44.3	2370	6	0.9	1.5	6.6
15	79.36	48.5	2170	7	3.6	2.6	21.6

Model evaluation

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 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ density levels from these cores to tune the densification parameterization in Eqs. () and (). 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.

Observed (black) and modelled (red) firn density profiles for 24 of the 62 firn cores on the Greenland Ice Sheet. The four lines of text in each panel show (1) the core name, (2) mean annual accumulation and (3) melt (in mm w.e. ${\mathrm{yr}}^{-\mathrm{1}}$) from RACMO2.3, and (4) the ratio of these fluxes ($R_{\text{ma}}$, dimensionless). Remaining profiles are shown in Appendix A as Figs. A1 and A2.

[Figure omitted. See PDF]

Figure  shows the observed and modelled density profiles for all core locations. Each panel includes the mean accumulation and melt from RACMO2.3 (in $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$) for 1960–2014, and the ratio $R_{\text{ma}}$ of these melt and accumulation averages.

The vertical resolution of the firn core data is typically 0.5–2.0 $\mathrm{m}$, thereby smoothing out the effect of ice lenses and higher-density layers. The model data in Fig.  are shown at full resolution, i.e. with layers of 5–10 $\mathrm{cm}$ thickness. The high-density layers usually represent thick layers of refrozen meltwater with a density close to that of solid ice.

Modelled vs. observed firn air content $F$ (in m) for 59 of the 62 firn cores (cores H2-1, H3-1, and H4-1 are not shown). The colour scale to the right indicates the ratio of modelled mean melt and accumulation fluxes $R_{\text{ma}}$, and the colour of each dot in the scatter plot corresponds to the value of $R_{\text{ma}}$ for that firn core location.

[Figure omitted. See PDF]

Up to an $R_{\text{ma}}$ value of 0.3–0.4, the agreement between the firn model and the observations is good. But for higher $R_{\text{ma}}$, the firn model starts to overestimate the density throughout the firn column. Figure  compares the observed and modelled firn air content $F$, showing $R_{\text{ma}}$ in colour. The model bias clearly increases for higher $R_{\text{ma}}$. This means that there are three possible causes for the misfit, which 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.

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${}^{\circ }$00${}^{\prime }$ N 47${}^{\circ }$01${}^{\prime }$ W, 1850 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) is indeed about 20 $\mathrm{\%}$ smaller than simulated by RACMO2.3 . Further north, find the equilibrium line around the EGIG line at about 1100–1200 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, while the equilibrium line altitude in RACMO2.3 is at $\sim 1650$ $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, about 45 $\mathrm{km}$ further inland. For firn cores H1-1 down to H5-1 (Fig. ), 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 . 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 run-off 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.

Upper left panel: map of elevation change points observed by the airborne laser altimeter ATM, with colour scale giving the 20-$\mathrm{year}$ trend in elevation ($\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$). The 15 panels (numbers in the upper left corner of each panel correspond to the numbered locations on the map) show surface elevation change from (black dots with 1$\mathit{\sigma }$ error bars) ATM lidar and (blue curves) firn thickness model. Green curves are the firn model results adjusted to the trend in residuals between the model and the observations.

[Figure omitted. See PDF]

Altimetry from the high-elevation interior

In the high interior of the Greenland Ice Sheet, horizontal surface ice velocities are low generally less than 10 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$;, and elevation changes resulting from ice-dynamical effects are expected to be small. Figure  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 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$ over the altimetry record (map in Fig. , Table ). 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. ) 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 modelled 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 modelled 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 . 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 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$, which accounts for more than the 6 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$ of observed increases, indicating thickening. At site 6, a negative trend in residuals of 1 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$ opposes the 2 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$ of observed surface raising, suggesting opposing contributions from dynamics 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 modelled surface elevations provides the linear contribution in ice dynamics plus the error, the error is then assessed by the root mean square (rms) of the detrended residuals (Table ). This is equivalent to adjusting the firn model time series to the ice-dynamical trend (shown as green curves in the panels of Fig. ) and computing the difference from the observations. The mean rms error is 17.4 $\mathrm{cm}$, which is close to the lidar observational uncertainty ($\sim 10$ $\mathrm{cm}$). Sites 3 and 4 have the largest errors, reaching 1.7 and 2.6 standard deviations, respectively. These sites are located at similar elevations (2930 $\mathrm{m}$) on the central eastern portion of the ice sheet, where altimetry shows steadily increasing elevations while the firn model predicts an initial increase in firn thickness until about 2005 and then a decrease thereafter.

Elevation change due to firn and SMB

Firn air content

Figure  shows the modelled firn air content $F$ on 1 September 2014. As noted in Sect. , these values are probably realistic in the dry interior and the upper part of the percolation area. In the lower percolation area, where the annual melt flux exceeds $\sim 30$ $\mathrm{\%}$ of the annual accumulation rate, the modelled firn air content is likely underestimated. Around the central dome, we find the highest values of $F$ of about 25 $\mathrm{m}$. 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.

Modelled firn air content $F$ on 1 September 2014 (in m). Dashed contour lines at 500 $\mathrm{m}$ height intervals.

[Figure omitted. See PDF]

Modelled average firn thickness change $\dot{h}$ (in $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$) for three periods: (a) 1980–2014, (b) 1980–1995 and (c) 1995–2014. The equilibrium line (according to the RACMO2.3 SMB) is shown as a thick black line in (a). Note that the colour scale is asymmetric around 0. In the ablation area, where no firn layer is present, ice thickness change by surface processes is presented.

[Figure omitted. See PDF]

Trends

By adding up all the velocity components in Eq. (), we find ${\dot{h}}_{\mathrm{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  shows ${\dot{h}}_{\mathrm{f}}$ (in $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$) for the periods 1980–2014, 1980–1995, and 1995–2014. 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 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$, or 34–170 $\mathrm{cm}$, over the entire 34-$\mathrm{year}$ period. The marginal thinning rates are much larger; they can be up to 20–50 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$, or 6–18 $\mathrm{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. b) 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–2014 then shows moderate thickening up to 5 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$ in the eastern and northern interior. Thinning occurs all around the margins, with the smallest rates (0–15 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$) in the northern and eastern coastal regions. The largest rates (exceeding 40 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$) occur in the southeast and in the western ablation area.

Decomposing the trends

The firn model allows for a decomposition of the ${\dot{h}}_{\mathrm{f}}$ signal into its velocity components (Eq. ). The upper panels in Fig.  show this decomposition for the period 1980–2014. The thickening in the eastern interior (Fig. a) can be almost fully ascribed to a positive accumulation anomaly (Fig. a), offset by a small increase in firn compaction due to this extra firn (Fig. g). In the lower accumulation area, the positive accumulation anomaly is offset by a significant increase in surface melt, giving zero or slightly negative ${\dot{h}}_{\mathrm{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.

Anomaly of vertical velocity components with respect to the reference period 1960–1979 (in $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$). (a–c) Accumulation velocity $v_{\text{acc}}$; (d–f) melt velocity $v_{\text{me}}$; (g–i) firn compaction velocity $v_{\text{fc}}$; (j–l) snowdrift sublimation, erosion and deposition velocities ($v_{\text{snd}}+v_{\text{er}}$). Note that the colour scales for $v_{\text{me}}$ is strongly asymmetric around 0.

[Figure omitted. See PDF]

For the period 1980–1995, the accumulation anomaly differed from the 1995–2014 period, as shown in panels b and c of Fig. . 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–2013 shows 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 ${\dot{h}}_{\mathrm{f}}$, melt and firn compaction, are negative by definition. To complement this, 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. . 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 $\mathrm{\%}$, 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.

Fraction of surface lowering (negative velocities only) during 1980–2014 caused by (a) melt; (b) firn compaction; (c) snowdrift processes; and (d) sublimation. Up to a fraction of 0.15, the colour scale is divided in steps of 0.025. Above 0.15, the step size is 0.05.

[Figure omitted. See PDF]

Error estimate

An important source of uncertainty in ${\dot{h}}_{\mathrm{f}}$ is the steady state assumed for the spin-up of the firn model. As explained in Sect. , 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 $\mathrm{years}$, and an accumulation increase by 12 $\mathrm{\%}$ over the period 1600–2009 . The mean reconstructed, ice-sheet-wide accumulation over this period is 782 $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$. For 1960–1979, it is 786 $\mathrm{Gt}{\mathrm{yr}}^{-\mathrm{1}}$, i.e. 0.5 $\mathrm{\%}$ larger than the long-term average. Other reconstructions show the 1960–1979 solid precipitation flux to be 0–10 $\mathrm{\%}$ different from the period between $\sim 1850$ and the present day . For the error analysis, we therefore assume that the 1960–1979 accumulation fluxes (precipitation minus evaporation/sublimation) can differ by up to 15 $\mathrm{\%}$ from the long-term accumulation history. We denote this accumulation uncertainty by ${\mathit{\sigma }}_{\dot{b}}$ (in $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$).

For snowmelt, we also assume a maximum deviation of 15 $\mathrm{\%}$ of the 1960–1979 melt flux compared to the long-term history. There is limited evidence for this, but according to reconstructions from the run-off flux for 1960–1979 differs about 5–10 $\mathrm{\%}$ from the 1870–2010 mean. This uncertainty is given as ${\mathit{\sigma }}_{\dot{m}}$ ($\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$).

For the suite of 62 firn core locations, we perform four sensitivity tests, in which we increase and decrease melt or accumulation by 15 $\mathrm{\%}$ after completion of the spin-up. The resulting drift in ${\dot{h}}_{\mathrm{f}}$ for 1960–2014 gives an error in the firn thickness change. The error in surface elevation change due to ${\mathit{\sigma }}_{\dot{b}}$ is written as ${\mathit{\sigma }}_{\dot{h},\dot{b}}$. Analogously, the error in surface elevation change due to ${\mathit{\sigma }}_{\dot{m}}$ is given as ${\mathit{\sigma }}_{\dot{h},\dot{m}}$. By relating ${\mathit{\sigma }}_{\dot{h},\dot{b}}$ and ${\mathit{\sigma }}_{\dot{h},\dot{m}}$ empirically to accumulation and melt at the core locations, we can expand the error product over the entire ice sheet. These relations are $$\begin{array}{}& & {\mathit{\sigma }}_{\dot{h},\dot{b}}={\mathit{\sigma }}_{\dot{b}}(0.107+3.609\times {10}^{\mathrm{4}}\dot{b}),& & {\mathit{\sigma }}_{\dot{h},\dot{m}}={\mathit{\sigma }}_{\dot{m}}(0.225+1.064\times {10}^{\mathrm{3}}\dot{m}),\end{array}$$ with $\dot{b}$ and $\dot{m}$ in millimetres per year.

The uncertainties ${\mathit{\sigma }}_{\dot{b}}$ and ${\mathit{\sigma }}_{\dot{m}}$ cannot be regarded as independent. The SMB module in RACMO2.3 contains interactions between accumulation and melt. We identify the melt-albedo 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 to be dependent, and we add them up linearly (${\mathit{\sigma }}_{\dot{h}}={\mathit{\sigma }}_{\dot{h},\dot{b}}+{\mathit{\sigma }}_{\dot{h},\dot{m}}$).

The resulting total error and its components are shown in Fig. . 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 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$ in the interior to 10–20 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$ in the lower percolation area. The largest total error (more than 40 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$) is found in the southeast, where high snowfall rates coincide with large amounts of melt.

Estimate of errors in firn thickness change ($\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$). (a) Error ${\mathit{\sigma }}_{\dot{h},\dot{b}}$ due to accumulation uncertainty ${\mathit{\sigma }}_{\dot{b}}$; (b) error ${\mathit{\sigma }}_{\dot{h},\dot{m}}$ due to melt uncertainty ${\mathit{\sigma }}_{\dot{m}}$; and (c) total error ${\mathit{\sigma }}_{\dot{h}}$. Note the non-linear colour scale.

[Figure omitted. See PDF]

For low and medium values of $\dot{b}$, the term between brackets in Eq. () 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. () 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 $\dot{b}>2427$ $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$, the term between brackets in Eq. () becomes larger than unity. There is no physical explanation for this behaviour: it is caused by the empirical nature of Eq. (). However, we decided not to cap the uncertainty, for the following reason: the densification rate in Eq. () depends on a 20-$\mathrm{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 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.

(a) Ice-sheet-integrated volume change for 1980–2014 (in ${\mathrm{km}}^{\mathrm{3}}$) due to firn and SMB processes. The thin red line is the volume change based on weekly output from the firn model; the thick red line is a 1-$\mathrm{year}$ running average; and the pink shaded area shows the uncertainty estimate. The blue line shows the volume change for the part of the ice-sheet surface elevated more than 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (b) Partitioning of the total volume change into the three main components: accumulation $v_{\text{acc}}$, melt $v_{\text{me}}$, and firn compaction $v_{\text{fc}}$. (c) Like (b) but for the ice sheet higher than 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$

[Figure omitted. See PDF]

Integrated volume change

Figure  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\pm 1030$ ${\mathrm{km}}^{\mathrm{3}}$ due to firn and SMB processes alone. Averaged over the ice sheet, this is a mean surface lowering of $1.96\pm 0.61$ $\mathrm{m}$. Almost all of this total volume loss took place in the part of the Greenland Ice Sheet where the surface is under 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$

Panels b and c of Fig.  show the partitioning of the volume change for the entire ice sheet and the part elevated above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, respectively. Over 1980–2014, the volume loss through melt was slightly over 4700 ${\mathrm{km}}^{\mathrm{3}}$. This loss was partly compensated for by an increase in snow accumulation of about 1500 ${\mathrm{km}}^{\mathrm{3}}$. On most of the ice sheet, firn compaction has accelerated (become more negative) due to an increase in accumulation. Above 2000 $\mathrm{m}$, the effect is clearly visible (Fig. c). Integrated over the ice sheet however, we see a small slowdown in firn compaction, corresponding to about 500 ${\mathrm{km}}^{\mathrm{3}}$ (Fig. b). 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. h and i).

Up to about 2005, firn volume change above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ was dominated by accumulation variability (consistent with e.g. ), and below 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{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 2005 however, the total firn volume above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ has started to decrease, mainly because surface melt has migrated inland e.g., but also because the accumulation increase, clearly visible between 1980 and 2000, stagnated 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 is clearly visible in the results of the firn model. A large drop in total volume of 1386 ${\mathrm{km}}^{\mathrm{3}}$ is seen in the summer of 2012, of which 1150 ${\mathrm{km}}^{\mathrm{3}}$ is contributed to melt. Melt in the part of the ice sheet above 2000 $\mathrm{m}$ contributed almost one third (371 ${\mathrm{km}}^{\mathrm{3}}$) to this volume anomaly. For perspective, the volume loss above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ in the summer of 2012 is equal to the volume gained by snowfall in the interior in the 16 $\mathrm{years}$ between 1980 and 1996.

Altimetry correction and mean density of firn-related mass loss

The present data set 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 $\dot{h}(t_{\mathrm{0}},t_{\mathrm{1}})\equiv \left(h\right(t_{\mathrm{1}})-h(t_{\mathrm{0}}\left)\right)/(t_{\mathrm{1}}-t_{\mathrm{0}})$ between initial time $t_{\mathrm{0}}$ and time $t_{\mathrm{1}}>t_{\mathrm{0}}$. The firn model computes a surface elevation change due to firn and SMB processes ${\dot{h}}_{\mathrm{f}}(t_{\mathrm{0}},t_{\mathrm{1}})\equiv h_{\mathrm{f}}\left(t_{\mathrm{1}}\right)-h_{\mathrm{f}}\left(t_{\mathrm{0}}\right)$. The supposed ice-dynamical contribution (neglecting vertical bed movement, but this could be included; ) is then ${\dot{h}}_{\mathrm{d}}=\dot{h}-{\dot{h}}_{\mathrm{f}}$. The associated mass change ${\dot{m}}_{\mathrm{d}}$ is simply computed as $${\dot{m}}_{\mathrm{d}}={\mathit{\rho }}_{\mathrm{i}}{\dot{h}}_{\mathrm{d}}.$$ The mass change ${\dot{m}}_{\mathrm{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 ice-sheet viscosity , and transient variations in ice flow due to long-term changes in accumulation. The mass change due to the SMB, ${\dot{m}}_{\mathrm{f}}$, is computed directly from the RACMO2.3 forcing, using SMB anomalies with respect to the appropriate reference period (1960–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 $\dot{m}$ is then computed as $$\dot{m}={\dot{m}}_{\mathrm{d}}+{\dot{m}}_{\mathrm{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, due not only to firn processes but also 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 2 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 $\mathrm{\%}$) 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 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{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 $\mathrm{cm}{\mathrm{yr}}^{-\mathrm{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.

During the period 1980–2014, the surface elevation increase in the interior shifted from the east and northeast towards the centre 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 (${\dot{h}}_{\mathrm{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.

Remaining vertical profiles of firn density

Here, we present the remaining 33 vertical profiles of firn density, in addition to the 24 firn cores shown in Fig. . They are shown in Figs. A1 and A2.

Observed (black) and modelled (red) firn density profiles for 24 of the 62 firn cores on the Greenland Ice Sheet. The four lines of text in each panel show (1) the core name, (2) mean annual accumulation and (3) melt (in mm w.e. ${\mathrm{yr}}^{-\mathrm{1}}$) from RACMO2.3, and (4) the ratio of these fluxes ($R_{\text{ma}}$, dimensionless). The other vertical profiles are shown in Figs.  and A2.

[Figure omitted. See PDF]

Observed (black) and modelled (red) firn density profiles for 9 of the 62 firn cores on the Greenland Ice Sheet. The four lines of text in each panel shows (1) the core name, (2) mean annual accumulation and (3) melt (in mm w.e. ${\mathrm{yr}}^{-\mathrm{1}}$) from RACMO2.3, and (4) the ratio of these fluxes ($R_{\text{ma}}$, dimensionless). The other vertical profiles are shown in Figs.  and A1.

[Figure omitted. See PDF]

P. Kuipers Munneke carried out the model experiments and analysis, and wrote the paper; S. R. M. Ligtenberg co-developed the firn model; B. P. Y. Noël provided the RACMO2.3 forcing data; I. M. Howat provided ATM altimetry data and contributed to the writing; J. E. Box contributed firn core data and historical accumulation data; E. Mosley-Thompson, J. R. McConnell, K. Steffen, J. T. Harper and S. B. Das provided firn core data; and M. R. van den Broeke conceived this study. All authors reviewed and commented on the manuscript.

Acknowledgements

P. Kuipers Munneke received financial support from the Netherlands Polar Programme (NPP) of the Netherlands Institute for Scientific Research (NWO). ECMWF at Reading (UK) is acknowledged for use of the Cray supercomputing system. Graphics were made using the NCAR Command Language (version 6.2.1, 2014). The J. E. Box contribution is supported by Det Frie Forskningsråd grant 4002-00234 and Geocenter Denmark. Edited by: J. L. Bamber