Estimating Greenland surface melt is hampered by melt induced dampening of temperature variability

Abstract

The positive degree-day (PDD) model provides a particularly simple approach to estimate surface melt from land ice based solely on air temperature. Here, we use a climate and snow pack simulation of the Greenland ice sheet (Modèle Atmosphérique Régional, MAR) as a reference, to analyze this scheme in three realizations that incorporate the sub-monthly temperature variability differently: (i) by local values, (ii) by local values that systematically overestimate the dampened variability associated with intense melting or (iii) by one constant value. Local calibrations reveal that incorporating local temperature variability, particularly resolving the dampened variability of melt areas, renders model parameters more temperature-dependent. This indicates that the negative feedback between surface melt and temperature variability introduces a non-linearity into the temperature – melt relation. To assess the skill of the individual realizations, we hindcast melt rates from MAR temperatures for each realization. For this purpose, we globally calibrate Greenland-wide, constant parameters. Realization (i) exhibits shortcomings in the spatial representation of surface melt unless temperature-dependent instead of constant parameters are calibrated. The other realizations perform comparatively well with constant parametrizations. The skill of the PDD model primarily depends, however, on the consistent calibration rather than on the specific representation of variability.

1. INTRODUCTION

The Greenland ice sheet (GrIS) is sensitive to a warming climate, as surface melting is a major element of its mass budget (Rignot and others, 2011). Ice volume reconstructions reveal that past climate changes affected the GrIS massively: Greenland's ice volume experienced changes of ~5–9 m sea-level equivalent (SLE) since the last interglacial's peak sea level (Vasskog and others, 2015). The total present-day volume of the GrIS, in comparison, is of the same order of magnitude (~7.4 m SLE).

To interpret historic and recent mass changes of the GrIS, it is necessary to quantify the components of the surface mass balance (SMB), which primarily are surface melt, accumulation and refreezing. Direct measurements of surface melt rates are only available at isolated locations and until recently were not designed to cover more than a few years. Over the last decade, however, a network of automatic mass-balance stations has been established to monitor the annual mass loss of the GrIS (Ahlstrom and others, 2008). Since 2002, the Gravity Recovery and Climate Experiment (GRACE) (Tapley and others, 2004) satellite mission allows to detect integral mass changes of the GrIS (Wouters and others, 2014), but additional data are necessary to separate and quantify the individual SMB components such as surface melt (Sasgen and others, 2012; Tedesco and Fettweis, 2012). In particular, high-resolution regional climate models such as the Regional Atmospheric Climate Model (RACMO) (Noel and others, 2015) or the Modèle Atmosphérique Régional (MAR) (Fettweis and others, 2017) allow for a separation of the SMB.

Where high-resolution modelling is not feasible, a common approach is to estimate surface melt rates, M with realizations of the positive degree-day (PDD) model formulated for example in Reeh (1989). The term PDD refers to the temporal integral of near surface temperatures T exceeding the melting point, sampled at daily or higher frequency. The PDD scheme linearly relates PDDs to snow and ice melt through respective empiric degree-day factors (DDF) in the form M = DDF × PDD. Following Braithwaite (1985), PDD values are commonly approximated from monthly mean temperatures (or, from annual mean and mean July temperatures based on a sinusoidal seasonal cycle, Reeh (1989)), assuming that on sub-monthly timescales temperatures are normally distributed around the monthly mean temperatures with a standard deviation, σ. One benefit of this approach is its simplicity, as it only requires information about near surface air temperature. Even spatially incomprehensive or coarse resolution temperature data still permit some rough estimate of the surface melt rates, as temperatures can be downscaled via lapse rate corrections. The PDD model has been thus applied to estimate surface melt rates from coarse resolution climate simulations (e.g. Charbit and others, 2013; Heinemann and others, 2014; Roche and others, 2014; Ziemen and others, 2014; Gierz and others, 2015) or, extending the instrumental record into the past, from temperature reconstructions deduced from ice cores (e.g. Box, 2013; Wilton and others, 2017). In many such applications σ is assumed to be constant in time and space and then σ commonly is chosen in the range 4.5–5.2°C. Fausto and others (2009) pointed out that the ablation zone of the GrIS exhibits a considerably smaller variability during the melting season with σ _T ranging between 1.7 and 2.9°C. Consequently, a growing number of studies approximate PDDs using spatially or temporally variable σ (e.g. Roche and others, 2014; Ziemen and others, 2014; Contoux and others, 2015; Wilton and others, 2017).

