Comparison of surface mass balance of ice sheets simulated by positive-degree-day method and energy balance approach

Glacial cycles of the late Quaternary are controlled by the asymmetrically varying mass balance of continental ice sheets in the Northern Hemisphere. Surface mass balance is governed by processes of ablation and accumulation. Here two ablation schemes, the positive-degree-day (PDD) method and the surface energy balance (SEB) approach, are compared in transient simulations of the last glacial cycle with the Earth system model of intermediate complexity CLIMBER-2. The standard version of the CLIMBER-2 model incorporates the SEB approach and simulates ice volume variations in reasonable agreement with paleoclimate reconstructions during the entire last glacial cycle. Using results from the standard CLIMBER-2 model version, we simulated ablation with the PDD method in offline mode by applying different combinations of three empirical parameters of the PDD scheme. We found that none of the parameter combinations allow us to simulate a surface mass balance of the American and European ice sheets that is similar to that obtained with the standard SEB method. The use of constant values for the empirical PDD parameters led either to too much ablation during the first phase of the last glacial cycle or too little ablation during the final phase. We then substituted the standard SEB scheme in CLIMBER-2 with the PDD scheme and performed a suite of fully interactive (online) simulations of the last glacial cycle with different combinations of PDD parameters. The results of these simulations confirmed the results of the offline simulations: no combination of PDD parameters realistically simulates the evolution of the ice sheets during the entire glacial cycle. The use of constant parameter values in the online simulations leads either to a buildup of too much ice volume at the end of glacial cycle or too little ice volume at the beginning. Even when the model correctly simulates global ice volume at the last glacial maximum (21 ka), it is unable to simulate complete deglaciation during the Holocene. According to our simulations, the SEB approach proves superior for simulations of glacial cycles.


