Quantifying the snowmelt–albedo feedback at Neumayer Station, East Antarctica

Abstract

We use 24 years (1992–2016) of high-quality meteorological observations at Neumayer Station, East Antarctica, to force a surface energy balance model. The modelled 24-year cumulative surface melt at Neumayer amounts to 1154 mm water equivalent (w.e.), with only a small uncertainty ($\pm \mathrm{3}$ mm w.e.) from random measurement errors. Results are more sensitive to the chosen value for the surface momentum roughness length and new snow density, yielding a range of 900–1220 mm w.e. Melt at Neumayer occurs only in the months November to February, with a summer average of 50 mm w.e. and large interannual variability ($\mathit{\sigma }=\mathrm{42}$ mm w.e.). This is a small value compared to an annual average (1992–2016) accumulation of $\mathrm{415}\pm \mathrm{86}$ mm w.e. Absorbed shortwave radiation is the dominant driver of temporal melt variability at Neumayer. To assess the importance of the snowmelt–albedo feedback we include and calibrate an albedo parameterisation in the surface energy balance model. We show that, without the snowmelt–albedo feedback, surface melt at Neumayer would be approximately 3 times weaker, demonstrating how important it is to correctly represent this feedback in model simulations of surface melt in Antarctica.

1 Introduction

The Antarctic ice sheet (AIS) contains more than 25 million km${}^{\mathrm{3}}$ of ice, sufficient to raise global mean sea level by almost 60 m if melted completely . Between 1992 and 2017, the AIS lost mass at an accelerated rate, contributing $\mathrm{7.6}\pm \mathrm{3.9}$ mm to global sea level . This mass loss is mainly observed in coastal West Antarctica and the Antarctic Peninsula (AP) and is caused by glaciers that accelerated after their buttressing ice shelves had thinned or disintegrated . The interaction between meltwater and firn, the intermediate product between snow and glacier ice, is hypothesised to play an important role in ice shelf disintegration . If the firn layer contains enough air, as is the case for most of the AIS, meltwater can percolate downwards and refreeze . If the storage capacity of the firn layer is reduced, surface meltwater can flow laterally towards the ice shelf edge , be stored englacially or form ponds on the ice shelf surface . In all cases, meltwater can accumulate in crevasses, thereby increasing the hydrostatic pressure in the crevasse tip, causing it to penetrate farther down. When a crevasse reaches the bottom of the ice shelf or a basal crevasse, part of the ice shelf disintegrates, a process called hydrofracturing . Hydrofracturing has been identified as a potential precursor of rapid loss of Antarctic ice, accelerating sea level rise . In combination with enhanced ocean swell under low sea-ice conditions , hydrofracturing likely caused the disintegration of the Larsen B ice shelf in the AP in 2002 . In July 2017, a large iceberg calved from the Larsen C ice shelf, but it is unclear whether this signifies a further southward progression of ice shelf destabilisation in the AP .

Improving our predictive capabilities of future ice shelf stability, AIS mass loss and associated sea level rise thus requires a thorough understanding of the surface melt process on Antarctic ice shelves. In contrast to meltwater occurrence, which is readily observed from space , observational estimates of surface melt rates on Antarctic ice shelves are rare; they have been obtained locally through explicit modelling of the surface energy balance (SEB) . In turn, these enabled continent-wide melt rate estimates using calibrated satellite products based on backscatter strength of radio waves . These studies invariably demonstrate that, in most parts of Antarctica, melt is currently a weak and intermittent process. In this melt regime, the positive snowmelt–albedo feedback (SMAF) plays a decisive role: when snow melts, meltwater may refreeze in the cold snowpack, resulting in considerably larger grains ($\sim \mathrm{1}$ mm) than new snow or snow that has been subjected to only dry compaction ($\sim \mathrm{0.1}$ mm). Larger snow grains reduce backward scattering of photons into the snowpack, increasing the probability of absorption and reducing the surface albedo, especially in the near-infrared . This further enhances absorption of solar radiation and melt. For pure, uncontaminated snow, the strength of the SMAF depends on multiple factors, e.g. the intensity and duration of the melt and the frequency and intensity of snowfall events, which provide new snow consisting of smaller grains. We therefore expect the SMAF to be spatially and temporally variable on Antarctic ice shelves.

Most studies on the SMAF address the removal of (seasonal) snow and the appearance of dark soil or water , leading to further warming of the air and water. These studies commonly express the melt–albedo feedback in terms of air and water temperature sensitivity. Our aim is to quantify the impact on the melt rate of the darkening but not the disappearance of snow, a process addressed by far fewer studies . To that end, we implement a snow albedo parameterisation in an SEB model, which is then calibrated using observations and used to study the sensitivity of melt rates to snow properties that influence snow albedo. We use 24 years of high-quality in situ observations from the German research station Neumayer (Fig. ) to calculate the SEB and melt rate. We investigate the effects of measurement uncertainties and model settings on the modelled cumulative amount of surface melt. We then analyse the main drivers of surface melt and the magnitude of the SMAF at Neumayer by switching the feedback process in the albedo parameterisation on and off.

Figure 1

Map of the Antarctic continent. The red cross indicates the location of Neumayer Station. Imagery (©) 2016 DigitalGlobe, Inc.

[Figure omitted. See PDF]

The SEB model is explained in Sect. , followed by a description of the albedo parameterisation in Sect. . The meteorological data used to force the SEB model are described in Sect. . The results section is split into two parts: in Sect.  we present and discuss the SEB and melt rate that are obtained using the observed albedo. In Sect.  the albedo parameterisation is used instead and the SMAF is quantified and discussed. Finally, the results are discussed in Sect. .

2 Methods

2.1 Surface energy balance model

The one-dimensional energy balance model is a further development of the models presented by , , and ; here only the main features are described. The energy balance of an infinitesimally thin surface layer (the “skin” layer) is defined as follows:

1$$M=\mathrm{SW}\downarrow +\mathrm{SW}\uparrow +\mathrm{LW}\downarrow +\mathrm{LW}\uparrow +Q_{\mathrm{S}}+Q_{\mathrm{L}}+Q_G,$$ where positive fluxes are defined to be directed towards the surface. $\mathrm{SW}\downarrow$ and $\mathrm{SW}\uparrow$ are the incoming and reflected shortwave radiation, $\mathrm{LW}\downarrow$ and $\mathrm{LW}\uparrow$ are the downward and upward longwave radiation, $Q_{\mathrm{S}}$ and $Q_{\mathrm{L}}$ the turbulent sensible and latent heat fluxes and $Q_{\mathrm{G}}$ is the conductive subsurface heat flux. We neglect latent energy from rain. $M$ is the energy used to melt snow or ice and is non-zero only when the surface has reached the melting point of ice ($T_{\mathrm{s}}=\mathrm{273.15}\mathrm{K}$). Throughout this paper, melt and accumulation amounts are expressed in terms of millimetre water equivalent (mm w.e.), which equals kg m${}^{-\mathrm{2}}$. In order to calculate $Q_{\mathrm{G}}$ and allow for densification, meltwater percolation and refreezing, a snow–firn model is used, initialised with 70 layers. The layer thickness varies from 1 cm at the top to 2 m at the bottom (25 m depth). We impose a no-energy flux boundary condition at the lowermost model level. New snow density is parameterised following the expression of , which relates it to the prevailing surface temperature ($T_{\mathrm{s}}$) and 10 m wind speed ($V_{\mathrm{10}\mathrm{m}}$) and imposes a lower limit of new snow density ${\mathit{\rho }}_{\mathrm{s},\mathrm{0}}$. Meltwater percolation is based on the tipping-bucket method e.g., allowing for immediate downward transport (within a single timestep of 10 s) of remaining water if a layer has attained its maximum capillary retention, as modelled using the expressions of . Meltwater refreezing increases the density and temperature of a layer. At the bottom of the firn layer, the meltwater is assumed to run off immediately, i.e. the model does not allow for slush/superimposed ice formation or lateral water movement. Turbulent fluxes are calculated following the “bulk” method, which is based on Monin–Obukhov similarity theory (see e.g. for relevant equations) between a single measurement level (2 m for temperature and humidity, 10 m for wind) and the surface, assuming the latter to be saturated with respect to ice and using the stability functions according to for unstable and for stable conditions.

Subsurface penetration of shortwave radiation is calculated using a spectral model , based on the parameterisation by , which is in turn based on the two-stream radiation model of . The impact on modelled melt and the quantification of the SMAF is discussed in the relevant sections.

The terms in Eq. () are either based on observations or can be expressed as a function of the skin temperature $T_{\mathrm{s}}$. The SEB is solved iteratively by looking for a value of $T_{\mathrm{s}}$ that closes the SEB to within 0.005 K between iterations: if $T_{\mathrm{s}}>\mathrm{273.15}\mathrm{K}$, it is reset to 273.15 K and excess energy $M$ is used for surface melt. To evaluate model performance, the modelled value of $T_{\mathrm{s}}$ is compared to observed $T_{\mathrm{s}}$ calculated from $\mathrm{LW}\uparrow$, using Stefan–Boltzmann's law for a longwave emissivity $\mathit{ϵ}=\mathrm{1}$: 2$$\mathrm{LW}\uparrow =\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}},$$ where $\mathit{\sigma }=\mathrm{5.67}\cdot {\mathrm{10}}^{-\mathrm{8}}$ W m${}^{-\mathrm{2}}$ K${}^{-\mathrm{4}}$ is the Stefan–Boltzmann constant.

Surface roughness lengths for momentum, heat and moisture are related through the expression of : 3$$\ln \left({\displaystyle \frac{z_{\mathrm{0}*}}{z_{\mathrm{0},\mathrm{m}}}}\right)=a_{\mathrm{1}}+a_{\mathrm{2}}\ln \left({\mathit{\text{Re}}}_{*}\right)+a_{\mathrm{3}}\ln {\left({\mathit{\text{Re}}}_{*}\right)}^{\mathrm{2}},$$ where $z_{\mathrm{0}*}$ represents either $z_{\mathrm{0},\mathrm{h}}$ or $z_{\mathrm{0},\mathrm{q}}$, the roughness lengths for heat and moisture respectively, and $a_{\mathrm{1}}$, $a_{\mathrm{2}}$ and $a_{\mathrm{3}}$ are coefficients determined by for various regimes of the roughness Reynolds number ${\mathit{\text{Re}}}_{*}=\frac{u_{*}{z}_{\mathrm{0},\mathrm{m}}}{\mathit{\nu }}$ with kinematic viscosity $\mathit{\nu }$ and friction velocity $u_{*}$.

2.2 Albedo parameterisation

Because the shortwave radiation sensor faces the sky and includes a significant direct component, measured $\mathrm{SW}\downarrow$ suffers from relatively large uncertainties owing to poor sensor cosine response, sensor tilt and/or rime formation . In order to improve the accuracy of observed net shortwave radiation used in the SEB calculations (Sect. ), we calculate ${\mathrm{SW}}_{\mathrm{net}}$ based on $\mathrm{SW}\uparrow$, which is diffuse and hence much less sensitive to these errors. To further decrease the impact of these errors, we use a 24 h moving average albedo, as described in . In Sect. , in which albedo is parameterised to study melt–albedo feedbacks, for consistency we use measured $\mathrm{SW}\uparrow$ in combination with parameterised albedo to estimate ${\mathrm{SW}}_{\mathrm{net}}$.

In Sect. , the parameterised surface albedo $\mathit{\alpha }$ is described as a base albedo ${\mathit{\alpha }}_S$, modified by perturbations describing the effect of changing solar zenith angle $\mathit{\theta }$ ($\mathrm{d}{\mathit{\alpha }}_u$), the cloud optical thickness $\mathit{\tau }$ ($\mathrm{d}{\mathit{\alpha }}_{\mathit{\tau }}$) and the concentration of black carbon in the snow ($\mathrm{d}{\mathit{\alpha }}_c$) :