In modelling studies, the DDFs of snow and ice can be treated as convenient tuning parameters, which can be constrained by observed mass changes from satellite data. Also, many applications of the PDD model in Earth system modelling adopt DDFs of 3 and 8 mm °C⁻¹ d⁻¹ for snow and ice, respectively, following an approach, which was applied in a simulation of the GrIS dynamics by Huybrechts and others (1991). Direct calibrations of DDFs, however, rely on the few available measurements of melt rates which often differ in sampling frequency and period. In a global compilation, Hock (2003) reviews DDF calibrations from ice sheets and glaciers, which exhibits a wide range of values for different locations indicating that DDFs strongly depend on parameters such as altitude, latitude or local conditions (e.g. debris cover, liquid water content, shading). For Greenland DDF estimates for ice range from 5.8 to 12 mm °C⁻¹ d⁻¹ for stations at altitudes lower than 1000 m a.s.l., while measurements at higher elevations yield values of up to 20 mm °C⁻¹ d⁻¹. According to a calibration to long-term measurements and energy-balance modelling from southern Greenland (Johannesson and others, 1995) most applications employ for snow DDF ~3 mm °C⁻¹ d⁻¹ for Northern Hemispheric ice sheets. However, a calibration using data from two stations at 1020 m and 1520 m a.s.l. in Southwest Greenland suggests DDF of snow well above 20 mm °C⁻¹ d⁻¹ (van den Broeke and others, 2010). In consequence, recent studies have proposed modifications of the PDD method, which make use of temperature-dependent DDFs (Tarasov and Peltier, 2002; Greve, 2005; Fausto and others, 2009).

Charbit and others (2013) find that ice sheet simulations of the last glacial period exhibit a strong sensitivity to the realization of the PDD model. Specifically Charbit and others (2013) compare schemes with constant (Reeh, 1989) or temperature dependant DDFs (Tarasov and Peltier, 2002; Fausto and others, 2009), using either PDD approximations based upon spatially variable (Fausto and others, 2009) or constant standard deviation (Reeh, 1989; Tarasov and Peltier, 2002). Rogozhina and Rau (2014) evaluate a similar selection of PDD model realizations with respect to their ability to reproduce Greenland's surface melt rates of the years 1958–2001 as simulated from the regional climate – snow pack model RACMO in combination with GRACE ice mass measurements. Rogozhina and Rau (2014) find a good agreement with the surface melt rates of the regional model if spatially variable standard deviations are combined with the parameterization of Fausto and others (2009). The other combinations exhibit a strong, but mostly temporally constant bias.

In this study, we aim to fathom principles behind heterogenous DDF calibrations. Direct observations are still spatially sparse, rarely cover more than a decade and may be influenced by local conditions which are not representative on a larger scale. For a comprehensive analysis we thus use a simulation of Greenlands climate of the years 1948–2016 with the state-of-the-art regional climate and snow pack model MAR (Fettweis and others, 2017) as a reference. We consider and evaluate three PDD model realizations, which differ in the way that temperature variability is represented, and calibrate these individually. We optimize the calibration on the basis of the first 25 years of the reference simulation. We first investigate the spatial properties of DDFs if locally calibrated. We then determine ideal Greenland-wide parameterization for each formulation and test their skill to reproduce MAR surface melt rates in their spatial characteristics and evolution, using the MAR temperatures as input.

2. METHODS

2.1. The modèle atmosphérique régional

We use outputs from version 3.5 of the MAR regional climate model applied over the GrIS for calibration and evaluation of different PDD-model realizations. The MAR, is a coupled land-atmosphere regional climate model featuring an atmospheric model (Fettweis, 2007; Reijmer and others, 2012) and including a snow pack component. The simulation considered here is analyzed in detail in Fettweis and others (2017). The simulation covers Greenland and surrounding areas at a 20 km spatial resolution, forced at the lateral boundaries and ocean surface with reanalysis data from the National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP) for the years 1948–2016 (Kalnay and others, 1996). MAR has been validated against satellite and in situ data (e.g. Lefebre and others, 2005; Fettweis and others, 2011) and has been used in multiple studies to simulate ice sheet and atmospheric parameters, such as SMB, surface melt and atmospheric circulation near Greenland (e.g. Fettweis and others, 2011; Franco and others, 2013; Tedesco and others, 2013).

2.2. The PDD-model

