The importance of insolation changes for paleo ice sheet modeling

Introduction Conclusions References


Introduction
Paleo ice sheet studies can help to constrain an ice sheet's evolution for a wide range of climatic conditions, leading to both a better understanding of past climatic changes and reduced uncertainty concerning potential future changes.The Eemian interglacial, for example, has been shown to have exhibited a considerably warmer Arctic than today (e.g.Bakker et al., 2013;Lunt et al., 2013), so while it is not a perfect ana-Correspondence to: A. Robinson (robinson@fis.ucm.es)logue for the future due to its different orbital configuration, it can still be used to learn about ice sheet sensitivity during warmer periods (Ganopolski and Robinson, 2011).Aside 35 from records obtained from several ice cores (e.g.Vinther et al., 2009;NEEM Community Members, 2013) and marine sediment cores (e.g.Colville et al., 2011), dynamic ice sheet modeling is the major source of information used to constrain past ice sheet conditions and evolution.

40
Transient coupled climate-ice sheet simulations at high spatial resolution -needed to explicitly model ice-climate interactions -are largely not feasible on glacial-interglacial time scales.Therefore, a range of climate modeling strategies is currently applied to improve our understanding of ice 45 sheet evolution during past climates.At one end of the spectrum are Regional Climate Models (RCMs) of high spatial resolution that include a full representation of the ice/snowatmosphere interface (e.g.Helsen et al., 2013).However, RCMs must be forced at the lateral boundaries by output 50 from past climate simulations with General Circulation Models (GCMs), which usually implies a fixed ice sheet surface elevation or strategies to circumvent this problem (e.g.Helsen et al., 2012;Edwards et al., 2013).Attempts have also been made to directly couple ice sheet models and GCMs 55 (e.g.Fyke et al., 2011;Lipscomb et al., 2013), but the high computational demands of both RCMs and GCMs currently limit long-term transient simulations to those using asynchronous coupling (Helsen et al., 2013) or other interpolation strategies (e.g.Stone et al., 2013).

60
At the other end of the spectrum are modeling strategies that allow for full glacial-interglacial simulations.Currently, this requires intermediate complexity and/or highly parameterized representations of some physical processes, in particular for the climate and the interface between the atmo-65 sphere/ocean and the ice sheet (e.g.Letreguilly et al., 1991;Ritz et al., 1997;Ganopolski and Calov, 2011;Robinson et al., 2011).Many paleo modeling studies parameterize icesheet surface melting as a function of temperature (e.g.Le-A. Robinson and H. Goelzer: Insolation in paleo ice sheet modeling treguilly et al., 1991;Ritz et al., 1997), by relying on the widely used Positive-Degree-Day (PDD) approach (Reeh, 1991).
In the PDD approach, insolation is not explicitly considered, which could lead to unwanted biases mainly during periods with higher than present-day insolation.This is due to the fact that PDD implicitly accounts for a certain contribution of insolation to the melt rate based on present-day insolation.In reality, when insolation increases, the proportion of melt due to insolation increases (and vice-versa), which implies that the PDD coefficients should change too.The relative importance of insolation changes during the Eemian has been estimated to be 45-50% of the total melt anomaly by applying either a simplified energy-balance equation that includes the effect of insolation (Robinson et al., 2010) or a full RCM (van de Berg et al., 2011) for peak Eemian conditions.
While these studies only consider peak insolation anomalies relative to today, it is clear that such a large contribution to melt should not be ignored.
For the present study, we extend the analysis of insolationrelated melt to the entire Eemian and Holocene (Sec. 2 and 3) to better understand how the contribution of insolation to surface melt changes for different orbital configurations.In doing so, we also derive a general parameterization that can be used to correct for the melt contribution due to insolation changes in melt models that lack that component (Sec.4 and 5).