4$$\mathit{\alpha }={\mathit{\alpha }}_S+\mathrm{d}{\mathit{\alpha }}_u+\mathrm{d}{\mathit{\alpha }}_{\mathit{\tau }}+\mathrm{d}{\mathit{\alpha }}_c.$$ For Antarctica, we neglect the impact of impurities in the snow ($\mathrm{d}{\mathit{\alpha }}_c=\mathrm{0}$); $\mathrm{d}{\mathit{\alpha }}_u$ and $\mathrm{d}{\mathit{\alpha }}_{\mathit{\tau }}$ both depend on the base albedo ${\mathit{\alpha }}_S$, $\mathrm{d}{\mathit{\alpha }}_u$ in addition depends on the solar zenith angle ($u=\cos \mathit{\theta }$) and $\mathrm{d}{\mathit{\alpha }}_{\mathit{\tau }}$ on the cloud optical thickness $\mathit{\tau }$: $$\begin{array}{}\text{5}& {\displaystyle }& {\displaystyle }\mathrm{d}{\mathit{\alpha }}_u=\mathrm{0.53}{\mathit{\alpha }}_S(\mathrm{1}-{\mathit{\alpha }}_S)(\mathrm{1}-\mathrm{0.64}x-(\mathrm{1}-x)u{)}^{\mathrm{1.2}},\text{6}& {\displaystyle }& {\displaystyle }\mathrm{d}{\mathit{\alpha }}_{\mathit{\tau }}={\displaystyle \frac{\mathrm{0.1}\mathit{\tau }({\mathit{\alpha }}_S+\mathrm{d}{\mathit{\alpha }}_c{)}^{\mathrm{1.3}}}{(\mathrm{1}+\mathrm{1.5}\mathit{\tau }{)}^{{\mathit{\alpha }}_S}}},\end{array}$$ where $x=min\left(\frac{\sqrt{\mathit{\tau }}}{\mathrm{3}u},\mathrm{1}\right)$. The base albedo depends on the snow grain size $r_{\mathrm{e}}$ (in metres): 7$${\mathit{\alpha }}_S=\mathrm{1.48}-\mathrm{1.27048}{r}_{\mathrm{e}}^{\mathrm{0.07}},$$ in which the snow grain size $r_{\mathrm{e}}$ on time step $t$ is parameterised as 8$$r_{\mathrm{e}}\left(t\right)=\left[r_{\mathrm{e}}(t-\mathrm{1})+\mathrm{d}{r}_{\mathrm{e},\mathrm{dry}}+\mathrm{d}{r}_{\mathrm{e},\mathrm{wet}}\right]f_{\mathrm{o}}+r_{\mathrm{e},\mathrm{0}}{f}_{\mathrm{n}}+r_{\mathrm{e},\mathrm{r}}{f}_{\mathrm{r}}.$$ Here, $\mathrm{d}{r}_{\mathrm{e},\mathrm{dry}}$ and $\mathrm{d}{r}_{\mathrm{e},\mathrm{wet}}$ describe the metamorphism of dry and wet snow respectively, $f_{\mathrm{o}}$, $f_{\mathrm{n}}$ and $f_{\mathrm{r}}$ are the fractions of old, new and refrozen snow, and $r_{\mathrm{e},\mathrm{0}}$ and $r_{\mathrm{e},\mathrm{r}}$ are the grain sizes of new and refrozen snow. Dry snow metamorphism is parameterised following : 9$${\displaystyle \frac{\mathrm{d}{r}_{\mathrm{e},\mathrm{dry}}}{\mathrm{d}t}}={\left({\displaystyle \frac{\mathrm{d}{r}_{\mathrm{e}}}{\mathrm{d}t}}\right)}_{\mathrm{0}}{\left({\displaystyle \frac{\mathit{\eta }}{(r_{\mathrm{e}}-r_{\mathrm{e},\mathrm{0}})+\mathit{\eta }}}\right)}^{\mathrm{1}/\mathit{\kappa }},$$ where $r_{\mathrm{e},\mathrm{0}}$ is the new snow grain size, and the coefficients ${\left(\frac{\mathrm{d}{r}_{\mathrm{e}}}{\mathrm{d}t}\right)}_{\mathrm{0}}$, $\mathit{\eta }$ and $\mathit{\kappa }$ are obtained from a look-up table. This look-up table is compiled based on simulations with the SNICAR model , which calculates the snow metamorphism resulting from temperature gradient metamorphism. $\mathrm{d}{r}_{\mathrm{e},\mathrm{wet}}$ is a function of the snow grain size $r_{\mathrm{e}}$ itself and the liquid water content $f_{\mathrm{liq}}$ : 10$${\displaystyle \frac{\mathrm{d}{r}_{\mathrm{e},\mathrm{wet}}}{\mathrm{d}t}}={\displaystyle \frac{C{f}_{\mathrm{liq}}^{\mathrm{3}}}{\mathrm{4}\mathit{\pi }{r}_{\mathrm{e}}^{\mathrm{2}}}},$$ where $C$ is a constant ($\mathrm{4.22}\cdot {\mathrm{10}}^{-\mathrm{13}}$ m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$).

The fractions $f_{\mathrm{o}}$, $f_{\mathrm{n}}$ and $f_{\mathrm{r}}$ are derived from the snow/firn model, and the grain sizes of new and refrozen snow are constants; the method for determining their values from a tuning exercise is described in Sect. .

Figure 2

(a) Downward longwave radiation vs. air temperature. The red lines are quadratic fits of the upper and lower 5 percentile boundaries. The longwave-equivalent cloud cover is determined by linear interpolation between these bounds. (b) Optical thickness vs. cloud cover. The red line resembles the best fit to a function $\mathit{\tau }=c_{\mathrm{1}}\left(e^{c_{\mathrm{2}}{N}_{\mathit{ϵ}}}-\mathrm{1}\right)$. The shaded area indicates the 95 % uncertainty range.

[Figure omitted. See PDF]

To determine cloud optical thickness $\mathit{\tau }$, an empirical relation between $\mathit{\tau }$ and the longwave-equivalent cloud cover $N_{\mathit{ϵ}}$ is used following : 11$$\mathit{\tau }=c_{\mathrm{1}}\left(\exp \left(c_{\mathrm{2}}{N}_{\mathit{ϵ}}\right)-\mathrm{1}\right),$$ with fitting parameters $c_{\mathrm{1}}$ and $c_{\mathrm{2}}$. $N_{\mathit{ϵ}}$ is determined using a method described by , which relates hourly values of downward longwave radiation $\mathrm{LW}\downarrow$ to near-surface air temperature $T_{\mathrm{2}\mathrm{m}}$ as illustrated in Fig. a. Red lines indicate quadratic fits through the upper and lower 5 percentiles of the data, assumed to represent fully cloudy and clear conditions, respectively. $N_{\mathit{ϵ}}$ is obtained by linearly interpolating between these upper and lower bounds, yielding values between 0 and 1. Hourly values for cloud cover are then used to obtain values for $\mathit{\tau }$ (Fig. b). The values used for the fit parameters $c_{\mathrm{1}}=\mathrm{5.404}$ and $c_{\mathrm{2}}=\mathrm{2.207}$ (both dimensionless) differ somewhat from , who used daily values for the fit.

2.3 Observational data

The SEB model is forced with data from the meteorological observatory at the German research station Neumayer, situated on the Ekström ice shelf . The observatory has been operational since 1981 and was relocated in 1992 and 2009. In 2016, its location was 70${}^{\circ }$40${}^{\prime }$ S, 8${}^{\circ }$16${}^{\prime }$ W (Fig. ). The observatory is one of only four Antarctic stations – and the only one situated on an ice shelf – that is part of the Baseline Surface Radiation Network (BSRN), a global network of stations with high-quality (artificially ventilated) radiation observations, coordinated by the Alfred Wegener Institute (AWI). We use hourly averages of 2 m temperature ($T_{\mathrm{2}\mathrm{m}}$) and specific humidity ($q_{\mathrm{2}\mathrm{m}}$), 10 m wind speed ($V_{\mathrm{10}\mathrm{m}}$), surface pressure ($p$) and radiation fluxes for the period April 1992–January 2016 (24 years) to force the SEB model; their uncertainty ranges are provided in Table . Approximately 4.1 % of the data points contained at least one missing variable, which mostly come from daily performed visual inspection of the data. To obtain a continuous data set, all missing data were replaced: pressure, relative humidity, wind speed, temperature and longwave radiation were simply linearly interpolated. In the case of shortwave radiation, the missing value was replaced by imitating the average daily cycle of the 2 preceding days. As the measurement station is visited and maintained every day, the impact of rime formation is limited, as is the tilt of the observation mast, resulting in a high-quality meteorological data set.