A PDD model (Braithwaite, 1985; Reeh, 1989) allows estimating the monthly surface melt rates from temperatures T sampled at daily or higher frequency. The PDD is here defined to be the cumulative temperature exceeding melting point T ⁺ (with T ⁺ = max(T, T ₀), T ₀ = 0°C) throughout a time interval [t ₀, t ₁]

(1)$$ PDD = \int_{t0}^{t1} {T^ +} (t) dt. $$

In practice, temperatures will be sampled or simulated with discrete time steps Δt and are often only available on monthly or longer timescales. Therefore, we focus on realizations of the PDD model that approximate PDD based on monthly statistics following an approach from Braithwaite (1985): assuming that temperatures are normally distributed around monthly or climatologic mean temperatures T̄, the PDD of a given month can be approximated as

(2)$$ {PDD}\left({\bar T},\sigma\right) = {\displaystyle{\Delta t} \over {\sigma\sqrt{2\pi}}} \int_0^\infty {Texp\left(-{\displaystyle{(T - {\bar T})^{2}} \over {2\sigma^{2}}}\right) dT}, $$

with Δt = 1 month.

Here, we consider surface melt estimates of the form

(3)$$ M = {\rm DDF} \times {PDD}({\bar T}, \sigma), $$

and distinguish between snow and ice melt by applying distinct constant DDF for snow and ice, respectively. The surface melt of a given month is considered in two phases: any surface layer of snow must be melted, before the subjacent ice is melted with the remaining PDD that is not used up by snow melt.

We consider different options to represent unresolved temperature variability by σ. We focus on the slush zone below the runoff line, which is here defined to be the region where mean surface melt rates M exceed refreezing rates. In Greenland the runoff line is usually positioned above the equilibrium line where ablation and accumulation balance.

In the first instance, we consider two approximations, PDD_VAR = PDD(T̄, σ_VAR(x, t)) and PDD_CONST = PDD(T̄, σ_CONST = 5°C) and compare these with PDD_6h, which is directly calculated from 6-hourly temperatures according to (1). Here, T̄ are monthly mean temperatures from the MAR simulation and σ _VAR represents its monthly temperature variability calculated at each grid point by combining the monthly mean of the daily temperature amplitude dT _day = 0.5 (T̄ _max − T̄ _min) and the standard deviation of the daily mean temperatures throughout a month, σ_month, according to:

(4)$$ \sigma_{\rm VAR} = \sqrt{(\sigma_{\rm month}^2 + (0.564 \; dT_{\rm day})^2)}. $$

The factor 0.564 in this equation is somewhat smaller than what would be assumed for a sinusoidal daily cycle and was determined from the MAR calibration period by optimizing the regression between PDD_VAR = PDD(T̄, σ_VAR) and the PDD_6h from 6-hourly temperature data. Additionally, we consider PDD_EFF = PDD(T̄, σ_EFF) with

(5)$$ \sigma_{\rm EFF} = \sqrt{(\sigma_{{\rm month}}^2 + (0.564 \; {\max}[dT_{\rm day},{5}^{\circ} {\rm C}])^{2})} $$

assuming that daily temperature amplitude never falls below 5°C. This approximation aims to eliminate a negative feedback that may influence the temperature–melt relation if mean temperatures approach melting point: in this case, not only do temperature variations induce surface melt, but, on the other hand, surface melt (or refreezing) dampens temperature variations. This feedback is particularly prominent in the diurnal variability. If surface melt occurs, the maximum (daytime) temperature will be lowered due to latent heat flux and likewise minimum (nighttime) temperatures will be increased due to refreezing meltwater. In the MAR simulation we find, that above melting surfaces daily temperature amplitude declines from ~5°C–1.5°C when monthly mean temperatures approach melting point from both directions (Fig. 1). σ _EFF can be understood as a rough estimate of an ‘effective’ temperature variability, freed from the dampening effect of phase transitions on the daily temperature amplitude. Unlike parameterizations that aim to explicitly resolve the low variability in the ablation zone (Fausto and others, 2011; Seguinot and Rogozhina, 2014; Wake and Marshall, 2015), this parameterization raises the low variability in the ablation zone to suppress the effect of negative feedbacks.

Fig. 1. Multiyear monthly mean daily temperature amplitude averaged over the first 25 years of the MAR simulation (calibration period) against respective multiyear monthly mean temperatures. The size of each point is scaled by the surface melt rate, colours reflect standard deviation of daily mean temperatures throughout July; grid points with melt rates M < 2 mm d⁻¹ are not represented.