Introduction
Glacial-interglacial cycles of the Quaternary are characterized by large fluctuations in the continental ice mass in the Northern Hemisphere (NH).These fluctuations result from the interplay of the processes of snow accumulation and surface ice ablation and the dynamic processes of calving and basal melt.The sum of gains and losses in ice mass constitutes the net mass balance.During the last glacial cycles, ice sheets have typically built up relatively slowly over roughly four precessional periods until glacial maximum, and thereafter they retreat rapidly over about 10 millennia.
The net surface mass balance is the volumetric change across an entire ice sheet and across a full accumulation and melt season.On existing ice sheets or glaciers, the surface mass balance can be obtained from local measurements of the amounts of snow accumulated in winter and snow and ice melted in summer.On long orbital timescales, the changing surface mass balance of the NH ice sheets is considered the main factor in ice sheet evolution during glacial cycles.
The net surface mass balance of ice sheets is equal to the difference between accumulation, which is controlled by the hydrological cycle, and ablation, which is determined by the surface energy balance (SEB).The SEB primarily depends on the absorption of insolation reaching the ice sheet surface and on air temperature.Numerical modelling suggests that both the accumulation and the ablation of the ma-Published by Copernicus Publications on behalf of the European Geosciences Union.
jor ice sheets in America and Europe vary in the range of 0.05 to 0.2 Sv (1 Sv = 10 6 m 3 s −1 ) for most of the glacial time (Ganopolski et al., 2010).This means that the surface mass balance is highly sensitive to small changes in accumulation and ablation, and a successful simulation of glacial cycles crucially depends on adequate descriptions of the accumulation and ablation processes.Difficulties in describing these processes arise from the non-linear nature of the climate system and from insufficient data which are needed to constrain the model parameters.
Two methods are widely used to simulate the surface mass balance of ice sheets.One method is the so-called positivedegree-day (PDD) method.This semi-empirical parameterization calculates only ablation and requires information about surface air temperature (usually, monthly mean values are used) and annual snow accumulation.This method is computationally fast and therefore widely used to compute the surface mass balance of ice sheets both in past (Tarasov andPeltier, 1999, 2002;Zweck andHuybrechts, 2003, 2005;Charbit et al., 2007;Abe-Ouchi et al., 2007;Lunt et al., 2008;Gregoire et al., 2012;Beghin et al., 2014;Liakka et al., 2016) and in future climate simulations (van de Wal and Oerlemans, 1997; Huybrechts and de Wolde, 1999;Greve, 2000;Huybrechts et al., 2004;Ridley et al., 2005;Charbit et al., 2008;Winkelmann et al., 2015).The PDD method can be calibrated with measurements from glacier surfaces, but different glaciers give different values for the PDD scaling parameters.
The other method is the physically based SEB method, which computes the melting of snow and ice from a surplus in the surface energy balance when the ice sheet surface temperature reaches melting point.This method requires calculations of all components of the energy balance (shortwave and long-wave radiation, sensible and latent heat fluxes), which in turn requires a complete set of meteorological conditions.This method is computationally much more demanding than the PDD method and was therefore used mostly in the framework of regional climate models for short-term climate predictions (Bougamont et al., 2006;Box et al., 2006Box et al., , 2012;;Fettweis, 2007Fettweis, , 2013;;Ettema et al., 2009).However, simulations with a comprehensive Earth system model have demonstrated that feedbacks between climate and ice sheets, which are not accounted for by the PDD method, are important for simulating the ice mass balance in future climate change scenarios (Vizcaino et al., 2010).
In spite of the obvious advantages of the PDD method for modelling the long-term interaction between the climate and the ice sheet, there is also a growing body of evidence that the PDD method is inadequate for modelling Quaternary glacial cycles.One obvious problem is that the PDD method does not explicitly account for the absorption of shortwave radiation, which represents the main energy component of the SEB.This can lead to a significant underestimation of the effect of the varying insolation on orbital timescales, which is the primary driver of the glacial cycles (Robinson et al., 2010;van de Berg, 2011;Ullman et al., 2015).Numerical parameters for the PDD method can only be derived from observations over the existing ice sheets, primarily in Greenland, and it is unclear a priori how different such parameters should be when the PDD method is applied to completely different climate conditions and different geographical distributions of ice sheets during glacial times.
Another semi-empirical approach, the ITM (insolationtemperature-melt) scheme, does explicitly account for absorption of insolation and reveals reasonable agreement with the SEB method in the simulation of Greenland ice sheet surface mass balance for the Eemian interglacial (Robinson et al., 2011;Robinson and Goelzer, 2014).However, ITM requires the prescription of "atmospheric transmissivity", which varies strongly in space and time and is not known for climates different from the present one.In addition, ITM does not account for the effect of dust deposition on surface albedo.This could be a serious disadvantage since paleoclimate data indicate significant increases in aeolian dust deposition during glacial times, especially along the southern margins of the NH ice sheets (Kohfeld and Harrison, 2001;Mahowald et al., 2006).Both theoretical analysis (Warren and Wiscombe, 1980;Aoki et al., 2011) and direct measurements (Painter et al., 2010(Painter et al., , 2012;;Skiles et al., 2012;Bryant et al., 2013;Doherty et al., 2013Doherty et al., , 2014;;Gautam et al., 2013) demonstrate that even small amounts of impurities affect the surface albedo significantly.In turn, the results from the SEB simulations show that these changes in albedo might significantly affect the surface mass balance of ice sheets during glacial times (Krinner et al., 2006;Ganopolski et al., 2010) and in future climate change scenarios (Dumont et al., 2014;Goelles et al., 2015).Charbit et al. (2013) discuss the effect of different PDD parameterizations on the NH ice sheet evolution, but a direct comparison between the PDD and SEB approaches in a transient simulation over the glacial cycle with a climate-icesheet model has not been performed yet.Using results from an ensemble of transient simulations of the last glacial cycle performed with an Earth system model of intermediate complexity, we undertake a systematic comparison of ice sheet surface mass balance simulated using the SEB and PDD approaches for different ice sheets and during different periods of the last glacial cycle.

Model set-up
The set-up of the Earth system model of intermediate complexity CLIMBER-2 (Petoukhov et al., 2000;Ganopolski et al., 2001) for simulations of glacial cycles and its performance are described in Calov et al. (2005) and Ganopolski et al. (2010).This model is designed to investigate processes and their interactions in the Earth climate system over long timescales, such as Quaternary glacial cy-Clim.Past, 13, 819-832, 2017 www.clim-past.net/13/819/2017/cles, which is achieved at the expense of spatial resolution and with increased complexity.The model has already been used to study the 100 ka climatic cyclicity of the Quaternary (Ganopolski and Calov, 2011), the mineral dust cycle (Bauer and Ganopolski, 2010), the climate response to dust radiative forcing (Bauer and Ganopolski, 2014) and the impact of permafrost on the simulation of glacial cycles (Willeit and Ganopolski, 2015).CLIMBER-2 consists of interactively coupled models of the atmosphere, the ocean, the land surface, the vegetation and the ice sheets.The atmospheric fields are computed on a longitude-latitude grid containing 7 × 18 grid cells.The 3-D polythermal ice sheet model SICOPOLIS operates on the NH between 21 and 85.5 • N on a longitude-latitude grid (x s , y s ) with a resolution of (1.5 • , 0.75 • ).Thus one atmospheric grid cell can overlap with more than 450 grid cells of the ice sheet model.CLIMBER-2 computes the atmospheric fields with a daily time step, the oceanic fields every 5 days and the vegetation distribution every year.SICOPOLIS computes the ice sheet evolution from losses and gains in ice mass over a 1year period.The climate component and SICOPOLIS are coupled once every 10 years through the interface module SEMI (Surface Energy and Mass balance Interface).SEMI performs physically based 3-D downscaling of climatological fields from a coarse atmospheric grid to the ice sheet model grid and computes the surface mass balance and the surface temperature using the SEB approach with a 3-day time step.Computed annual fields of surface ice sheet mass balance and surface temperature are used in SICOPOLIS.In turn, SICOPOLIS feeds back into the climate component: the ice sheet elevation, the fraction of land area covered by ice sheets, the sea level and the freshwater flux into the ocean from the ablation of ice sheets and ice calving.

Surface energy and mass balance interface (SEMI)
The interface module SEMI computes the surface mass balance on the SICOPOLIS grid (Calov et al., 2005).The surface mass balance F SEB (x s , y s ) is defined by where P (x s , y s ) is the snow accumulation and A SEB (x s , y s ) is the surface ablation (positively defined), which is hereafter called SEB-derived ablation.In SEMI, prognostic equations for ice surface temperature and snow layer thickness are solved based on the surface energy balance.The SEB is comprised of shortwave and long-wave radiative fluxes and turbulent energy fluxes.These fluxes are calculated through the horizontal and vertical interpolation of the climatological fields computed by the coarse-resolution atmospheric component.
A SEB is computed from a surplus in SEB values, which contain an explicit dependence on snow albedo.Here, snow albedo refers to broadband albedo composed of contributions from visible and near-infrared bands.Snow albedo is a function of snow aging and dust deposition (Warren and Wiscombe, 1980).Effective snow age in CLIMBER-2 is a function of temperature and snowfall.Dust deposition is composed of aeolian dust transported from remote desert regions and glaciogenic dust from the glacial erosion processes (Ganopolski et al., 2010).The computation of P includes the elevation-desert effect, which causes decreasing P with increasing ice sheet elevation, and the elevation-slope effect, which causes increasing P with the increasing slope of the ice sheet surface.The slope effects also depend on the wind direction (Calov et al., 2005).Sublimation is neglected.

Positive-degree-day (PDD) method
The PDD method is based on the reasoning that ablation is driven by the annual sum of positive temperature values, which is seen as a proxy for melt energy (Braithwaite, 1984;Braithwaite and Olsen, 1989;Reeh, 1991).The semiempirical PDD method is represented by a linear relationship using PDD values and proportionality factors for snowmelt and ice melt.The PDD value (in • C d) is defined as an excess of daily surface air temperature above the melting point accumulated over 1 year.Most implementations of the PDD method take daily temperature values from interpolated monthly mean climatological data.To account for the missing diurnal cycle and synoptic variability, a temperature variability term is included because the short-term temperature variability may implicate melt occurrences, even if the mean temperature is negative.
The PDD value is computed as the integral over time t with where t = 1 year, T (in • C) is the climatological mean surface air temperature, erfc(x) is the complementary error function and σ is the standard deviation of daily temperature from the climatological mean value (Calov and Greve, 2005).Usually, σ is prescribed in the range of 4.5-5.5 • C (Reeh, 1991;Ritz et al., 1997;Tarasov andPeltier, 1999, 2002;Greve, 2005).Fausto et al. (2009) analysed observations and showed that σ for the Greenland ice sheet may increase from 1.6 to 5.2 • C for an altitude increase from 0 up to 3000 m.The PDD-derived ablation is defined analogous to Eq. (1): where the snow accumulation P (x s , y s ) in Eqs.
(1) and ( 3) is computed in SEMI.The surface mass balance F PDD (in mm yr −1 ) is calculated by www.clim-past.net/13/819/2017/Clim.Past, 13, 819-832, 2017 where α S and α I (in mm • C −1 d −1 ) are the melt factors of snow and ice, respectively, and r S = 0.3 is a constant refreezing factor.This factor is introduced for the nocturnal refreezing of snow and causes a slowdown of the snowmelt.The factor Q (in • C d yr −1 ) is the actual remainder of PDD per year t: where PDD S is which represents the PDD value required to melt the annual accumulated snow P .The sign of Q determines the sign of the surface mass balance F PDD .When the PDD value (Eq.2) is too small to melt the available snow, then the remaining snow at the end of the year builds ice mass and F PDD is positive.Conversely, when the PDD value is large enough to melt all snow in the grid cell, then surface ice is melted and F PDD is negative.
The values of the melt factors in the PDD scheme which are suitable for a realistic simulation of ice sheets over the entire glacial cycle are not known (Hock, 2003).In the following, we attempted to find a unique set of three empirical parameters of the PDD scheme which are optimal for this task.To this end, we used the PDD scheme to simulate ablation in the "offline" mode and then in the "online" mode.In the first case, the PDD scheme is used to calculate the annual ablation rate in parallel with the standard SEB scheme employed in SEMI.Ablation simulated with the PDD scheme does not affect ice sheet evolution and is only used for comparison with the standard SEB scheme.Note that this approach is fully equivalent to the standard "offline" technique, in which temperature and precipitation fields are stored in the process of simulations with the standard CLIMBER-2 model and then only used to simulate surface mass balance with the PDD scheme.In the online mode, the SEB scheme of the SEMI module is disabled and ablation is computed with the PDD scheme.Note that in both cases (online and offline) accumulation is computed in the same way but precipitation fields are not the same for these two methods because precipitation also depends on ice sheet distribution and elevation, which are not the same in online and offline simulations.The offline and online modes are both useful in comparing different ablation schemes because in offline mode both schemes are forced by identical climate fields; however, it does not indicate how much differences in simulated ablation would affect ice sheet evolution.In online simulation, the comparison of the two ablation schemes is complicated by strong non-linearity in the climate-cryosphere system; even small differences in the forcings can lead to dramatic differences in the system response on long timescales.

Reference simulation of the last glacial cycle
The reference simulation of the last glacial cycle is driven by the insolation calculated from the varying orbital parameters (Berger, 1978) and the time-varying concentration of greenhouse gases (Fig. 1a) expressed as equivalent CO 2 concentration (Ganopolski et al., 2010).The initial condition is the equilibrium climate state computed with the greenhouse gas concentration and orbital forcing of the pre-industrial period with the Greenland ice sheet as the only NH ice sheet.The shortwave radiative forcing by aeolian dust and the deposition of desert dust on the snow of ice sheets are computed by using time slice simulations from a general circulation model.The horizontal fields of the time slices are transformed to temporally varying fields by scaling the time slices with the simulated ice volume (Calov et al., 2005).The dust deposition on ice sheets includes further dust from internally simulated sediments produced by glacial erosion (Ganopolski et al., 2010).Note that in online simulations, both the dust radiative forcing and the snow albedo differ from those in offline experiments.
Figure 1b  5.5 • C within about 10 kyr until the early Holocene.The global precipitation is thermodynamically controlled and varies in close relationship to T (Fig. 1b). Figure 1c shows the mean sea level variations computed from the NH ice volume (assuming a constant ocean surface area and an additional 10 % contribution from the Antarctic ice sheet) in comparison to the global mean sea level reconstruction (Waelbroeck et al., 2002).
Figure 2 shows the characteristics of the NH ice sheets by comparing NH total values with values from the American and the European ice sheets, which respectively represent all ice sheets in North America and Eurasia up to 120 • E. Note that the Greenland ice sheet is not included in the selections but contributes to globally averaged values.Up to 70 % of the total ice-covered area occurs in America and less than 20 % occurs in Europe (Fig. 2a).The total ice volume, given in meter sea level equivalent (msle), varies about proportionally to the total ice sheet area (Fig. 2b).The area and the volume vary in parallel with the precession-and obliquitydriven variations in the northern summer insolation.
Figure 2c shows areal averages of the time-varying ice sheet thickness.During the interglacial periods, the relatively high average ice thickness over the NH is related to the persisting Greenland ice sheet.In the initial millennia of glacial inception, the drop in the average ice thickness results from the fast spreading of the ice sheet area (Calov et al., 2005).Thereafter the average thickness of the American ice sheet increases, stays high beyond the LGM and drops rapidly toward the beginning of the Holocene.The European ice sheet thickness starts to grow at glacial inception a few millennia before the American ice thickness.Around the LGM, the European ice thickness increases by about 30 %, which is ac-companied by an extra cooling over the northern Atlantic.The lead of the thinning of the European ice sheet compared to the American ice sheet at glacial termination is attributed to the lower elevation of the European ice sheet, which facilitates the ice melt.Yet, during glacial termination the thinning of the American ice sheet occurs more rapidly than the European ice sheet.
The time series of snow accumulation (Fig. 2d) and surface ablation (Fig. 2e) vary in comparable ranges.P is well correlated with the ice sheet area and varies with the precessional period in a rather linear manner.A SEB varies in response to different driving factors, such as insolation and surface ice area exposed to temperature above the melting point and albedo of the snow surface.The maximum ablation after the LGM occurs in America some millennia earlier than in Europe.The lead of the maximum ablation in America is related to the larger perimeter exposed to melt conditions and the more southerly extent of the American ice sheet.The resulting surface mass balance (Fig. 2f) is positive and exceeds the calving rate (not shown) during most of the last glacial cycle, leading to the buildup of large NH ice sheets at the LGM.

Mass balance computed by PDD method in offline simulation
Over the last glacial cycle, empirical data needed to calibrate the PDD scheme are absent.We therefore considered the results of the standard, SEB-based mass balance simulations as the target for the PDD scheme.We compare ablation simulated by the PDD method with different empirical parameters with the results of the standard model version.We performed  a large set of offline simulations with the PDD scheme for which the standard deviation for temperature σ (Eq.2) and melt factors α S and α I (Eq.4) are considered tunable parameters.Each simulation was run with constant parameter values over the entire glacial cycle.

Selection of PDD parameter values
The PDD value computed with Eq. ( 2) depends on a prescribed standard deviation for temperature.In this study we used two different values, σ = 3 • C and σ = 5 • C. Figure 3 shows time series of T and the corresponding PDD values as areal averages over the ice sheets.After the Eemian at about 120 ka, the temperature averaged over the NH ice sheet area decreases by 13 • C (from −16 to −29 • C) in a time interval of nearly 100 kyr and then T recovers rapidly within about 10 kyr (Fig. 3a).The PDD values are closely correlated with T , showing a progressive decrease after glacial inception and a rapid increase during glacial termination.The averages of the PDD value for the total ice sheet lie in the ranges of 10-70 and 20-120 • C d with σ = 3 and 5 • C, respectively (Fig. 3a).The glacial cycle asymmetry is substantiated by the massive and widespread ice sheet in America, which shows a temperature evolution from −16 to −27 • C (Fig. 3b).The temperature of the smaller European ice sheet fluctuates more strongly, from −10 to −29 • C.These fluctuations are connected with changes in the sea-ice albedo effect in the northern Atlantic and changes in the Atlantic meridional overturning circulation.The PDD values for Europe range from 10 to 260 • C d with σ = 3 • C and from 30 to 370 • C d with σ = 5 • C (Fig. 3c).
Previous climate model studies have often used σ = 5 • C and the so-called standard melt factors for snow and ice, which are (α S , α I ) = (3, 8) mm • C −1 d −1 as derived from the measurements on the Greenland ice sheet (Huybrechts and de Wolde, 1999;Tarasov andPeltier, 1999, 2000).However, other observations show that melt factors may vary with latitude and the height of the glaciers.Hock (2003) summarized worldwide measurements during the melt season of glaciers and snow-covered basins and showed that (α S , α I ) may vary in the ranges of ([2.5-11.6],[5.4-20]) mm • C −1 d −1 .The ranges of the melt factors are relatively wide because they are obtained from different environments and incorporate variations in space and time, insolation, ice sheet elevation, sensible heat flux and surface albedo.Here, we consider two σ values to test the possible effect on A PDD from unresolved space-time variations in the modelled temperature.For each σ value, the values for α S and α I are varied in wide ranges.For σ = 3 • C, (α S , α I ) are varied in the ranges of ([3-10] Thus the offline PDD-derived ablation can capture the entire variability in the simulated SEBderived ablation during the glacial cycle.

Ablation time series for ice sheets over glacial cycle
In an attempt to find PDD parameter values which produce the best fit to A SEB , we calculate as measures of agreement the mean anomaly m and the RMSE r from time series of ablation averaged over ice sheets.Figure 4 shows contour plots of m and r as a function of α S and α I calculated from 130 kyr series of A PDD and A SEB for all NH ice sheets.For both σ values, no unique pair of (α S , α I ) values exists at the minimum in m (Fig. 4a, c), while the minimum in r can be associated with a specific pair of (α S , α I ) values (Fig. 4b, d).However, the minimum RMSE of about 0.025 Sv is large and amounts to more than 50 % of the peak value in A SEB simulated at 15 ka.In another attempt, we try to find optimal PDD parameter values separately for the American and European ice sheets.The contour plots of the RMSE as a function of α S and α I (Fig. 5) show that very different values of (α S , α I ) are optimal for the American and European ice sheets (Table 1).Overall, r for the American ice sheet (Fig. 5a, c) is about a factor 3 larger than for the European ice sheet (Fig. 5b, d) in both σ sets.
Figure 6 shows the PDD-derived ablation evolution for the American and European ice sheets for the entire ensemble together with the SEB-derived ablation.The agreement between the series is much lower for the American than for the European ice sheets irrespective of the σ value.Typically, if A PDD and A SEB agree better during glacial inception, then A PDD underestimates the peak in A SEB at glacial termination; conversely, if A PDD reproduces the peak in A SEB at glacial Clim.Past, 13, 819-832, 2017 www.clim-past.net/13/819/2017/termination, then A PDD overestimates A SEB at glacial inception.Hence, smaller melt factors would be needed for glacial inception than for glacial termination.So, we divide the time series into the intervals 130-30 ka and 30-0 ka and determine for each sub-interval the PDD parameter values which minimize the RMSE (Table 1).Nonetheless, ablation series fitted for the American ice sheet and for the sub-intervals diverge repeatedly from the reference series (Fig. 6a, c).In particu-   (5; 3, 16).See Table 1 for a summary of the PDD parameter values at the minima of r.
Table 1.Summary of the PDD parameters inducing minimum RM-SEs between series of offline A PDD and A SEB for the American and European ice sheets covering the entire glacial cycle (see Fig. 5), glacial phase and glacial termination (see Fig. 6).lar, the enhanced ablation from the American ice sheet during MIS 4 (ca.75-60 ka) is difficult to reproduce with the PDD method.Otherwise, the PDD-derived ablation series fitted for the European ice sheet and for the sub-intervals agree quite well with the reference, and the mismatch at 30 ka is small (Fig. 6b, d).

Geographically resolved ablation rates at 15 ka
At 15 ka, the total SEB-derived ablation reaches its maximum of 0.41 Sv (Table 2).This is why we choose this time slice to analyse the geographical distribution of ablation sim-Table 2. Ablation (in Sv) from NH total and the American and European ice sheets at glacial termination (16-14 ka) where maximum A SEB at 15 ka for NH is closely reproduced with the offline PDD method using σ = 3 • C and (α S , α I ) = (6, 19) in mm • C −1 d −1 (see Fig. 7).Maxima in ablation (bold) occur 1 millennium earlier in A SEB than in A PDD for NH total and the American ice sheets.Note that while the total ablation rates at 15 ka from both methods are close, A SEB in America is underestimated and A SEB in Europe is overestimated by the PDD method.ulated with the PDD scheme versus the standard SEB approach.The ensemble member which produces similar total ablation to the reference simulation at 15 ka is obtained with (α S , α I ) = (9, 16) mm • C −1 d −1 and σ = 3 • C (Fig. 7).
Figure 8 compares the spatial patterns of the ablation rates simulated with both methods using the ensemble which produces the same NH total ablation as the reference simulation at 15 ka (Fig. 7).The scatter diagram shows that the PDD method tends to overestimate large ablation rates and to underestimate low ablation rates.The PDD-derived American melt rates overestimate the reference melt rates larger than ∼ 10 mm d −1 but underestimate the American ice melt rates less than ∼ 8 mm d −1 (Fig. 8a).The PDD-derived European melt rates are overestimated mainly for ablation rates larger than ∼ 6 mm d −1 (Fig. 8b).The largest ablation rates occur naturally at the ice sheet margins, and the largest differences between the two methods used here also occur at the ice sheet margins (Fig. 9).

Glacial cycle simulations with online PDD method
Above, we evaluated the PDD-derived ablation from offline simulations against the SEB-derived ablation.In doing so we explicitly assumed that the latter provides a realistic spatial and temporal distribution of ablation because in the reference simulation ice sheet evolution during the last glacial cycle is in reasonably good agreement with paleoclimate reconstructions.In offline simulations we found that ablation simulated with the PDD scheme generally deviates from that simulated with the standard SEB approach.To assess how these differences will influence ice sheet evolution during the last glacial cycle, we performed a set of PDD online simulations, in which the PDD scheme for ablation replaces the standard SEB scheme.Note that the accumulation scheme remains the same in these simulations.We evaluated the PDD online simulations by comparing their results with the reconstructed global sea level and climate characteristics from the reference simulation.

Selection of PDD parameter values
A few dozen glacial cycle simulations were performed with the online application of the PDD method.In these experiments, we tested how well the evolution of ice sheets and climate can be simulated with constant PDD parameter values.It appears that the values of three PDD parameters can be tuned adequately for certain time periods but not for the www.clim-past.net/13/819/2017/Clim.Past, 13, 819-832, 2017 entire glacial cycle.The PDD online simulations can be split into two clusters (Fig. 10).In the first, the sea level is reasonably simulated during glacial inception (from 120 ka until about 110 ka), but it diverges dramatically from the paleoclimate reconstruction for the rest of glacial cycle.In particular, all of these simulations fail to produce deglaciation toward the end of the Holocene.Simulations from the second cluster are able to reproduce complete deglaciation before the end of the experiment but significantly underestimate ice sheet volume during most of the glacial cycle.In the following, we show representative simulations from the two clusters with the parameter values given in Table 3.The target pe-riods, inception and termination, are seen to impose a rather weak constraint for selecting the PDD parameter values.In contrast, the target period LGM (21 ka) emerged as a rather strong empirical constraint.Only one specific pair of melt factor values for each σ value (Table 3) is found suitable to simulate the LGM climate with the online PDD method.

Target periods: glacial inception and termination
During glacial inception (from about 120 until 110 ka) PDD online simulations I3 and I5 (Table 3; blue line in Fig. 10) closely reproduce the global temperature (Fig. 10a, c) and the sea level (Fig. 10b, d).In this time interval, the ice sheet area Clim.Past, 13, 819-832, 2017 www.clim-past.net/13/819/2017/grows sufficiently fast together with accumulation.Thereafter the ice volume grows too fast in combination with amplified snow accumulation, and the simulations drift into an excessively cold climate state.At 21 ka in these experiments, the ice volume is about twice as large as reconstructed (Table 3), and simulations I3 and I5 fail to terminate the glacial climate state.Contrary to the experiments described above, the temperature and the sea level in simulations T3 and T5 (red line in Fig. 10) successfully simulate complete deglaciation after a weak glacial phase (Table 3).The global cooling after inception is about in phase with the reference temperature, though the cooling in the PDD online simulations is substantially underestimated (Fig. 10a, c).The sea level drop in simulation T3 is about half as large as reconstructed over the glacial phase (Fig. 10b), and in simulation T5, the maximum sea level drop of 40 m occurs after the LGM (Fig. 10d).From 38 to 20 ka, the cooling rate in both simulations T3 and T5 intensifies and thereby the ice volume grows continuously beyond 21 ka until around 18 ka.Therefore, both simulations T3 and T5 substantially undershoot the buildup of the ice volume.

Target period:
LGM PDD online simulations L3 and L5 (green line in Fig. 10) reproduce the reconstructed sea level at 21 ka (Table 3) reasonably well.In the initial phase of the glacial cycle, simulation L3 produces a weaker cooling and less ice volume than the reference experiment, but in the time interval 40-21 ka the agreement is close (Fig. 10a, b).Simulation L5 generates growing ice sheets over the entire glacial phase, which agrees well within uncertainties inferred from the reference and the reconstructed sea level (Fig. 10c, d).However, in simulations L3 and L5, the ice volume continues to grow beyond 21 ka.
Consequently, glacial termination is delayed and complete deglaciation is not achieved.
The geographic distribution of the ice sheet thickness at 21 ka from PDD online simulation L3 agrees closely with the reference simulation (Fig. 11).Simulation L3 reproduces the maximum thickness of 3500 m in America as simulated by the reference, but in simulation L3 the ice sheet spreads slightly more southward beyond the American Great Lakes, and the ice sheet in the European Arctic and in northeastern Asia is slightly thinner.Also, simulation L5 produces an ice sheet distribution similar to the reference, although the maximum thickness in America is only 3300 m at the LGM.Both PDD online simulations L3 and L5 simulate a sea level of −120 m at the LGM, but thereafter their mass balances remain more positive than in the reference experiment.

Conclusions
In this study, we compared the simple and computationally efficient PDD scheme with the much more complex and computationally demanding SEB-based scheme implemented in the Earth system model of intermediate complexity CLIMBER-2.To this end, we performed a large set of experiments, first in offline and then in online mode.In the first case, the PDD index was computed using surface air temperature simulated with the standard CLIMBER-2 model version.The ablation computed with the PDD method was not used to simulate ice sheet evolution, and therefore the differences between the ablation simulated with the SEB and the PDD methods in the offline mode are only due to the differences between these two schemes.By comparing the ablation simulated by the different methods, we found that there is no single set of melt factors which allows the PDD scheme to simulate ablation similar to that simulated with the SEB scheme for the entire glacial cycle.Our analysis showed that for a realistic simulation of glacial termination, significantly larger PDD melt factors are required than for simulating glacial inception.Additionally, we found that larger melt factors are required for the American ice sheet compared to the European one.In general, it appears that the mass balance of the European ice sheet correlates better with the PDD than the mass balance of the American ice sheet.This suggests that the evolution of the American ice sheet is more strongly influenced by changes in absorbed shortwave radiation and surface albedo.
Similar to the offline simulations, the results of the online simulations show that no universal PDD parameter values exist for which ice volume evolution during the entire glacial cycle is satisfactorily simulated.We found that different PDD melt parameter values are required for reproducing the reconstructed ice volume during glacial inception and glacial termination.The sets of melt factors in the PDD scheme which generated a realistic ice volume rise during glacial inception led to a strong overestimation of ice volume at the LGM www.clim-past.net/13/819/2017/Clim.Past, 13, 819-832, 2017 and failure to simulate complete deglaciation.At the same time, the model versions with the PDD melt parameters that simulated correct timing of deglaciation during the Holocene strongly underestimated ice volume during the entire glacial cycle.Lastly, the PDD melt factors which allowed us to simulate correct LGM ice volume led to an underestimation of ice volume during glacial inception and a too-late onset of deglaciation.
In summary, the results of our offline and online simulations demonstrate that the PDD scheme cannot reproduce the results of the physically based SEB scheme with a constant set of model parameters.Hence, the use of the PDD method for large climate changes and geographically varying continental ice sheets, as reconstructed during glacial cycles, is found to be problematic.At the same time, the climate component of the model used in this study has a very coarse spatial resolution and the SEB-based scheme includes a number of tunable parameters; not all of them are well constrained by empirical data.Therefore, our comparison between the SEB and PDD approaches should be considered as tentative, and using higher-resolution climate models would be desirable before drawing a final conclusion.

Figure 1 .
Figure 1.Reference simulation of last glacial cycle with the CLIMBER-2 model coupled with the SICOPOLIS model via the SEB approach.(a) Driving equivalent CO 2 concentration, (b, red) global mean surface air temperature, (b, blue) global mean precipitation and (c) sea level are shown by the green line from simulated ice volume variation and by the black dashed line from reconstructions by Waelbroeck et al. (2002).

Figure 2 .
Figure 2. Glacial cycle series from the reference simulation for NH total (green lines) and the American (red lines) and European (blue lines) ice sheets showing (a) ice-covered area, (b) ice sheet volume, (c) average ice sheet thickness, (d) accumulation, (e) SEB-derived ablation and (f) surface ice mass balance.Note that the European ice sheets represent all ice sheets in northern Eurasia up to 120 • E.

Figure 3 .
Figure 3. Glacial cycle series averaged over (a) NH total and the (b) American and (c) European ice sheets showing on the left axes (red) surface air temperature and on the right axes (black) PDD values (Eq.2) computed with σ = 3 • C (dashed lines) and with σ = 5 • C (continuous lines).

Figure 4 .
Figure 4. Bivariate distributions (a, c) of mean anomaly m and (b, d) of RMSE r (in Sv) from the 130 kyr NH ablation series as a function of α S and α I using ensemble simulations of A PDD (offline) relative to A SEB .A PDD simulations (a, b) use σ = 3 • C and (c, d) σ = 5 • C; the former involves larger values for (α S , α I ).See Table 1 for PDD parameter values at a minimum of RMSE (b, d).

Figure 5 .
Figure 5. Bivariate distributions of RMSE r (in Sv) the from 130 kyr ablation series as in Fig. 4 but separately for (a, c) the American ice sheet (a, c) and for the European ice sheet (b, d).A PDD simulation (a, b) with σ = 3 • C and (c, d) with σ = 5 • C. See Table1for PDD parameter values at a minimum of RMSE in each panel.

Figure 7 .Figure 8 .
Figure 7. Ablation series from interval 30-0 ka for NH total (green lines) and for the American (red lines) and European (blue lines) ice sheets showing A SEB of the reference simulation by thick lines and offline A PDD by thin lines.A PDD with σ = 3 • C and (α S α I ) = (9, 16) mm • C −1 d −1 is compatible with A SEB at 15 ka.See Table2for peak ablation values.

Figure 9 .
Figure 9. Geographic distribution of ablation differences (in mm d −1 ) obtained from the PDD offline simulation using σ = 3 • C and (α S α I ) = (9, 16) mm • C −1 d −1 relative to the reference simulation at 15 ka; NH total ablation from the PDD and SEB methods agree closely (see Fig. 7).The thin black lines are present-day coastlines.
Table 1 for PDD parameter values at a minimum of RMSE in each panel.

Table 3 .
Global surface air temperature (T ) and sea level (sl) at 21 ka (LGM) and 0 ka (MOD) from the reference simulation (RS) compared with the PDD online simulations using σ in • C and (α S , α I ) in mm • C −1 d −1 .PDD online simulations are selected to fulfil the target windows glacial inception (I3, I5), glacial termination (T3, T5) and LGM (L3, L5) as shown in Fig.10.Note that simulation I5 uses standard PDD parameter values (bold).