Table 1

Listing of used measurement variables and their associated measurement uncertainties.

Variable	Uncertainty range
$V_{\mathrm{10}\mathrm{m}}$	max (0.5 m s${}^{-\mathrm{1}}$, 5 %)
$\mathrm{SW}\downarrow$	5 W m${}^{-\mathrm{2}}$
$\mathrm{SW}\uparrow$	5 W m${}^{-\mathrm{2}}$
$\mathrm{LW}\downarrow$	5 W m${}^{-\mathrm{2}}$
$\mathrm{LW}\uparrow$	5 W m${}^{-\mathrm{2}}$
$T_{\mathrm{2}\mathrm{m}}$	0.1 ${}^{\circ }$C
${\mathrm{RH}}_{\mathrm{2}\mathrm{m}}$	5 %
$p$	0.5 hPa

Accumulation observations are only available from stake measurements, provided by AWI, which were performed weekly for the period April 1992–January 2009. As timing of precipitation is important for correctly simulating the effects of new snow on snow albedo, we combined the stake observations with precipitation predicted by the regional atmospheric climate model RACMO2.3p2 to obtain realistic timing of precipitation in between stake observations, as well as for the post-2009 period. The amount of precipitation modelled by RACMO2 was scaled such that the modelled surface height changes agree with stake measurements; this required a 15.3 % upward adjustment of the modelled precipitation flux.

Figure 3

(a) Seasonal cycles of 2 m temperature (red, left axis), 10 m wind speed (green, right axis) and 2 m specific humidity (blue, right axis). Shaded areas indicate the standard deviations of monthly means. (b) Same as panel (a) for melt (red), net shortwave radiation (blue), net longwave radiation (orange), sensible heat (black), latent heat (magenta) and ground heat (green).

[Figure omitted. See PDF]

Local near-surface climate

Neumayer station is located on an ice shelf $\sim \mathrm{20}$ km from Halvfarryggen ice rise to the south-east, $\sim \mathrm{100}$ km from the ice shelf break (grounding line) to the south, $\sim \mathrm{20}$ km from open water and sea ice to the north and $\sim \mathrm{5}$ km to open water and sea ice to the east. As a result, Neumayer experiences relatively mild conditions without significant impact from katabatic winds but with a pronounced influence of synoptic low-pressure systems passing mainly from west to east in the South Atlantic Ocean to the north of the station. The seasonal cycles of 2 m temperature, 10 m wind and 2 m specific humidity are presented in Fig. a. Summer temperatures around $-\mathrm{4}$ ${}^{\circ }$C and winter temperatures around $-\mathrm{25}$ ${}^{\circ }$C imply a substantial ($>\mathrm{20}$ K) seasonal temperature amplitude based on monthly mean values. This is in line with the formation of a surface-based temperature inversion in winter, a phenomenon that is representative for the flat ice shelves as well as the interior ice domes and in contrast to the topographically steeper escarpment zone, where the quasi-continuous mixing by katabatic flow limits the formation of such an inversion . As expected from the strong link to the air temperature through the Clausius–Clapeyron relation and a high annual mean relative humidity of 82 % (relative to either water or ice, depending on the air temperature), because of the proximity of a saturated snow surface and the ocean, the seasonal cycle of $q_{\mathrm{2}\mathrm{m}}$ closely follows that of temperature.

3 Results: surface energy balance and melt

3.1 SEB model performance and uncertainties

There are several SEB model parameters for which the exact values or formulations are unknown, e.g. the surface roughness length for momentum $z_{\mathrm{0},\mathrm{m}}$, the density of new snow ${\mathit{\rho }}_s$, the stability functions (required to calculate the turbulent scales) and the effective conductivity, which couples the magnitude of $Q_{\mathrm{G}}$ to the temperature gradient in the snow. We estimated the impact of observational and model uncertainties on modelled melt by running the model 600 times while randomly varying all hourly observations within the specified measurement uncertainty ranges (Table ) and using multiple expressions for the heat conductivity and stability functions. Model performance is quantified by comparing modelled with observed $T_{\mathrm{s}}$ and assessing the changes in modelled 24-year cumulative melt. Note that in this section, the albedo based on observations is used to obtain ${\mathrm{SW}}_{\mathrm{net}}$.

The choice of expressions for the stability functions and heat conductivity did not significantly impact the modelled amount of melt (total within 30 mm w.e. or 2.7 %, not shown). The model outcome is more sensitive to the choice of surface roughness length for momentum $z_{\mathrm{0},\mathrm{m}}$ and the lower limit of density of new snow ${\mathit{\rho }}_{\mathrm{s},\mathrm{0}}$: when $z_{\mathrm{0},\mathrm{m}}$ is varied between 0.5 and 50 mm and ${\mathit{\rho }}_{\mathrm{s},\mathrm{0}}$ between 150 and 500 kg m${}^{-\mathrm{3}}$, the cumulative amount of surface melt over the 24-year period varies between 900 and 1220 mm w.e., with higher melt values for smaller values of $z_{\mathrm{0},\mathrm{m}}$ and ${\mathit{\rho }}_{\mathrm{s},\mathrm{0}}$. Optimal values in terms of simulated $T_{\mathrm{s}}$ are $z_{\mathrm{0},\mathrm{m}}=\mathrm{1.65}\mathrm{mm}$ and ${\mathit{\rho }}_{\mathrm{s},\mathrm{0}}=\mathrm{280}$ kg m${}^{-\mathrm{3}}$, resulting in a $T_{\mathrm{s}}$ bias of 0.01 K and an RMSD of 0.79 K (Fig. ). We use these values in the remainder of this study. Figure a and b show modelled 24-year cumulative melt and annual melt (March–February) at Neumayer, combined with uncertainties associated with model parameters. The annual values for year $X$ are obtained by summing monthly values for March of year $X$ until February of year $X+\mathrm{1}$. The total melt amounts to 1154 mm w.e., with a small uncertainty associated with measurement uncertainties ($\mathrm{1}\mathit{\sigma }\approx \mathrm{3}\mathrm{mm}\mathrm{w}.\mathrm{e}.$, i.e. 0.3 %). The method adopted to estimate this uncertainty has its limitations, as measurement errors are probably autocorrelated: if a measurement at one time is disturbed in some way, it is probably disturbed in a similar way at the next time step. Therefore, this result could be interpreted as a lower bound of the uncertainty range, which is supported by the larger uncertainty estimates ($\sim \mathrm{15}$ %) by , who applied a constant systematic error which can be interpreted as an upper bound on the modelled uncertainty range. This also explains why the model outcome is much more sensitive to different values of $z_{\mathrm{0},\mathrm{m}}$, as these runs effectively introduce a systematic error between the true (unknown) value and the chosen value. Furthermore, this approach assumes the true value to be constant, which likely is an oversimplification .