Parameterizing variability by =5°C, as applied in PDD_CONST, strongly overestimate variability. In summer, σ _VAR rather takes values ~2°C–3.5°C below the runoff line (Fig. 2, left). Applying a lower bound on daily temperature amplitudes (PDD_EFF) essentially increases σ _EFF by 1–1.6°Cd⁻¹ below the runoff-line compared with σ _VAR (Fig. 2, right). At high elevation or outside of the melt season σ _EFF and σ _VAR are almost identical .

Fig. 2. July σ _VAR (left) and σ _EFF (right) averaged over the first 25 years of the MAR simulation (calibration period). Black contour indicates the runoff line.

In consequence, PDD_VAR agrees well with PDD_6h, whereas PDD_CONST reveals a strong bias relative to PDD_6h which is particularly pronounced for small values and dominantly positive for PDD_6h > 2°C d⁻¹. The bias in PDD_EFF is less pronounced but also in the main positive (Fig. 3).

Fig. 3. PDD_VAR (a), PDD_CONST (b) and PDD_EFF (c) as functions of PDD_6h. Identity is displayed as a grey line in all panels for comparison.

2.3. Calibration of the PDD model

In this study, we first calibrate DDFs for each grid point individually before a Greenland-wide calibration is considered in later sections. We calibrate DDFs for snow and ice and for each PDD approximation (DDF_{VAR, snow}, DDF_{CONST, snow}, DDF_{EFF, snow}, DDF_{VAR, ice}, DDF_{CONST, ice}, DDF_{EFF, ice}). For calibration, we use monthly means of the different approximations of PDD together with the local snow height (SNH), snow fall (SF) and melt rates, M _MAR from the MAR simulation. We consider a 25-year calibration period a suitable choice. Relative to parameters obtained from a 69-year calibration period (i.e. the length of the simulation), errors of <10, 5, 3 and 2% are obtained for periods longer than 3, 6, 16 and 54 years (not shown). Thus, a 25-year period offers a reasonable accuracy, while the remaining 44 years allow to assess the temporal stability of the individual realizations. We determine snow melt M _{MAR, snow} as

(6)$$ M_{\rm MAR, snow} = {\min} \left({\displaystyle{\rm SNH} \over {dt}} + {\rm SF},M_{\rm MAR}\right),\quad dt = 1 \, {\rm month} $$

and ice melt as

(7)$$ M_{\rm MAR, ice} = M_{\rm MAR} - M_{\rm MAR,\rm snow}. $$

To determine DDF of snow locally at a grid point i, we identify those months, ms, in which snow melt is predominant (M _{MAR, ice} ≈ 0, PDD_,snow ≈ PDD) and use these to obtain a local calibration,

(8)$$ {\rm DDF}_{,{\rm snow}}(i) = {\displaystyle{\sum \nolimits_{ms} M_{\rm MAR}(i,ms)} \over {\sum \nolimits_{ms} {\rm PDD}(i,ms)}}, $$

for all glaciated grid points i. As the calibration becomes inaccurate where substantial surface melting does not occur on a regular basis, we reject the local calibration as meaningless where PDDs accumulated over the 25-year period do not exceed 10°C d⁻¹. For a Greenland-wide calibration we modify this method as

(9)$$ {\rm DDF}_{,{\rm snow}} = {\displaystyle{\sum \nolimits_{i,ms} M_{\rm MAR}(i,ms)} \over {\sum \nolimits_{i,ms} {\rm PDD}(i,ms)}}. $$

In a second step, we can consider all months m and all grid points i to calibrate DDF of ice. Where melt rates are large enough to remove the snow cover and melt subjacent ice, we determine the proportion of PDD, which is available for ice melt as

(10)$$ {\rm PDD}_{,{\rm ice}}(i,m) = {\rm PDD}(i,m)-{\rm PDD}_{,{\rm snow}} $$

with

(11)$$ {\rm PDD}_{,{\rm snow}}(i,m) = {\displaystyle{M_{\rm MAR, snow}}(i,m) \over {\rm DDF}_{,{\rm snow}}(i)}, $$

DDF of ice is then analogously to DDF_,snow calibrated locally as

(12)$$ {\rm DDF}_{,{\rm ice}}(i) = {{\displaystyle{\sum \nolimits_{m} M_{\rm MAR, ice}}(i,m)} \over {\sum \nolimits_{m} {\rm PDD}_{,{\rm ice}}(i,m)}}, $$