Melt equations
In order to be able to derive analytical estimates of the contribution of insolation anomalies to total melt, we use a simplified energy-balance equation for melt that includes contributions from insolation and temperature (van den Berg et al., 2008;Robinson et al., 2010).This model has been validated with simulations of the Eurasian ice sheet (van den Berg et al., 2008) and the Greenland ice sheet (Robinson et al., 2010(Robinson et al., , 2011;;Fitzgerald et al., 2012).The insolationtemperature melt (ITM) equation approximates the melt rate, M (m • s −1 ), as In the first term, S is the incoming solar radiation at the top of the atmosphere, τ a is the atmospheric transmissivity and α s is the surface albedo.The second term c + λT is a linear parameterization of the long-wave radiation and turbulent heat fluxes , where T is the daily mean temperature and λ and c are constants.The parameter λ = 10 is chosen such that it corresponds to a melt rate of 3 mm w.e. per degree, analogous to the typical choice of degree day factor for snow in the PDD approach (Reeh, 1991).c is used as a free parameter that depends on the region being modeled (van den Berg et al., 2008) and the representation of albedo (Robinson et al., 2010).ρ w is the density of water and L m is the latent heat of melting.We calculate the daily melt rate by multiply-110 ing M by the seconds in a day.This melt equation applies locally at each grid point, thus both the insolation S and the temperature T are spatially and temporally varying terms.
The only two forcing variables in Eq. 1 are insolation and temperature, which is a great simplification of a full energy 115 balance at the surface.Nonetheless, in most climatic scenarios at the continental scale, it can be assumed that this forcing carries the information necessary to represent the energybalance well to the first order.In this context, such an equation is convenient because we can easily assess how a change 120 in either insolation or temperature would affect the calculated melt not accounting for second-order effects on atmospheric transmissivity.For example, by examination, it is clear that if insolation is held constant (e.g. at present-day levels), then its relative contribution to total melt will decrease asymp-125 totically as temperatures become higher.Conversely, at low temperatures melt is determined more strongly by insolation.
Given equation Eq. 1 for the melt rate, which is a function of insolation, temperature, transmissivity and surface albedo, it is possible to isolate the influence of changing insolation on 130 the estimated surface melt for any local conditions: where S 0 is the spatially and seasonally explicit insolation 135 fixed to present-day levels.If the surface albedo and atmospheric transmissivity are known, the value of a ( m 3 W s ) can be determined analytically: a = 1 ρwLm τ a (1 − α s ).However, all of these variables are typically not available in simple models.Assuming that over the ice sheet, τ a is in the range of 140 0.5-0.7 (Robinson et al., 2010) and α s is in the range of 0.4 (ice) to 0.8 (dry snow), a reasonable range for the coefficient a is 3-13•10 −10 m 3 W s .For a 1 W/m 2 insolation anomaly over one month, the additional melt would be between 0.8-3.3mm w.e, not counting any feedbacks.For comparison, 145 summer insolation anomalies over Greenland during glacialinterglacial cycles can range from -50 to +80 W m −2 .

Assessing the insolation-anomaly melt contribution
We performed two sets of experiments to assess the importance of insolation anomalies on melt -one set using the 150 standard ITM model (referred to as the "full model") and one set using the ITM model with insolation fixed to its present-day values (referred to as the "reduced model").For these experiments, we held the Greenland ice sheet topography constant to its present-day distribution (no ice sheet 155 model).Each transient simulation consisted of calculation of the surface melt over the Eemian ( 135-115 ka BP) and the Holocene ( 15-0 ka BP) interglacial periods.We also applied constant temperature anomalies of 1, 3 and 5 • C throughout the simulations.These idealized experiments help us to quantify each component of melt, whereas later we will show results from a realistic experimental setup with changing ice sheet geometry and transient temperature anomalies.
Figure 1 shows the time series of annual simulated Greenland ice sheet melt for several experiments.During peak Eemian conditions (in terms of summer insolation anomaly) and for a 3 • C temperature anomaly, the melt predicted by the full model (M (S)) is ca.30% higher than the melt predicted when insolation is fixed to present day values (M (S 0 )).However, this percentage depends strongly on the temperature anomaly, which is not well constrained for this time period.For 1 and 5 • C temperature anomalies, the percentage of melt due to insolation changes is 50% and 20%, respectively.
The insolation-anomaly melt contribution has a non-linear effect on the total melt via the albedo-melt feedback.To quantify the pure contribution of insolation changes to melt without this feedback, in each simulation using the full ITM model, we also simultaneously calculated the melt that would be predicted given present-day insolation.In this way, the melt in both cases is calculated given the same surface conditions (albedo, snowpack thickness) that evolve as determined by the full model.The difference between M (S) and this synthetic melt (M (S 0 ) full ) is equal to the analytical contribution of insolation to melt (M insol ) from Eq. 2. Through this decomposition, we find that the albedo-melt feedback accounts for up to 30% of the additional melt due to insolation, or 10% of the total melt anomaly, in case of a 3 • C temperature anomaly.
Figure 2 shows the spatial distribution of melt anomalies at 126 ka BP given a 3 • C temperature anomaly.As expected, higher melt anomalies are found at lower surface elevations near the margins.However, the percentage of melt due to the insolation anomaly increases towards higher elevations, where temperatures are lower.In the interior of the ice sheet where temperatures are well below the freezing point, any small amount of additional melt that occurs is due to the insolation anomaly.In intermediate regions, particularly near the equilibrium line, the insolation anomaly plays an important role in increasing melt and thus, in further decreasing the albedo.It is there that the insolation-related melt anomaly has the largest absolute value.

Applying and testing a correction term
Having estimated the melt contribution due to insolation changes, we now derive a parameterization that may be used to correct melt models that do not already account for it.The additional contribution is formulated to vanish as the insolation anomaly approaches zero, which enables the use of a classical melt model without modification for the present day and the next few millennia, when insolation changes are 210 small.This is useful because historically many paleo icesheet simulations rely on a PDD approach and could be corrected in a way consistent with previous simulations, rather than employing a different melt model (such as ITM) all together.

215
Our parameterization follows from Eq. 2, in that the insolation related melt M insol is a function of the insolation anomaly multiplied by an empirical coefficient a.The coefficient a to the first order could be considered as a constant, representing average local conditions of surface albedo 220 and atmospheric transmissivity for the ice sheet.In practice, however, such an approach leads to significant errors.Instead, we find that the surface conditions must be taken into account in the formulation of a.
By solving for a in Eq. 2, we were able to diagnose the emprical model value of a

225
for each point on the ice sheet and each month.We obtained the melt terms M (S) and M (S 0 ) full above by performing melt calculations with the full ITM model and, with the same surface conditions, the reduced ITM model (with present-day insolation).Figure 3 shows the local values of a 230 calculated from this approach for each month as a function of temperature.These results indicate that a cannot be taken as a constant, but must contain some seasonal variation.Furthermore, a high correlation between this variation and local temperature exists because the near-surface temperature acts 235 as a first-order proxy of local surface conditions (albedo), and thus of how much insolation can be absorbed at different times in the year.
To better capture this local variability without introducing additional variables, we parameterize the now spatially explicit coefficient a as a function of the local mean daily temperature (T ) and the day of the year (d): where a(T,d) is a piecewise linear function of T: ( The values of a max and T max were chosen to be constant for simplicity.T min was parameterized as a non-linear sinusoidal function of the time of year between T max in winter and T min,sum in summer (Figure 4): T min,sum and the exponent p were tuned to best match the diagnosed values of a shown in Fig. 3. See Table 1 for the parameter values used here.Defining a as a function of temperature and season allows us to capture a large part of its variability (Fig. 3).The relationship is most robust at the start of the melt season, which is more important for the initiation of melt and thus for propagating albedo-melt feedbacks throughout the year.It can be seen that a reaches values as high as 13•10 −10 as expected from the analytical calculations.For higher temperatures, it levels out to a more intermediate value of 8.3•10 −10 , which corresponds to an average surface state with the albedo of melted snow in the model (α s = 0.6).The parameterization shows deficiencies in the later part of the melt season (e.g., July and August) when the clear relationship between a and temperature breaks down somewhat.During this time, the surface in different locations has had time to evolve in a more complex response to the forcing than can be captured by the parameterization.In effect, the coefficient a is acting as a substitute for a full surface albedo model, based on the idea that a simple melt model would not have albedo available.A more complex approach that better captures the evolution of the surface conditions would be outside the scope of this first order correction term.
We would like to know if such a correction term can be applied in a melt model that does not account for changes in insolation, such as the PDD approach.To be able to compare our results consistently, we assume that the reduced ITM model, in which insolation is fixed to present day, can be used to represent melt from such an approach.This is justified because this model has previously been shown to have a similar sensitivity to climatic changes as the PDD approach given present-day insolation (Robinson et al., 2010).Based on the above derivation, we use M insol from Eq. 4 to provide a first-order correction to the melt estimated by the reduced-ITM model (M (S 0 )) for different orbital conditions, such that the corrected melt is estimated as Figure 5 shows the transient melt anomaly analogous to the middle panels of Fig. 1, now including results using the corrected melt estimate from Eq. 7. Overall, the melt correction functions as expected.Both entering and leaving the Eemian interglacial, the corrected melt is slightly overestimated compared with the full model.This implies that our formulation of a is too simplistic to fully capture the surface conditions during intermediate transition periods.Nonetheless, the corrected melt is able to reproduce the full model's melt very well.Particularly, the rate of increase in melt as the interglacial intensifies, as well as the timing and magnitude of peak melt, are consistent with the full model.There 290 is more uncertainty for peak interglacial conditions, as this is when the model has the most freedom to evolve.
We performed several simulations with different values of T min,sum as a sensitivity study of the correction term.We only modified T min,sum as it was found to have the most 295 influence on the results, while the remaining tunable parameters (a max , T max , p) have only a secondary effect (not shown).The main impact of T min,sum is to modulate the magnitude of the correction.Because the insolation correction term is only nonzero above the value of T min , as 300 T min,sum increases, the term's contribution to total melt decreases.We find that a value of T min,sum in the range of -15 • C to -12 • C approximates the contribution of insolation well, with an optimal value of T min,sum =-14 • C.
Figure 6  error in melt appears in the intermediate elevations, as insolation becomes less important at lower elevations with higher temperatures.Given a total melt anomaly of between 500-2000 mm (see Fig. 2) in these regions, it is clear that a model that does not account for insolation changes would signifi-315 cantly underestimate melt during this time period.Furthermore, it can be seen that the reduced model shows the equilibrium line elevation to be much lower than that of the full model in several regions.This means that the full model has a significantly wider ablation zone.

320
Conversely, the corrected melt model reduces this error to well below 50 mm in most regions.Only a few spurious points are not well corrected, which is a result of a small equilibrium line shift compared to the full model.In regions that do experience a significant amount of melt, at interme-325 diate elevations and towards the margins, the corrected melt model compares very well with the full model.While this may be expected from the analytical derivation, the result is encouraging since it shows that the proposed formulation of a captures the effect of changing surface conditions suffi-330 ciently well.

Coupled transient simulations
To understand the impact of insolation changes on the evolution of the ice sheet geometry, as well as to further test the ability of the new parameterization, we performed tran-335 sient simulations using the fully coupled model REMBO-SICOPOLIS (Robinson et al., 2010(Robinson et al., , 2011)).In this setup, the climate model REMBO is driven in the same way as above, but now by a transient time series of temperature anomalies as well (Fig. 7).The temperature anomalies were obtained from a transient global climate simulation of the last several glacial cycles (Ganopolski and Calov, 2011).The SMB calculated in REMBO is used to drive the ice sheet model SICOPOLIS (Greve, 1997).As before, the full version of the melt model is used (M (S)), along with the reduced version (M (S 0 )) with insolation fixed to present day and the corrected version (M corr ).
In all simulations, the ice sheet steadily grows entering the Eemian (Fig. 7) due to increased precipitation and low temperature anomalies.Its growth slows as the temperature anomaly becomes positive, and it begins to lose mass after peak climatic conditions are reached.It slowly regrows several kiloyears later under lower temperatures, with the timing of recovery depending on the amount of volume lost.Using the full melt model, the ice-sheet volume decreases about 25% relative to present day, with the minimum volume at 122 ka BP.In contrast, the reduced model gives a volume reduction of just 8% and faster regrowth, reaching its minimum volume earlier between 123-124 ka BP.There is less variation during the Holocene, largely due to the near-zero temperature anomalies applied in these simulations.Still the full model produces somewhat more melt as peak Holocene conditions are reached, leading to a smaller volume of the ice sheet throughout the Holocene as well.
In these fully coupled simulations, the corrected melt model produces ice loss almost indistinguishable from that of the full model.The additional melt provided from the insolation correction term results in significantly more ice loss throughout both the Eemian and the Holocene compared to the reduced model.Particularly, the magnitude of ice loss and the timing of the minimum and regrowth are well reproduced by the corrected model.Different values of T min,sum perform as expected, with higher parameter values corresponding to less volume reduction throughout the Eemian.The uncertainty is reduced significantly during the Holocene, given the lack of such extreme climatic forcing.The corrected model also simulates a spatial distribution of ice loss on Greenland that is very similar to that of the full model (Fig. 8).When the minimum volume is reached during the Eemian, the corrected model shows the same pattern of ice loss everywhere, and particularly in the more sensitive Northwest and Southwest regions.The consistency of these results further strengthens our confidence in the ability of our formulation of a to properly account for changing surface conditions, since in a coupled setup, small differences can compound into larger errors.

Discussion
Using the ITM equation, we were able to derive an analytical expression for the insolation-anomaly contribution given local albedo, transmissivity, temperature and insolation.As these variables are not always available in reduced complex-ity melt models, we further parameterized the coefficient a (that encapsulates the contribution of local conditions) as a function of temperature.This is, of course, a great simplification compared to a full surface albedo model.Because 395 the local temperature provides a relatively consistent (albeit crude) picture of the surface conditions, it works well to represent the albedo-melt feedback in the correction term.
Our approach should be valid over all paleoclimatic conditions, since it accounts for the local insolation anomaly and 400 the local surface conditions (via a).It has only been tested during periods warmer than today, however during glacial periods, melt on the ice sheet reduces to zero, making such a test inconsequential.Still, our formulation allows for a negative contribution from insolation when it is lower than 405 today.As seen for the latter part of the Eemian in Fig. 5, when the local insolation anomaly decreases faster than the temperature anomaly, including the negative contribution to melt from insolation makes the overall melt anomaly decline more rapidly.Such transient behavior could have important 410 implications concerning the regrowth of ice sheets after an interglacial period.
As the ITM model cannot capture the full complexity of the surface energy balance of the ice sheet, our analysis can only provide a first order estimate of insolation's role in paleo 415 melting.However, given the inherent simplicity of reduced melt models such as PDD, the estimates here should be good enough to allow improvement to be made.The application of full energy-balance models in transient experiments would allow us to refine our understanding further of how ice-sheet 420 melting may have differed during paleoclimatic conditions.
Nonetheless, it is clear that the contribution of insolation changes to melt is intimately linked to the temperature anomaly experienced at the same time.This link is tightly coupled, because an increase in either term causes a reduc-425 tion in albedo and therefore leads to a increase in melt.As found by previous studies (e.g., Fitzgerald et al., 2012), the representation of albedo is a critical component of melt and surface mass balance modeling.Here we also find that it is necessary to incorporate at least a simple treatment of albedo, 430 in order to properly simulate the changes in the snowpack.
Ignoring the insolation effect is a systematic error that is especially important when modeling the transient glacialinterglacial cycle and when quantitatively comparing the Eemian and the present day.Our analysis shows that this 435 effect is indeed very important, both for equilibrium estimates of surface mass balance and especially for simulating the transient evolution of the ice sheet.
If the correction term presented here is incorporated into reduced complexity models, such as PDD, it should then be 440 possible to tune the insolation-independent parameters from present-day data, while relying on the insolation term to handle changes in paleo conditions.Such an approach could help to better constrain the range of melt model parameters to realistic values.Incorporating the effect of insolation changes 445 into a PDD approach would however increase the number of parameters used in the model.This implies that a more direct approach could be to simply calculate insolation changes within the model and implement a simple representation of albedo.

Conclusions
We have quantified the relative importance of insolation changes to surface melt and its evolution during Eemian and Holocene conditions using a simplified energy-balance model.Its impact is highest for peak Eemian and Holocene conditions and decreases with decreasing insolation.Our results are consistent with those of Robinson et al. (2011) and van de Berg et al. (2011), who found a contribution around 50% due to higher insolation during peak Eemian conditions.However, its relative impact on total melt depends on the temperature anomaly.Locally and for a given insolation, its contribution to melt is relatively more important for colder surface conditions.
Additionally, we find that the contribution of insolation to melt is non-linearly dependent on temperature, as feedbacks between temperature, melting and albedo modify the ice sheet surface conditions.Spatially, the relative importance of insolation is largest in the cold accumulation zone and decreases towards the ice sheet margin, where the effect of temperature is dominant.The largest impact on the snowpack can be seen near the equilibrium line, as this is where an additional melt contribution can most strongly modify the surface conditions through the albedo-melt feedback.In simulations where the surface mass balance is coupled with a dynamic ice sheet model, the additional melt due to insolation and associated feedbacks can cause a significantly larger reduction in ice sheet volume during the Eemian.This highlights the need to account for insolation changes in paleo ice sheet modeling.
Based on this analysis, we have presented a simple parameterization that can be used to correct for the effect of insolation changes in models that do not explicitly account for them.The provided parameterization assumes an additional melt contribution due to local insolation anomalies as a function of temperature and season.It also implies that a given melt model will not be modified for present-day insolation, and thus consistency can be maintained with previous work.
The new parameterization has been shown to give good results when used to correct the melt calculations of a fully coupled ice-climate model simulation that does not explicitly account for insolation changes.
) Fig. 3. Monthly diagnostic values of the coefficient a versus temperature (shading represents density of points), along with the parameterization given by Eq. 6 (magenta line).Month Tmin (°C) q q q q q q q q q q q q Fig. 4. Seasonal value of Tmin estimated from a linear fit to the diagnostic values of a in the range of 0 < a < amax (black points) and compared to the parameterization given by Eq. 7 (magenta line) A.
and H. Goelzer: Insolation in paleo ice sheet modeling This formulation follows from the non-linear relationship of local conditions to temperature during different seasons.Figure 4 shows how T min changes over the year according to Eq. 6.This is compared with the T -intercept of a linear fit to each month of values in the range 0 < a < a max , which indicates that our parameterization reproduces the seasonal change in the value of T min well.The values of a max , T max ,

Fig. 5 .Fig. 6 .
Fig.5.Greenland ice sheet melt anomalies for the Eemian (left) and the Holocene (right) given a 3 • C summer temperature anomaly relative to today, calculated with the full model (M (S), thick black line), the reduced model (M (S0) thin black line) and an ensemble of insolationcorrected model versions (M (S0)corr) with the parameter Tmin set to -15, -14, -13 and -12 • C (cyan to magenta lines, respectively).

Fig. 7 .
Fig. 7. Transient simulations of the Eemian (left) and the Holocene (right) with a fully coupled climate-ice sheet model.Top: Applied summer temperature anomaly forcing.Bottom: Ice volume from simulations in which melt is calculated with the full model (thick black line), the reduced model forced with present-day insolation (thin black line) and an ensemble of insolation-corrected model versions with the parameter Tmin set to -15, -14, -13 and -12 • C (cyan to magenta lines, respectively).

Fig. 8 .
Fig. 8. Distribution and timing of the minimum ice sheet volume during the Eemian for coupled simulations using the full model (left), the reduced model (middle) and the insolation-corrected model (right).
compares the performance of the reduced melt 305 model with that of the corrected model for peak Eemian insolation.While the reduced model is able to reproduce the full model results in the interior of the ice sheet (where melt is essentially zero), it underestimates the annual melt by up to 400 mm further out towards the margins.The maximum