Figure 4

Daily values of modelled vs. measured $T_{\mathrm{s}}$ for the parameter settings used in this study: $z_{\mathrm{0},\mathrm{m}}=\mathrm{1.65}\mathrm{mm}$, ${\mathit{\rho }}_{\mathrm{s},\mathrm{0}}=\mathrm{280}$ kg m${}^{-\mathrm{3}}$.

[Figure omitted. See PDF]

Figure 5

Effect of model uncertainties on (a) cumulative melt and (b) seasonal melt. The shaded area indicates the 1$\mathit{\sigma }$ range due to model uncertainties (changing $z_{\mathrm{0},\mathrm{m}}$ and ${\mathit{\rho }}_{\mathrm{s},\mathrm{0}}$ between their respective values). The vertical grey patches in panel (a) indicate November–February of each season. Note that panel (b) ends earlier than panel (a) because the observations do not cover the 2015–2016 melt season entirely.

[Figure omitted. See PDF]

The sensitivity of modelled cumulative melt to $z_{\mathrm{0},\mathrm{m}}$ is somewhat unexpected. Following Eq. () both $z_{\mathrm{0},\mathrm{h}}$ and $z_{\mathrm{0},\mathrm{q}}$ decrease for increasing $z_{\mathrm{0},\mathrm{m}}$; in combination with the bulk method this acts to dampen the effect of $z_{\mathrm{0},\mathrm{m}}$ on the magnitude of the turbulent fluxes. Our interpretation of this result is that decreasing $z_{\mathrm{0},\mathrm{m}}$ and ${\mathit{\rho }}_{\mathrm{s},\mathrm{0}}$ lead to a decrease in the turbulent fluxes as well as the ground heat flux $Q_{\mathrm{G}}$. This reduces the efficiency with which heat is removed from the surface, in turn allowing more energy to be invested in melt. The obtained value of $z_{\mathrm{0},\mathrm{m}}=\mathrm{1.65}\mathrm{mm}$ is high compared to the average value of $z_{\mathrm{0},\mathrm{m}}=\mathrm{0.1}\mathrm{mm}$ found during a field campaign at Neumayer in 1982 but it is not uncommon for snow surfaces .

Measured values of $T_{\mathrm{s}}$ in excess of the melting point in Fig.  only occurred in the first six seasons; from 1998–1999 onwards they were removed by additional post-processing. These measurements mainly reflect uncertainties in the adopted unit value of longwave emissivity and in measured $\mathrm{LW}\uparrow$, e.g. from sensor window heating and the fact that the downward-facing radiation sensor also measures longwave radiation emitted by the relatively warm air between the surface and the sensor.

3.2 Surface energy balance

Annual (March–February) mean values of near-surface meteorological quantities and SEB components are presented in Table , with seasonal cycles of SEB components presented in Fig. b. These show that the summertime SEB is dominated by the radiation fluxes; despite the high albedo of the snow surface, ${\mathrm{SW}}_{\mathrm{net}}$ is the dominant heat source for the skin layer, whereas ${\mathrm{LW}}_{\mathrm{net}}$ extracts energy from the surface, most efficiently so in summer, when the surface is heated by the sun. In summer, $Q_{\mathrm{L}}$ becomes a significant source of heat loss in the SEB (sublimation), preventing strong negative $Q_{\mathrm{S}}$ (convection). The seasonal cycle of $Q_{\mathrm{G}}$ is small, indicating a small net transport of heat away from the surface in summer and towards the surface in winter. The net annually integrated amount is less than zero as a result of the refreezing of meltwater, warming the subsurface snow layers.

Table 2

Mean annual values and interannual variability (calculated as standard deviations of monthly means) of meteorological variables and SEB components. For precipitation and melt, total annual values are given.

Variable	Yearly mean	Variability
$T_{\mathrm{2}\mathrm{m}}$ (K)	257.1	0.7
$T_{\mathrm{s}}$ (K)	256.0	0.8
$q_{\mathrm{2}\mathrm{m}}$ (g kg${}^{-\mathrm{1}}$)	1.1	0.1
$V_{\mathrm{10}\mathrm{m}}$ (m s${}^{-\mathrm{1}}$)	8.9	0.6
$p$ (hPa)	981.6	2.0
${\mathrm{SW}}_{\mathrm{net}}$ (W m${}^{-\mathrm{2}}$)	20	2
$\mathrm{SW}\downarrow$ (W m${}^{-\mathrm{2}}$)	127	3
$\mathrm{SW}\uparrow$ (W m${}^{-\mathrm{2}}$)	107	2
${\mathrm{LW}}_{\mathrm{net}}$ (W m${}^{-\mathrm{2}}$)	$-\mathrm{28}$	3
$\mathrm{LW}\downarrow$ (W m${}^{-\mathrm{2}}$)	218	5
$\mathrm{LW}\uparrow$ (W m${}^{-\mathrm{2}}$)	246	4
$Q_{\mathrm{S}}$ (W m${}^{-\mathrm{2}}$)	14.5	2.7
$Q_{\mathrm{L}}$ (W m${}^{-\mathrm{2}}$)	$-\mathrm{6.3}$	1.2
$Q_{\mathrm{G}}$ (W m${}^{-\mathrm{2}}$)	0.7	0.4
$M$ (W m${}^{-\mathrm{2}}$)	0.5	0.4
Precipitation (mm w.e.)	415	86
Melt (mm w.e.)	50	42