for all grid points that meet the condition that PDD_,ice integrated over the 25-year period exceeds 10°Cd⁻¹. Analogously, for the Greenland-wide calibration this method is modified as in Eqn (9).

3. RESULTS

3.1. Influence of temperature variability on locally-calibrated parameters

As outlined above, we locally calibrate DDFs for snow and ice based on PDD_CONST, PDD_VAR and PDD_EFF. We find a strong sensitivity of DDFs on the individual PDD approximation (Fig. 4).

Fig. 4. Spatial distribution of DDF_,snow (upper panels) and DDF_,ice (lower panels) as calibrated locally based on PDD_VAR (left panels), PDD_CONST (center panels) and PDD_EFF (right panels). The black contour line indicates the runoff line which limits the slush zone as defined in the text.

Using PDD_VAR for calibrations yields the strongest spatial heterogeneity with DDF_{VAR, snow} and DDF_ice mostly ranging from 3 to 40 mm °C⁻¹ d⁻¹ in the slush zone. In contrast, the calibration based on PDD_EFF and PDD_CONST yields more homogenous distribution in the slush zone. DDF_{EFF, snow} and DDF_{EFF, ice} mostly remain in the range between 2.5–15 mm °C⁻¹ d⁻¹ and DDF_{CONST, snow} and DDF_{CONST, ice} in the range from 2 to 10 mm °C⁻¹ d⁻¹.

DDFs calibrated for the realization with constant standard deviation (DDF_CONST) exhibit primarily a latitudinal dependence (Figs 4, 5). If the realization incorporates the local standard deviations (DDF_VAR), calibrations primarily appear to be strongly temperature-dependent with higher DDFs at lower temperatures, while some latitudinal dependence persists. Parameters in the realization with lower bounded daily temperature variability (DDF_EFF) also exhibit a combination of temperature dependence, which is not as pronounced as in the DDF_VAR, and latitudinal dependence.

Fig. 5. DDF_,snow (upper panels) and DDF_,ice (lower panels) as calibrated locally based on PDD_VAR (left panels), PDD_CONST (center panels) and PDD_EFF (right panels) as functions of the local climatologic July temperature, T _July, of the calibration period. Each scatter point represents one grid point in the slush zone, colors reflect latitude.

3.2. Optimal parameters for Greenland-wide applications

To address the skill of different realizations of the PDD model we conduct a Greenland-wide calibration using all grid points that are within the domain of GrIS. Assuming a constant standard deviation we obtain DDF_{CONST, snow} = 5.1 mm °C⁻¹ d⁻¹ and DDF_{CONST, ice} = 5.4 mm °C⁻¹ d⁻¹ which is in the range of observation-based calibrations of high latitude ice sheets and glaciers. DDF_EFF and DDF_VAR are calibrated to 6.4 mm °C⁻¹ d⁻¹ and 10.8 mm °C⁻¹ d⁻¹ for snow and to 6.1 mm °C⁻¹ d⁻¹ and 8.1 mm °C⁻¹ d⁻¹ for ice, respectively. DDF_snow is not substantially smaller than DDF_ice, as it is in the local calibrations, because ice melt occurs largely at low elevations with relatively warm summer temperatures, while snow melt happens to largest parts under colder temperatures as in spring or higher elevations. The Greenland-wide calibration even results in DDF_,snow ≥ DDF_,ice for formulation VAR and EFF. An overview of the Greenland-wide DDFs is presented in Table 1.

Table 1. Overview of PDD model realizations. The evaluated realizations are based on approximations of PDD which differ in the representation of sub-monthly variability σ and in greenland-wide calibrated degree-day factors DDF_,snow and DDF_,ice. RMSE_time is the root mean square error of the 1948–2016 total Greenland surface melt predictions (Fig. 7) relative to the MAR simulation. RMSE_space is the root mean square error of predicted mean surface melt rates in the slush zone during the calibration period (Fig. 9) relative to the MAR simulation

Additional to the three optimally calibrated realization we define the fourth realization which is also based on PDD_VAR but which include temperature-dependent instead of constant Greenland-wide parameters to better match the local calibrations (Fig. 6). Using an analogue approach to the parameterization proposed by Fausto (2009) we find a good fit of the locally calibrated DDF_,snow with

(13)$$\eqalign{& {\rm DD}{\rm F}_{{\rm VARF},{\rm snow}} = \cr & \left\{ {\matrix{ {\beta _{{\rm snow}}^w } \hfill & {T_{{\rm July}} \ge T_w} \cr {\beta _{{\rm snow}}^w + \displaystyle{{\beta _{{\rm snow}}^c - \beta _{{\rm snow}}^w } \over {(T_w - T_c)}}(T_{{\rm July}} - T_c)} & {T_c \le T_{{\rm July}}{\rm < }T_w} \cr {\beta _{{\rm snow}}^c } \hfill & {T_{{\rm July}}{\rm < }T_c{\rm }} \cr } } \right.} $$

and DDF_,ice with the third-order polynomial:

(14)$$\eqalign{& {\rm DD}{\rm F}_{{\rm VARF},{\rm ice}} = \cr & \left\{ {\matrix{ {\beta _{{\rm ice}}^w ,} \hfill & {T_{{\rm July}} \ge T_w} \cr {\beta _{{\rm ice}}^w + \displaystyle{{\beta _{{\rm ice}}^c - \beta _{{\rm snow}}^w } \over {{(T_w - T_{{\rm July}})}^3}} * ({(T_w - T_{{\rm July}})}^3),} & {T_c \le T_{{\rm July}}{\rm < }T_w} \cr {\beta _{{\rm ice}}^c ,} \hfill & {T_{{\rm July}}{\rm < }T_c} \cr } } \right.} $$

with T _July being the mean July temperature and $\beta ^{w}_{\rm snow} = 5\,{\rm mm}\,^{\circ}{\rm C}^{-1}\,{\rm d}^{-1}$, $\beta ^{c}_{\rm snow} = 14\,{\rm mm}\,^{\circ}{\rm C}^{-1}\,{\rm d}^{-1}$, $\beta ^{w}_{{\rm ice}} =6\,{\rm mm}\,^{\circ}{\rm C}^{-1}\,{\rm d}^{-1},\, \beta ^{c}_{{\rm ice}} = 20\,{\rm mm}\,^\circ {\rm C}^{-1}\,{\rm d}^{-1}$, T ^w = 6°C and T _c = −1°C (Fig. 6) compared with the parameters $\beta ^{w}_{\rm snow} = \beta ^{w}_{\rm snow} = 3\,{\rm mm}\,^\circ {\rm C}^{-1}\,{\rm d}^{-1}$ and $\beta ^{w}_{{\rm ice}} = 7\,{\rm mm}\,^\circ {\rm C}^{-1}\,{\rm d}^{-1},\, \beta ^{c}_{{\rm ice}} = 15\,{\rm mm}\,^\circ {\rm C}^{-1}\,{\rm d}^{-1}$, T ^w = 10°C and T ^c = −1°C in Fausto (2009).

Fig. 6. Temperature dependency of DDF_{VARF, snow} (left, red) and DDF_{VARF, ice} (right, red) in comparison with the respective local calibration DDF_VAR(i) (black) for all grid points i in the slush zone (Fig. 5, left panels). The size of scatter points reflects annual mean surface melt.

We thus consider four realizations of the PDD model,

(15)$$\eqalign{& M_{\rm VAR} = {\rm DD}{\rm F}_{\rm VAR} \times {\rm PD}{\rm D}_{\rm VAR}, \cr & M_{\rm CONST} = {\rm DD}{\rm F}_{\rm CONST} \times {\rm PD}{\rm D}_{\rm CONST}, \cr & M_{\rm EFF} = {\rm DD}{\rm F}_{\rm EFF} \times {\rm PD}{\rm D}_{\rm EFF}, \cr & M_{\rm VARF} = {\rm DD}{\rm F}_{\rm VARF} \times {\rm PD}{\rm D}_{\rm VAR},} $$

which allow to estimate monthly surface melt rates from monthly climate variables of the MAR simulation, i.e. from monthly mean temperatures T̄, monthly standard deviation of daily mean temperatures σ _month, monthly mean daily temperature amplitude, dT, snow height, SNH and snow fall, SF.

3.3. Evaluation of Greenland-wide applications

The MAR simulation covers the years 1948–2016 and, averaged over this 69-year period, hindcasts a total annual Greenland surface melt of 489 Gt. In the first 50 years of the simulation, the total surface melt never exceeds 600 Gt (Fig. 7a). Towards the end of the simulation, surface melt intensifies and exceeds 600 Gt in 7 out of the last 19 years. As this intensification is not covered by the calibration period, total annual Greenland surface melt during the first 25 years amounts to 467 Gt on average (Fig. 8). The year 2012 of the simulation exhibits an exceptional melt of 987 Gt and exceeds all other years of the simulation by more than 200 Gt. Together with snow fall and refreezing this corresponds to a net surface mass loss of 427 Gt over the summer months. By comparison, the GRACE measurements indicate a total summer mass loss (including marine discharge) of 628 ± 96 Gt (Tedesco and others, 2013). In observations, the 12 of July 2012 stands out due to an extreme melt event which affected almost the entire surface of the GrIS (Nghiem and others, 2012).

Fig. 7. (a) Total Greenland surface melt from 1948 to 2016 as simulated by MAR (M _MAR, black) and predicted from M _VAR (green), M _CONST (red), M _EFF (blue) and M _VARF (cyan). (b) relative yearly bias of M _VAR (green), M _CONST (red), M _EFF (blue) and M _VARF (cyan) 1948–2016 relative to M _MAR .

Fig. 8. Mean yearly Greenland surface melt of the calibration period (years 1948–1972), as simulated by MAR. The black contour indicates the upper boundary of the slush zone. Numbers refer to the respective spatially-integrated yearly mean surface melt in the four sectors which are separated by the dotted lines and the mean total surface melt (lower right corner) in Gt.

We find that all four realizations generally capture the yearly evolution of the total yearly Greenland surface melt. Within the 69 years of the simulation the relative bias exceeds ± 10% in 5, 6 and 4 years for M _VAR, M _CONST, M _EFF, respectively, but in 20 years for M _VARF (Fig. 7b). M _CONST particularly fails to reproduce the extreme melt in year 2012 with a bias of 180 Gt. There is no indication, however, that M _CONST becomes increasingly unreliable with the intensifying surface melt towards the end of the simulation. The root mean square error, RMSE_time of the surface melt hindcast is smallest for M _EFF and largest for M _VARF (Table 1). In general, all four realizations reproduce the melt evolution relatively well. We want to stress, however, that surface melt estimates will be heavily offset if the calibration is not adapted to the PDD formulation: DDF_VAR × PDD_CONST or DDF_CONST × PDD_VAR produce a relative bias of ≈ 100% and −50 %, respectively throughout the simulation (not shown).

Throughout the calibration period of the MAR simulation, surface melt M _MAR is characterized by a narrow slush zone at the lower part of the GrIS with melt rates exceeding 3 m a^–1 near the western and south-eastern coast (Fig. 8). The spatial distribution of surface melt is reasonably well represented in all four realizations. Figure 9 shows the respective mean local bias of surface melt rates with respect to surface melt rate from the MAR simulation (Fig. 8) averaged over the calibration period. Inherent to their design, the bias in total Greenland-wide surface melt hindcast vanishes for M _VAR, M _EFF and M _CONST over the calibration period (Fig. 9); for M _VARF the total bias does not vanish but is small (≈1 Gt). Outside of the calibration period the spatial pattern remains generally the same while the total bias depends on the selected evaluation period (not shown). All surface melt estimates, especially M _CONST, tend to overestimate surface melt in northern Greenland and underestimate surface melt in southern Greenland, which is in line with the latitude-dependence found in the locally-calibrated DDFs (Fig. 5). A negative bias is particularly prominent in the ablation zone of south-western Greenland. In agreement with the temperature dependence found for the locally-calibrated DDF_VARs (Fig. 5), M _VAR tends to overestimate surface melt at low elevations and underestimate at high elevations. This is also reflected in a relatively large spatial root mean square of the surface melt distribuion in the slush zone averaged during the calibration period, RMSE_space (Table 1). The MAR simulation exhibits some residual surface melt rates in the interior of southern Greenland which usually do not exceed the refreezing rates and are thus negligible for the SMB (Fig. 8). It should be yet noted that only the realization incorporating local variability and temperature-dependent parameters (M _VARF) is able to reproduce this melt at high altitudes due to its temperature-dependent parameterization. The extreme melt of the year 2012, which is characterized by intensified surface melt at high altitudes is thus better captured in M _VARF than in the other realizations. Overall M _VARF provides the best representation of local melt rates which is also reflected in the smallest RMSE_space as presented in Table 1.

Fig. 9. (a) Bias of M _VAR prediction of Greenland surface melt to MAR simulation during the calibration period. (b) Same as (a) but for M _CONST. (c) Same as (a) but for M _EFF. (d) Same as (a) but for M _VARF. The grey contour indicates the upper boundary of he slush zone. Numbers refer to the bias in the four sectors which are seperated by the dotted lines and the mean total surface melt (lower right corner) in Gt.

4. DISCUSSION AND CONCLUSION

The PDD model is used to estimate the surface melt of ice sheets by assuming that the relation between PDD and surface melt is linear. On a continental scale, however, this assumption appears not to be justified as model parameters were found to be spatially incoherent (Hock, 2003) and in particular temperature-dependent (Tarasov and Peltier, 2002). Furthermore, previous studies have demonstrated that surface melt estimates derived from the PDD model are highly sensitive to the representation of sub-monthly temperature variability, for a given parameterization (Charbit and others, 2013; Seguinot, 2013; Rogozhina and Rau, 2014). Consequently, the suitability of the originally simple PDD model has been contested and multiple modifications have been proposed (Tarasov and Peltier, 2002; Greve, 2005; Fausto and others, 2009).

In this study, we used the regional climate model MAR as a reference to investigate how the representation of temperature variability influences characteristics and accuracy of different realizations of the PDD model.

Local calibration of a realization that incorporates sub-monthly temperature variability at each grid point yields a wide range of parameters. These parameters are primarily temperature-dependent and cover a range, which is in agreement with Hock (2003) and Fausto and others (2009). The temperature dependence is found to be reduced in two realizations which particularly overestimate the dampened variability in the ablation zone, by either applying a lower bound on diurnal variability or by representing sub-monthly variability by one constant value. Based on this observation, we propose that it is the negative feedback between melt and temperature variability that introduces non-linearity into the melt–PDD relation.

To evaluate their skill, we used the different realizations of the PDD model to hindcast MAR surface melt from MAR temperatures. We also evaluated a realization which makes use of temperature-dependent parameters instead of constant model parameters (inspired by Fausto and others, 2009). If calibrated individually (instead of using the parameters canonically used in Earth system modelling which follow Huybrechts and others, 1991) all formulations are found to generally reproduce the melt rates from the MAR simulation. The formulation inspired by Fausto and others (2009) yields the best agreement of the spatial distribution of multi-year mean surface melt but exhibits the strongest biases in the comparison of year to year total Greenland surface melt. A realization that partly neglects the dampened variability due to surface melt by restricting the daily temperature variability by a lower limit appears to be a more convenient and yet similarly skillful alternative. Our study suggests that domain-wide parameterization is suitable for applications that neglect or suppress the dampened variability of melting surfaces; otherwise calibration needs to consider temperature-dependent parameters as in Fausto and others (2009).

A regional model such as MAR includes parameterizations of important processes which are poorly constrained (e.g. snow ageing, meltwater percolation). Also local conditions such as debris cover or shading cannot be represented. Nevertheless, we expect that the feedback between melt and temperature variability and its effect on PDD schemes is robust. The presented calibrations, however, may be model-dependent. The skill of all realizations seems to remain stable throughout the reference simulation while surface melt intensifies towards the end of the 1948–2016 simulation. Further research is needed, however, to evaluate the PDD model under entirely different climates or for other ice sheets such as the North American ice sheets during the last glacial period.

Seguinot (2013) and Rogozhina and Rau (2014) demonstrate that using a PDD model that incorporates a constant mean standard deviation instead of acounting for the locally reduced temperature variability results in strongly overestimated surface melt. This is however not inconsistent with our findings. In the above studies the different realizations are not calibrated individually, but in all cases adopt the parameters canonically used in Earth system modelling (Huybrechts and others, 1991). The individually calibrated parameters presented here differ from the canonical parameters. In particular, the calibration of the realization that uses a constant standard deviation yields a parameter for ice melt that is substantially smaller than the canonical value. As inconsistent calibration was found to result in substantial errors, the disagreement between the different realizations in the above studies may well be a result of a insufficient calibration.

Existing realizations of the PDD model may be improved in the light of this study: alternative to restricting the diurnal cycle by a lower limit, it may be possible to directly diagnose the ‘effective’ variability (compensating or neglecting the dampening effect of phase transitions in applications using climate model output) or by considering temperature variability from non-glaciated, neighbouring locations (applications based on observations and reconstructions). Indeed, Shea and others (2009) present calibrations based on temperatures extrapolated from non-glaciated stations for nine Canadian glaciers and find these to be relatively consistent across the region and from year to year.Above all it appears vital to use a PDD formulation which is consistent to the one used for calibration of the PDD model. We expect that a consistent calibration based on regional models or on surface melt and temperature records such as van den Broeke and others (2010) will strongly improve the accuracy of the PDD model and reduce the sensitivity of ice sheet models to the incorporated PDD scheme.