Statistically significant and previously unreported trends over the full 24-year period (not shown) are detected in $\mathrm{LW}\uparrow$ ($-\mathrm{0.28}\pm \mathrm{0.14}$ W m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$) and $Q_{\mathrm{S}}$ ($+\mathrm{0.21}\pm \mathrm{0.07}$ W m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$). Both of these are a result of wintertime trends. $\mathrm{LW}\uparrow$ is linked directly to $T_{\mathrm{s}}$, which shows a statistically insignificant negative trend ($-\mathrm{0.029}\pm \mathrm{0.026}$ K yr${}^{-\mathrm{1}}$), which in magnitude exceeds the negative trend in $T_{\mathrm{2}\mathrm{m}}$ ($-\mathrm{0.0045}\pm \mathrm{0.02}$ K yr${}^{-\mathrm{1}}$; assuming a normal distribution, the probability that the negative trend in $T_{\mathrm{s}}$ is greater in magnitude than the trend in $T_{\mathrm{2}\mathrm{m}}$ is 0.76). As a result, the air temperature gradient near the surface has increased, enhancing $Q_{\mathrm{S}}$. The negative trend in $T_{\mathrm{s}}$ originates from a decrease in $\mathrm{LW}\downarrow$ ($-\mathrm{0.26}\pm \mathrm{0.17}$ W m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$), which is in turn driven by a slight decrease in cloud cover ($-\mathrm{0.003}\pm \mathrm{0.001}$ yr${}^{-\mathrm{1}}$). This is suggested independently by the decrease in average winter humidity ($-\mathrm{0.004}\pm \mathrm{0.002}$ g kg${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$). These findings agree with and , who determined from satellite observations that summer cloud cover has decreased over that part of coastal Antarctica in the period 1979–2011.

3.3 Melt season

Melt occurs at Neumayer from November to February (Fig. ) but is highly variable from year to year. The mean annual amount of melt is 50 mm w.e. with an interannual variability of 42 mm w.e. and a range of 2 mm w.e. in 1999–2000 to 176 mm w.e. in 2012–2013. Most melt occurs in December and January and the surface only sporadically reaches melting point in February. Only in 2007 did melt occur in November, and no melt occurs outside these 4 months. The cumulative melt occurring at Neumayer shows stepwise increases (Fig. a), which represent the peaked melt seasons, in which melt occurs on average on $\mathrm{18}\pm \mathrm{10}$ d. The uncertainty in the number of melt days due to the chosen values of $z_{\mathrm{0},\mathrm{m}}$ and ${\mathit{\rho }}_{\mathrm{s},\mathrm{0}}$ is relatively small compared to the interannual variability in melt totals (Fig. ), implying that this choice does not significantly affect the modelled melt duration, but it does affect the total melt.

Figure 6

Average number of melt days per month at Neumayer. The inner error bars (with larger caps) indicate the 1$\mathit{\sigma }$ uncertainty range resulting from the runs performed with different settings for roughness length $z_{\mathrm{0}}$ and lower limit of new snow density ${\mathit{\rho }}_{\mathrm{s},\mathrm{0}}$ (Sect. ). The outer error bars (with smaller caps) indicate the 1$\mathit{\sigma }$ range of the interannual variability.

[Figure omitted. See PDF]

To investigate the link between melt and climate, we compare the two summers with the highest (2003–2004 and 2012–2013, on average 145 mm w.e.) and lowest (1999–2000 and 2014–2015, on average 4 mm w.e.) melt amounts. Figure  shows the meteorological and SEB components for these years, averaged over December and January. The largest differences are found in $T_{\mathrm{2}\mathrm{m}}$ ($+\mathrm{2.3}$ K) and ${\mathrm{SW}}_{\mathrm{net}}$ ($+\mathrm{17}$ W m${}^{-\mathrm{2}}$); based on the measurement uncertainties (Table ), these differences are significant. In cold summers, the low $T_{\mathrm{2}\mathrm{m}}$ corresponds to a stronger temperature inversion ($T_{\mathrm{2}\mathrm{m}}-T_{\mathrm{s}}$), more longwave cooling, less sublimation and a larger $Q_{\mathrm{S}}$. $\mathrm{SW}\downarrow$ and $\mathrm{LW}\downarrow$ show almost no difference between high and low melt seasons; therefore, the difference in ${\mathrm{SW}}_{\mathrm{net}}$ cannot be caused by a change in cloud cover and is likely caused solely by surface albedo, which suggests an important role for the SMAF. This will be elaborated upon in the next section. Finally, the direction of $Q_{\mathrm{G}}$ is reversed: in high melt years, the surface is warmed from below, while in low melt years the surface loses heat to the subsurface. More refreezing of meltwater in high melt years warms the near surface snow layers, which in turn leads to a conductive heat flux towards the surface.

Figure 7

Average values of some SEB components (a, b) and some meteorological variables (c) for December and January in the years with the highest (2003–2004 and 2012–2013, in light grey) and lowest (1999–2000 and 2014–2015, in dark grey) amount of melt, as identified in Sect. . Note that $\mathrm{SW}\downarrow$, $\mathrm{SW}\uparrow$, $\mathrm{LW}\downarrow$ and $\mathrm{LW}\uparrow$ are scaled by a factor of 10 in panel (a) for clarification.

[Figure omitted. See PDF]

Using the subsurface radiation model of , the influence of subsurface penetration of shortwave radiation is estimated. Its inclusion increases the modelled cumulative amount of melt by 13 %, from 1154 to 1326 mm w.e. The absorbed shortwave radiation heats the subsurface layers, but the heat cannot be transported away as effectively as would happen at the surface by turbulent fluxes and longwave radiation. This leads to an increase in total melt.

The findings presented in this section are in good agreement with , who used a similar approach to calculate the SEB at Neumayer but used a lower value for $z_{\mathrm{0},\mathrm{m}}=\mathrm{0.32}\mathrm{mm}$ and a higher snow density that was assumed constant with depth (420 kg m${}^{-\mathrm{3}}$ in their study vs. 280 kg m${}^{-\mathrm{3}}$ in this study). Compared to melt estimates from the Larsen C ice shelf, obtained through a similar modelling approach by , melt at Neumayer is weak. Owing to its more northerly location, on the Larsen C ice shelf an annual (2009–2011) average melt energy of 2.8 W m${}^{-\mathrm{2}}$ is obtained, compared to the 2009–2011 annual average of 0.7 W m${}^{-\mathrm{2}}$ obtained at Neumayer. Furthermore, in November and February melt occurs much more frequently on the Larsen C ice shelf.

4 Results: the snowmelt–albedo feedback

The SMAF is a well-known phenomenon but has not before been quantified for Antarctica. The feedback occurs after the rapid growth of snow grains when meltwater penetrates into the subsurface and refreezes. Because a photon travels farther through snow with large particles than in new snow with smaller particles on average, the probability of it being absorbed is increased, effectively lowering the surface albedo . Even without melt, albedo decreases when snow ages, following grain growth from dry snow metamorphism, but this is a much slower process which mainly depends on temperature gradients in the snow, favouring moisture transport onto larger grains. Precipitation of new, fine-grained snow has been shown to inhibit the albedo decrease by metamorphism on the Antarctic plateau .

To quantify the SMAF at Neumayer, we need to be able to switch on and off the albedo dependency on melt-driven grain growth. To that end, we implemented an albedo parameterisation in the SEB model, as described in Sect. . Because no data on grain size are available from Neumayer, we optimised the albedo model performance by maximising the correspondence between (1) modelled and observed hourly $\mathrm{SW}\uparrow$ and (2) the total melt obtained from the calculations based on observed albedo (Sect. ). We compare $\mathrm{SW}\uparrow$ instead of the albedo itself because by doing so the hourly values are naturally weighted with its contribution to ${\mathrm{SW}}_{\mathrm{net}}$ and hence its importance for the SEB. We then perform several runs with different processes switched on and off affecting the surface albedo to investigate the importance of the SMAF for melt at Neumayer (Sect. ).

4.1 Optimising the albedo parameterisation

The albedo parameterisation, and especially the expression for snow grain size (Eq. ), contains several parameters that are not well constrained, such as new snow grain size $r_{\mathrm{e},\mathrm{0}}$ and refrozen snow grain size $r_{\mathrm{e},\mathrm{r}}$. These parameters were varied within reasonable ranges to optimise the results: new snow grain sizes between 0.04 and 0.3 mm and refrozen snow grain sizes between 0.1 and 10 mm. The best comparison with observed albedo was achieved when using the look-up table for dry snow metamorphism, $\mathrm{d}{r}_{\mathrm{e},\mathrm{dry}}$, corresponding to a grain size of $\mathrm{0.055}\mathrm{mm}$.

The first step in optimising the parameterisation was to split the summer season into two parts, the “dry” and the “wet” season. The respective starts of the dry and wet seasons are the first day on which the sun rises more than 15${}^{\circ }$ above the horizon and the first day that surface melt occurs. The wet season ends when the sun no longer rises higher than 15${}^{\circ }$. For the dry season, we varied the dry snow metamorphism factor and the new snow grain size to best match observed $\mathrm{SW}\uparrow$. This resulted in a new snow grain size of $\mathrm{0.25}\mathrm{mm}$. This value is then used in the second step, in which the refrozen snow grain size $r_{\mathrm{e},\mathrm{r}}$ is varied to best match the modelled cumulative melt using observed albedo. This was achieved for a refrozen snow grain size of 1.45 mm.

This value for refrozen snow grain size is compatible with the typical largest grains in dry metamorphosed snow of O(1 mm) and which used as a lower limit for refrozen snow grains. and present observations of snow grain sizes on the Antarctic plateau during field campaigns in 2012–2013 and 2013–2014 as well as estimates from satellite observations. On the plateau, summer temperatures are comparable to Neumayer winter temperatures. report summertime snow grain size estimates of approximately 0.11 mm (Fig. 6 in their study, reported as a specific surface area $\mathrm{SSA}=\frac{\mathrm{3}}{{\mathit{\rho }}_{\mathrm{i}}{r}_{\mathrm{e}}}$, where ${\mathit{\rho }}_{\mathrm{i}}$ is the density of ice and $r_{\mathrm{e}}$ is the snow grain size). In our study, wintertime snow grain sizes approach 0.21 mm. The difference is expected as the plateau is generally much colder than Neumayer. The seasonal cycle of modelled average specific surface area in the upper 7 cm (Fig. ) is comparable to the one presented in , although the wintertime values are probably too low. For the purpose of this study, however, the accurate representation of surface albedo during winter is less relevant as there is no shortwave radiation in winter.

Figure 8

Seasonal cycle of modelled average grain size in the upper 7 cm for the period 2000–2014. The grain size is expressed in terms of specific surface area ($\mathrm{SSA}=\frac{\mathrm{3}}{{\mathit{\rho }}_{\mathrm{i}}{r}_{\mathrm{e}}}$) rather than grain size itself to allow for a comparison with Fig. 6 of . The vertical grey patches indicate November–February of each season.

[Figure omitted. See PDF]

When the adopted albedo values are combined with the observations of $\mathrm{SW}\uparrow$, the model adequately reproduces the incoming shortwave radiation (Fig. , $\mathrm{bias}=+\mathrm{0.93}$ W m${}^{-\mathrm{2}}$, $\mathrm{RMSD}=\mathrm{7.3}$ W m${}^{-\mathrm{2}}$), providing confidence in the modelled albedo.

Figure 9

Measured vs. modelled daily average incoming shortwave radiation ($\mathrm{SW}\downarrow$). The modelled $\mathrm{SW}\downarrow$ is obtained by dividing the hourly measured $\mathrm{SW}\uparrow$ by the hourly modelled albedo.

[Figure omitted. See PDF]

4.2 Magnitude of the snowmelt–albedo feedback

Three experiments with the SEB model were carried out in addition to the original run ($R_{\mathrm{0}}$), which uses the measured albedo:

• $R_{\mathrm{1}}$: the average measured albedo (0.84, determined by adding all $\mathrm{SW}\downarrow$ and $\mathrm{SW}\uparrow$ for all measurements when the sun is higher than 15${}^{\circ }$ above the horizon and taking the ratio between the two) is prescribed for the entire period.

• $R_{\mathrm{2}}$: the full albedo parameterisation is used.

• $R_{\mathrm{3}}$: refrozen snow does not contribute to the changing snow characteristics, i.e. $f_{\mathrm{r}}=\mathrm{0}$ in Eq. ().

Figure a and b show time series of modelled cumulative and seasonal surface melt for the four experiments. Experiment $R_{\mathrm{1}}$ underpredicts melt in most seasons, yielding a mean annual amount of surface melt of $\mathrm{39}\pm \mathrm{27}$ mm w.e. yr${}^{-\mathrm{1}}$ (compared to $\mathrm{50}\pm \mathrm{42}$ mm w.e. yr${}^{-\mathrm{1}}$ for experiment $R_{\mathrm{0}}$). More melt was modelled in the 1995–1996 melt season, which was characterised by frequent precipitation events and cloudy conditions, keeping observed albedo higher than the long-term mean. Because the albedo parameterisation (used in experiment $R_{\mathrm{2}}$) has been calibrated to match observed albedo, experiment $R_{\mathrm{2}}$ adequately reproduces the amount of seasonal melt ($\mathrm{50}\pm \mathrm{34}$ mm w.e. yr${}^{-\mathrm{1}}$), although melt, e.g. in the 2012 melt season, is underestimated. Run $R_{\mathrm{3}}$ represents the situation in which the SMAF has been switched off, leading to significantly underpredicted melt ($\mathrm{21}\pm \mathrm{16}$ mm w.e. yr${}^{-\mathrm{1}}$).

Figure 10

(a) Time series of the modelled cumulative amount of melt for the run with measured albedo ($R_{\mathrm{0}}$, blue), a constant albedo of 0.84 ($R_{\mathrm{1}}$, red), a run in which refrozen snow impacts snow grain size ($R_{\mathrm{2}}$, yellow) and a run in which snow grain size is not influenced by refrozen snow ($R_{\mathrm{3}}$, purple). (b) Same as panel (a) but for seasonal amount of melt. (c) Ratio of modelled surface melt between yellow and purple lines in panels (a) and (b) (runs $R_{\mathrm{2}}$ and $R_{\mathrm{3}}$ respectively). The grey area indicates the uncertainty coming from the uncertainty in the determination of $\mathit{\tau }$ (Fig. b), $\pm \mathrm{5}$ W m${}^{-\mathrm{2}}$ measurement uncertainty in $\mathrm{SW}\uparrow$ and the inclusion of shortwave radiation penetration.

[Figure omitted. See PDF]

Defining the strength of the SMAF as the ratio between the total seasonal surface melt in experiments $R_{\mathrm{2}}$ and $R_{\mathrm{3}}$, we obtain an average value of 2.6, with a range of 1.3 (1996–1997) to 4.8 (1993–1994; see Fig. c). The effect of subsurface penetration of shortwave radiation on this result is estimated by repeating the above experiments with an inclusion of the radiation penetration model of . This yielded an average SMAF of 2.3, ranging from 1.5 (2005–2006) to 3.2 (2002–2003). The main difference between the two experiments is the reduced interannual variability: including penetration of shortwave radiation does not yield SMAF values larger than 3.5. Shortwave radiation penetration heats the subsurface, causing subsurface melt which is less affected by the SMAF because the radiative flux is smaller in the subsurface. Therefore, the “extreme” years in the sense of SMAF are less distinct in the experiment with shortwave radiation penetration. The effect of shortwave radiation penetration is included in the uncertainties indicated in Fig. c. Combining this with the uncertainties in observed $\mathrm{SW}\uparrow$ and the determination of $\mathit{\tau }$ (Fig. b) leads to uncertainties in the determination of the SMAF of typically 15 %, with a range of 4 % (1995–1996) to 32 % (1993–1994).

A weak positive correlation was found between SMAF and $\mathrm{SW}\downarrow$ ($R^{\mathrm{2}}=\mathrm{0.15}$, $p=\mathrm{0.07}$): if $\mathrm{SW}\downarrow$ increases, more energy is available at the surface for melting, which is then in turn further intensified by SMAF. Another weak negative correlation was found between SMAF and summer precipitation ($R^{\mathrm{2}}=\mathrm{0.13}$, $p=\mathrm{0.1}$): snowfall inhibits SMAF as it effectively “resets” the surface albedo as was also shown by in a dry region.

Only few studies report on the SMAF concerning the darkening of snow rather than disappearance of it. provide relationships between anomalies of seasonal $T_{\mathrm{2}\mathrm{m}}$ and ${\mathrm{SW}}_{\mathrm{net}}$ (Figs. 5 and 12 of ). They find a negative relationship for accumulation regions, i.e. lower 2 m temperatures are associated with smaller ${\mathrm{SW}}_{\mathrm{net}}$. No such relationship is found for Neumayer (not shown).

5 Conclusions

In this study, we used 24 years of high-quality meteorological and radiation observations from the BSRN station Neumayer, situated on the Ekström ice shelf, East Antarctica, to force a surface energy balance model. The primary goal was to calculate the amount of melt at Neumayer and to investigate the importance of the snowmelt–albedo feedback (SMAF). Model performance was evaluated based on the difference between modelled and measured surface temperature, and the modelled melt was tested for measurement and model parameter uncertainties. We found that measurement uncertainties, when considered random in time, do not significantly impact modelled melt at Neumayer over the full 24-year period ($<\mathrm{0.5}$ % difference). However, melt amount and model performance are sensitive to the values chosen for the surface roughness length for momentum $z_{\mathrm{0},\mathrm{m}}$ and lower limit of new snow density ${\mathit{\rho }}_{\mathrm{s},\mathrm{0}}$; thus accurate measurements of these values would further improve future modelling studies. Our results confirm that melt at Neumayer is an intermittent process, occurring on average on only 18 d each summer, totalling 50 mm w.e. and with an interannual variability of 42 mm w.e. Melt occurs mainly in December and January, sporadically in February and only once melt was modelled in November. Significant and previously unreported trends were found in the net longwave radiation (decreasing) and the sensible heat flux (increasing), but these are unrelated to the melt at Neumayer as they mainly occur in winter and are attributed to a decrease in cloud cover.

The main difference between high and low melt years was found to be surface albedo, implying an important role for the SMAF. We quantified SMAF by implementing and tuning an albedo parameterisation in the SEB model, which includes the effects of snowfall and wet and dry snow metamorphism on albedo. The albedo parameterisation adequately reproduces the seasonal variability in snow grain size, compared to measurements on the Antarctic Plateau . Our derived wintertime snow grain sizes at Neumayer are somewhat smaller than the satellite-derived summertime snow grain sizes at the Antarctic Plateau owing to the lower temperatures on the plateau. Our main finding is that SMAF on average enhances surface melt at Neumayer by a factor of $\mathrm{2.6}\pm \mathrm{0.8}$.

Weak correlations were found of SMAF with summertime $\mathrm{SW}\downarrow$ and precipitation ($\mathrm{0.1}<R^{\mathrm{2}}<\mathrm{0.2}$). To assess how the importance of the snowmelt–albedo feedback varies spatially and temporally, the next step in this research will be applying this method to other sites in Antarctica and a regional climate model .

Code and data availability

The Neumayer data are available upon request via the website of AWI (https://bsrn.awi.de/data/data-retrieval-via-pangaea/, last access: 9 June 2016). The model output is available upon request by the authors.

Author contributions

CLJ performed the study and wrote the manuscript. PKM assisted with the implementation of the albedo parameterisation. GKL was in charge of the Neumayer data. CHR, PKM, GKL and MRvdB have commented on the manuscript.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

We would like to thank the AWI for maintaining the station and the Baseline Surface Radiation Network (BSRN) for providing the data, with special thanks to Amelie Driemel for creating a citation reference and Holger Schmithüsen for helping us to interpret the data. Michiel R. van den Broeke acknowledges support from the Netherlands Earth System Science Centre (NESSC). We would like to thank Ghislain Picard and Achim Heilig for their constructive comments.

Financial support

This research has been supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (grant no. 866.15.204).

Review statement

This paper was edited by Mark Flanner and reviewed by Ghislain Picard and Achim Heilig.