Quantitative and qualitative constraints on hind-casting the formation of multiyear lake-ice covers at Lake El’gygytgyn

Analysis of the 3.6 Ma, 318 m long sediment core from Lake El’gygytgyn suggests that the lake was covered by ice for millennia at a time for much of its history and therefore this paper uses a suite of existing, simple, empirical degree-day models of lake-ice growth and decay to place quantitative constraints on air temperatures needed to maintain a permanent ice cover on the lake. We also provide an overview of the modern climatological and physical processes that relate to lake-ice growth and decay as a basis for evaluating past climate and environmental conditions. Our modeling results indicate that modern annual mean air temperature would only have to be reduced by 3.3 C± 0.9C to initiate a multiyear ice cover and a temperature reduction of at least 5.5C± 1.0C is likely needed to completely eliminate direct air–water exchange of oxygen, conditions that have been inferred at Lake El’gygytgyn from the analysis of sediment cores. Once formed, a temperature reduction of only 1–3C relative to modern may be all that is required to maintain multiyear ice. We also found that formation of multiyear ice covers requires that positive degree days are reduced by about half the modern mean, from about +608 to +322. A multiyear ice cover can persist even with summer temperatures sufficient for a two-month long thawing period, including a month above +4C. Thus, it is likely that many summer biological processes and some lake-water warming and mixing may still occur beneath multiyear ice-covers even if air–water exchange of oxygen is severely restricted.


Introduction
Sediment cores from Lake El'gygytgyn in northeastern Siberia (67.29 • N, 172.5 • E; Fig. 1) provide the only currently available, continuous record of Quaternary climate and environmental change from the terrestrial arctic.Records from the lake inform our overall understanding of Quaternary climate change in the Arctic and provide a critical terrestrial record against which more commonly available marine records can be examined (Melles et al., 2012).Analysis of the sedimentary record thus far has identified evidence of significant climate and environmental variability and indicated the occurrence of extended periods of anoxia within Lake El'gygytgyn in the past (Melles et al., 2007).Through its fundamental influence on air-water exchange of oxygen, heat fluxes, mixing, and biogeochemical processes, variations in lake-ice cover are considered a dominant control on lakewater oxygenation.Consequently, understanding the sensitivity of lake ice to climatic conditions is fundamental to the interpretation of sedimentary evidence of lake-water anoxia in terms of past climate variability.Nolan et al. (2003), analyzed modern lake-ice cover in detail using four years of synthetic aperture radar (SAR) and Landsat remote sensing data in combination with a physically-based numerical model of lake-ice growth and decay.Modeled ice growth and decay closely matched the timing observed in the remote sensing data.This model included a surface-energy balance model that determined the surface temperature and energy available for freezing, a lakemixing model to determine water temperature and stratification, a snow-transport model to determine snowpack evolution, and an ice growth model that tracked both lake ice and snow ice (formed by the freezing of water saturated snow caused by snow overburden) (Liston and Hall, 1995a,b;Liston and Sturm, 1998).While this model successfully reproduced modern ice growth and decay processes, the complexity of the model and the large number of input variables needed to run the model make it poorly suited to the evaluation of past changes in ice cover at times when climatic conditions are less well constrained.Consequently, this study aims to develop the capability to evaluate past changes in ice cover with simpler models that utilize the most reliable and readily available paleoclimate data, namely air temperature, to help constrain the climatic conditions needed to maintain a permanent ice cover on Lake El'gygytgyn.We also discuss implications of multiyear ice covers on light levels in the lake water, the duration of open water conditions, and the possibility of the lake drying out completely.
Our main goals in this paper are to place some quantitative constraints on the conditions necessary to form and maintain a multiyear lake-ice cover at Lake El'gygytgyn, to provide simple lake-ice modeling tools and insights for the paleoclimatologists working there to use in climate reconstructions, and to use our modern data to provide a baseline against which to compare paleoclimate reconstructions.The 175 m deep, 12 km wide lake lies inside a 3.6 million year old meteorite impact crater that is 18 km wide and has presumably been filling with sediments since the time of the impact.The impact, the crater, the sediment core, and the physical environment are all well-described in a series of prior papers (Belya and Chershnev, 1993;Glushkova, 1993;Layer, 2000;Nowaczyk et al., 2002;Nolan et al., 2003;Brigham-Grette et al., 2007;Cherapanova et al., 2007;Glushkova and Smirnov, 2007;Melles et al., 2007;Minyuk et al., 2007;Nolan and Brigham-Grette, 2007;Brigham-Grette, 2009;Swann et al., 2010) and those found in this special issue (Melles et al., 2012;Minyuk et al., 2013).
We begin by describing the ice growth, melt, and sublimation processes in Sect. 2 and presenting several simple models that were previously developed and tested by others.In Sect. 3 we describe the input and validation data, which includes 7 yr of our local automated weather station (AWS) measurements, 48 yr of NCEP/NCAR reanalysis (NNR; Kalnay et al., 1996) data, several years of satellite remote sensing, and one season of ice growth measurements.Unfortunately, the AWS measurements and ice thickness measurements do not overlap, and we discuss the uncertainty and constraints imposed by this in Sect. 4. Section 4 also describes our selection of model parameters specific to the lake and a sensitivity analysis of them.Section 5 implements the models to learn more about lake ice dynamics at Lake El'gygytgyn and we conclude with Sect.6's discussion of the relevancy of these results to the overall paleoclimate reconstruction effort.

Methods
This study employs simple, empirical degree-day models to evaluate lake-ice growth, ice melt, and ice sublimation.Degree days, divided into positive degree days (PDD) that drive ice melt and negative degree days (NDD) that drive ice growth, simply sum the average daily temperatures above or below 0 • C over the course of a year for PDD or NDD, respectively.For example, a week of weather at +5 • C adds 35 degree days to the annual PDD.Prior research and our own have shown that the differences between degree-day models and more sophisticated ones are small, usually less than 15 cm at the end of winter (Ashton, 1986;USACE, 2002).Because these simplified models track reality so well, especially in static growth conditions like lakes with no river flux, they are a standard tool used in ice engineering.The fundamental physical dynamic that must be captured in any lake-ice growth model is that as the ice gets thicker, it insulates more and slows growth at the lower interface, and this dynamic seems to be captured well by the empirical models used in this study.The surface-energy flux balance is largely what drives ice growth (usually the water temperature at the lower ice interface is at or near zero degrees and stratified, causing latent heat to escape mostly upwards), but because air temperature typically tracks the surface-energy flux balance, it makes a suitable proxy.Similarly, ice melt is also driven by the surface-energy balance.While other processes also come into play in lakes, these processes all tend to follow air temperatures and thus are captured by the degree-day models even though air temperature is not actually driving growth or melt.Because degree-day models are not physically based, they reveal nothing about the physics of melt; for our purposes here, we are only interested in ice growth and melt rates, so this poses no limitations on our goals.
Nonetheless, it is still important to understand the physical processes being modeled to help ascertain the validity of the models in different circumstances in the past, so the sections below include a discussion of these processes.

Ice growth
Ice growth on arctic lakes and rivers can be modeled reasonably successfully using freezing degree days and local tuning factors (USACE, 2004;Arp et al., 2010), in a much simpler way than our prior methods (Nolan et al., 2003).The model we use here (Eq. 1) is a simplified empirical form of the Stefan equation as derived in USACE ( 2002) along with related assumptions and simplifications.In short, it assumes vertical, one-dimensional, quasi-steady heat conduction through a uniform, horizontal ice layer, which is only growing at the ice-water interface.
where NDD are negative degree days (expressed as a positive number) in  2002).β is a newly defined parameter described below.We determine the best α and β for modern Lake El'gygytgyn in Sect. 4. Equation ( 1) can be used to track growth using NDD accumulating as a function of time by using daily accumulated NDD or to simply calculate final ice thickness using the annually accumulated NDD.
We have added a new parameterization to this model, β , to delay the start of freeze-up until a specified number of negative degree days have accumulated.This parameter accounts for the action of wind in breaking up thin ice packs that form at the start of the freezing season and other similar processes.A number of factors lead to the formation of ice, as described nicely by Locke (1990), but typically sub-freezing temperatures, cold skies, and calm winds are required to initiate a stable ice cover.The ice must then grow thick enough to prevent wind action from destroying the ice and starting the process over.For example, we noticed ice skims forming in August 2000 on calm clear nights at Lake El'gygytgyn, but wave action quickly destroyed them the next day.Thus, depending on the windiness of location, the latitude, and other factors, the value of β is site-and year-specific.Our choice of this delay factor is described in Sect. 4 and is driven by our local observations.
The addition of snow slows heat transfer further.Because ice grows fastest when it is thinnest, all else being equal, the timing of snowfall (and its insulating value) can strongly affect ice growth, especially in the early season.However, Eq. (1) does not account for interannual changes in snowpack (as could our prior model, if measurements were available to drive it), but rather parameterizes nearly all local effects in the tuning coefficient α , and therefore all conclusions here are based on snowpack dynamics remaining constant through time.In the next section on weather data, we describe why we believe this is a reasonable assumption in the modern environment, given that wind limits the snow thickness regardless of variations in the amount of snowfall, and include a discussion on the validity of this assumption in the past.

Ice melt modeling
In terms of initiating and maintaining a multiyear ice cover, predicting how ice grows is not as important as predicting how lake ice melts, as we describe later, and this unfortunately is a much more complicated process as it is more site dependent.Our prior remote sensing from 1997-2001 indicates that snowmelt on the lake surface begins in mid-May, a moat forms at the edge of the lake ice in the last week of June, and complete ice melt occurs in mid-July.Once the moat opens up in June, the ice pack is subjected to substantially higher mechanical erosion due to wind shove, and our remote sensing data show (Nolan et al., 2003) that soon after moat formation, the ice pack forms major leads across the lake at deltas or outcrops that jut into the lake (e.g., at Stream 12 on the west side, Nolan and Brigham-Grette, 2007).During this time, about half the melt energy is supplied by solar absorption in the ice pack, and crystal grain boundaries weaken, forming candle ice which lacks the cohesion of the winter ice pack (Locke, 1990).Once separate floes form, mechanical erosion increases further and soon the diminishing ice cover begins sloshing from one side of the lake to the other based on wind direction, causing massive ground-shove features as it becomes grounded (e.g., Fig. 1), and turning the pack into a jumble of bits rather than cohesive floes.The ice typically disappears completely within a few days of separate floes forming through wind action, which is almost always present.
Thus, there are at least three important variables in predicting ice melt in summer here: simple ice melt from summer sun and warm air, moat formation, and wind action.Unfortunately, physically modeling moat formation requires a full 3-D model with more input data than we have.The modern moat is likely formed by a combination of snow meltwater accumulation from the ice surface and surrounding landscape plus, and perhaps more importantly, from solar heating of adjacent land surfaces and shallow water shelves (< 5 m deep), which warm considerably from solar inputs through the snow-free lake ice.A thermistor string we deployed on a shelf here in 2000 (Nolan et al., 2003) showed that the water can reach 5 • C in spring.The shelves are likely former beaches, formed when the lake level was lower (Juschus et al., 2011).If lake water levels were 20 m lower than present, the area of lake bed that could absorb significant solar radiation (the littoral zone) would be significantly reduced, as it would be confined to a narrow strip of the steeply sloping lakebed (10-40 found today in the substantial lagoons and shelves would likely not occur under conditions of lower lake level.Thus in glacial times, when precipitation inputs were likely much lower due to the global ice sheet dominance, moats might have been much smaller and harder to form.Without moats, mechanical erosion caused by wind would also have to have been smaller, and the effect of wind itself would have decreased as the lake surface would be more deeply recessed into the lake basin.Again, these are affects beyond the scope of this paper to try to quantify, but for these reasons we believe it reasonable to assume that lower lake levels would restrict seasonal decay of lake ice, all else being equal. Given such large uncertainties in the past and complexities in the current environment, we decided to use another simple degree-day melt model to predict ice melt rates and potential.The physically-based 1-D energy balance model used previously (Nolan et al., 2003) to predict ice melt rates for the winter of 1999-2000 using NNR driving fields found that melting a 1.8 m ice pack takes about 4 weeks by simple melt alone, starting when air temperature rose above freezing and ending at about the summer peak in temperatures (12-15 • C).This model also ignores all 2-D and 3-D effects, such as melting at grain boundaries, causing candle ice formation and the influence of moats and wind.However, the predicted ice-out date due to 1-D melt alone was within a few days of actual, so we conclude that degree days are a reasonable proxy for the combined actions of all influences.In this paper we use positive degree days corrected by a local tuning factor following (Bilello, 1964): where Z max is the end of winter maximum ice thickness in centimeters, α is a tuning factor, PDD are positive degree days, and β is a delay factor similar to β that prevents ice from forming until sufficient PDD has accumulated to prevent brief late-winter thaws from melting ice or throwing off the date of melt onset.We determine the local value for α and β in Sect. 4. Equation (2) can be used to track daily ice thickness using daily accumulated PDD or, by letting Z melt go negative, Eq. ( 2) will calculate the melt potential using annually accumulated PDD.That is, melt potential is the ice thickness that could be melted if it were present at the start of winter.Similar degree-day models have been used with great success to model glacier ice melt (Hock, 2003(Hock, , 2005;;Hock and Holmgren, 2005).In Sect.4, we discuss a Z min that sets a minimum ice "thickness" that parameterizes how much ice needs to be left over at the start of winter to consider an ice pack as perennial.Note that while Eq.(2) calculates a thickness, Z min , we are not using this value as a physical thickness, just a parameterization of all of the processes that lead to ice breakup or maintenance into winter.That is, even though mechanical breakup might occur when the ice is still 50 cm thick, this parameterization might indicate 0.2 or 1.5 m.Given our satellite and field observations of lake ice, breakup usually occurs within several weeks of moat formation, as this gives the ice room for wind shove to move it and mechanically break it up, as well as time for weakening of the crystal boundaries and thus ice pack competency.Therefore choice of this parameter is guided by local, site-specific observations and intuition.

Ice sublimation models
Once a multiyear ice cover forms in glacial times, could ice sublimation rates exceed precipitation or groundwater inputs and dry out the lake?Unfortunately we do not have enough information to attempt to model the answer to this question.So rather than present a coded model here, in this section we provide some discussion on rates and controls of sublimation that paleoclimatologists can use to constrain sublimation rates when interpreting core proxies and processes related to multiyear ice covers.
The most widely studied modern examples of permanently ice-covered lakes are those in the Dry Valleys of Antarctica.Most of the sublimation there occurs during summer when relatively warm, dry winds pass over the lake surface.Ice sublimation largely balances any ice growth, maintaining a fairly constant ice thickness.The ablation rates of these lakes range from 10 cm to 1.5 m per year (including sublimation and evaporation of landlocked melt) with 35 cm per year being typical (Clow et al., 1988;Squyres et al., 1991;Doran et al., 1994Doran et al., , 1996)), noting that the occurrence of warm katabatic winds from the nearby ice sheet are in part what is responsible for the high rates of sublimation here.As a first approximation, therefore, it seems reasonable to assume that sublimation rates of multiyear ice at Lake El'gygytgyn would be within this wide range.These estimates of sublimation are within the the range of modern precipitation at the lake, thus if some glacial periods were characterized by less precipitation as proxies have indicated (Melles et al., 2007(Melles et al., , 2012)), it is within reason to suspect that sublimation rates may have exceeded precipitation rates.If sublimation exceeded precipitation by 1 to 50 cm, which is within reason, the 175 m deep lake would have dried out in 17 500 to 350 yr, not accounting for groundwater inputs.Given that glacial epochs often last 50 000 yr or more, it seems that there may have been multiple opportunities for complete drying out of the lake.However, given the variability in precipitation, temperature, and ice growth, it is also conceivable that a few centuries or even decades of warmer/wetter weather could balance out several thousand years of water loss, and there seems to be little bathymetric or stratigraphic evidence for submerged paleo-shorelines beyond the current shelves.Further, because the lake's watershed is 2.6 times the area of the lake, land runoff into the lake requires lake-ice sublimation rates to be that much higher than precipitation rates for drying to occur.Completely unaccounted for is the possibility of groundwater infiltration from beneath the thick permafrost layer through the talik that likely exists beneath the lake, as has been suggested previously (Nolan et al., 2003), or its variability over 3.6 M years.Therefore, while we do not have enough information to model sublimation beyond these considerations or determine whether the lake dried out and caused temporal gaps in the proxy record, we can state that this is within the realm of possibility and caution paleoclimatologists to be alert to this possibility in their interpretations and usage of the degree-day models.

Data: input and validation measurements for the models
Here we present our local weather, ice measurements, and the NNR data that we use to select our local empirical tuning factors and to help define the "modern" period for comparison with the paleo-regime.It may be important to keep in mind that our goal in this paper is not to validate these models, as this has already been done in numerous other studies (Bilello, 1964;Ashton, 1986;Hock, 1999Hock, , 2005;;USACE, 2002USACE, , 2004;;Arp et al., 2010), but rather to validate our choice of local empirical tuning factors.This validation includes sensitivity studies to determine the robustness of these models for our applications.Consider, for example, an N-th order polynomial that can be made to exactly match N validation points, but that this polynomial function will have no predictive capability in other circumstances, whereas we demonstrate that any choice of tuning factors within a wide and reasonable range in these 1st order models will predict ice onset and breakup dates within a few days or ice thicknesses within a few 10s of centimeters, which is adequate for our purposes.

Local weather measurements
We installed a weather station near the outlet of the lake in July 2001 (67.446517 • N, 172.186200 • E), as shown in Fig. 1.The station utilized a Campbell Scientific Inc. (CSI) CR10x datalogger, powered by a 100 Ah deep cycle battery recharged by a small solar panel.Measurements were logged hourly, with air temperatures and wind speeds being measured every minute and averaged.There was no telemetry and all data were stored in the logger and backup memory card.The instruments were mostly all standard, sold through CSI: -two Vaisala HMP45c measured air temperature and relative humidity, at 1.0 and 3.0 m above ground surface.
-a CSI SR50 sonic ranger was used to measure snow depth, recording values subtracting its own height so that the measurement read height above ground.
-a Young AQ Wind Monitor at 4 m.
-six thermistors on a calibrated string installed in a pit to depth of 60 cm.
-eight CSI SM615 soil moisture probes installed in a pit to depth of 60 cm.
-a tipping bucket (Texas Instruments) mounted on the ground surface (0.01 tips).
The station had little maintenance from the time it was installed in 2000 to when it was dismantled in March 2009.The day after leaving the site in 2000, the station was shot with a rifle by a local drunk, destroying the logger and the battery.In July 2001, the logger and battery were replaced.
In July 2003, the station was downloaded and one loose wire in the 3m air temperature sensor tightened and level checked on appropriate instruments.Battery power began to fail in late winter 2008, causing some gaps, but functioned fine once sun returned to this Arctic site, until it again failed in early 2009 during the deep drilling project.It was still recording data when dismantled, though the battery failure earlier that spring had caused the date and time to be reset to a random setting.Other than these gaps, the remaining years had at most a few missing time steps caused by static, low power, or some undetermined cause.Calendar years 2002-2007 therefore had better than 99.99 % recovery, and 2008 92.8 %, in terms of hourly time-steps logged.Individual sensors as well as the logging system itself did not necessarily fare well over the recording period.Fortunately, the data from 2001-2003 were described already in detail in Nolan and Brigham-Grette (2007).This paper utilizes air temperature and precipitation data, so we describe those records further in this section, and utilize the data in the next section.
The soil moisture and soil temperature data from this station were used in Federov et al. (2012) to better understand the local water balance.

Air temperature and relative humidity (AT/RH)
The 1 m AT/RH sensor performed without any issues for the duration of the installation.The 3 m AT/RH failed from 2002 Day 206 at 08:00 LT (local time) apparently caused by the socket handle coming loose.This was retightened on 2003 Day 159 at 1800.The only other gaps were caused by logger power-failure in early 2008.None of the probes were replaced or recalibrated during the study period.Comparison plots of the two AT sensors in winters (that is, when there was no solar radiation loading) shows they were always within 0.1 • C of each other with a cross-correlation value of over 99 %, which likely means that drift was minimal or both drifted the same amount (likely the former).Summer values (that is, when solar loading was appreciable) did show differences of up to several degrees, but much of this difference is likely real and the remainder due to the differences in solar loading of the radiation shields and differences in wind cooling rates.The lowest recorded temperatures are −40 • C due to the limitation of the sensor itself, but inspection of the hourly data reveals that such temperatures were uncommon and thus clipping was minimal and did not affect annual means.The maximum recorded range for both RH sensors decreased with time (that is, they no longer reached 100 %) as is typical for these sensors.The range of values is within reasonable expectations, with the largest variability in summer when higher air temperatures substantially affect the absolute amount of moisture the air can carry.

Tipping bucket and sonic ranger
The tipping bucket seemed to remain functional until June 2008.Values from 2008 are suspect, as they are anomalously low, with no tips shown throughout July and August even though the soil moisture probes were clearly recording rainfall during this time.The unit was not re-leveled during the study period, so some bias or drift may have occurred, but it was still firmly anchored upon removal and no noticeable tilt was recorded at that time.Figure 2a shows annual cumulative rainfall for 2002-2007, revealing nearly a factor of 3 interannual difference (73 to 200 mm), though the range is consistent with measurements elsewhere in the Arctic (Kane and Yang, 2004).The tipping bucket had no wind shielding, thus these values are probably low by 10-50 %, possibly more (Yang et al., 2005).The sonic ranger began performing erratically after 2002, but still provided useful (though noisy) information through 2005, with maximums depths of about 50 cm.After this it appears to be not useful at all.Likely the shielding deteriorated over time and was interfering with the sensor or the cable or its contacts may also have degraded.
In any case, though there is data there that may be useful to other studies, we decided not to use the sonic ranger data in this study as the only high quality data, from 2001-2002, was reviewed thoroughly in Nolan et al. (2003) who found an end of winter snowpack thickness of 40 cm and water equivalent of 11.0 cm.

Local ice thickness measurements
The only systematic measurements we have of ice onset, decay, and moat formation are from our prior work (Nolan et al., 2003) using satellite remote sensing to develop average dates of onset and breakup of lake ice from 1998-2001.These breakup data are not thickness, but rather percentage of ice cover, so the only time we know the actual ice thickness is when it reaches zero and we know this to within about a week in these years.These data are used primarily to tune the model melt parameter α and freeze-up delay β later in Sect. 4. The only measurements of ice thickness that we have at Lake El'gygytgyn occurred in June 2003 (Melles et al., 2005) and during the lead-up to the major drilling effort that occurred in spring 2009 (EBA, 2009) to recover the deep sediment cores.Snowmelt on the lake ice began on about 18 May The NNR fails to capture daily events, but does better in capturing annual trends, though 2-3 times higher in magnitude.In (b), annual precipitation from NNR, showing that winter and summer magnitudes are roughly equal, and dry in general.and ended about 8 June in 2003, at which point the moat began forming around the lake.In late May, ice thickness of 190 and 200 cm were measured near the center of the lake.From 10-14 June, a few days after snowmelt had completed, a north-south transect was made across the 12 km lake, measuring thickness at 12 locations.Mean thickness at these locations was 166 cm with a standard deviation of 17 cm, indicating substantial variability in melt rates, likely due to quite localized dynamics of snowmelt, albedo, and water mixing.The moat reached 10 and 30 m wide by 17 June and 8 July, respectively, and the lake was clear of ice by 19 July.No measurements were made of ice formation the previous fall, but visual inspection of the ice pack clearly indicated that flows had formed, been pushed around and refrozen before the final ice pack had formed.In 2009, ice roads were made, drilling pads thickened, and ice conditions were monitored for safety (EBA, 2009).In these engineering efforts and plans for drilling, much use was made of our prior work's ice growth curve and the thickness measurements followed that curve within several centimeters for most of the drilling effort (EBA, 2009).Note that measurements were made at various parts of the lake but we use only measurements from those locations not affected by the ice road or pad thickening.The drilling stopped in late-April and no further ice measurements were made, such as peak thickness or melt dates.These data are used to tune the model to growth rate parameter α in Eq. ( 1) later in Sect. 4.

NNR compared to AWS and proxy for modern weather
Our goal in this section is to determine how well NNR (Kalnay et al., 1996) captures our AWS air temperatures because (1) our AWS data does not overlap with our ice thickness validation data so we later use NNR to drive our models, (2) paleodata input to ice models will be derived from models similar to NNR, and (3) we will use a subset of the NNR over the past 50 yr to characterize the "modern" period against which we can make paleoclimate comparisons.We extracted the NNR parameters for the grid node encompassing the lake from 1961-2009 encompassing the highest quality time-range of this record.Of primary interest to the lake-ice modeling here are air temperatures and precipitation fields.We are interested in NNR air temperatures primarily because air temperature drives our models.We are interested in NNR precipitation because snowpack affects ice growth and because rainfall affects recharge, especially when sublimation of multiyear ice packs is occurring, so understanding the modern range of variability will help us put limits on modern ice dynamics for later comparisons with the paleoenvironment.

Air temperature comparisons
Air temperatures compared well between NNR and AWS.Mean annual air temperature (MAAT) and its trend compare well between datasets (Fig. 3a) -especially considering the coarseness of the reanalysis grid (about 200 km), with the AWS recording 0.1 to 2.0 • C colder than NNR, with a mean difference of 1.1 • C from 2002-2007.Also shown in Fig. 3a are PDD and NDD for both, which indicates that winters are largely responsible for the differences in MAAT. Figure 3b compares average daily temperatures for 2003, a year with one of the largest MAAT and NDD misfits.Though there are minor differences in values, the daily trends are all captured well, with a cross-correlation coefficient of 0.96, which was typical for all years.The diurnal fluctuations revealed by the hourly data (not shown) suggest that there are swings in temperature throughout the day that are larger than the mean misfit.It is difficult to say why the match is better in summer, but likely much of the difference comes from topographic effects in the crater bowl that is not captured by coarse NNR grid and therefore it cannot reproduce cold-trap inversions as well.The most important point here is that daily and annual trends are captured by the global model, and the offset, though variable, has a mean of 1.1 • C; later we show that these differences have only a small effect on ice dynamics.Because these lake-ice models are not physically based, their best use is as a comparison between two different climatic regimes, in this case between the modern air temperatures and air temperatures of some period during the past 3.6 Ma.Though there are several options to use to represent modern conditions, a multiyear mean is satisfactory because it is easy to calculate, it naturally smooths interdaily variability, it typically creates single transitions across 0 • C, and different time intervals are easily tested.Because there was a strong trend in modern NNR air temperatures, we explore the trend in more detail here to determine a suitable interval to use as the modern era. Figure 3a indicates a strong warming trend in NNR MAAT beginning about 1988. From 1960to 1988, the mean NNR MAAT was −12.1 • C with a standard deviation of 0.79 • C and a slight cooling trend.From 1989 to 2009, the mean MAAT was −10.1 • C, more than 2 standard deviations higher than the mean of the previous period.Prior research here (Nolan and Brigham-Grette, 2007) noted the start of this trend (up to 2002) and attributed it to winter warming, with negative degree days (NDD) rising more than the magnitude of the long-term mean of positive degree days (PDD) and a sharp reduction in the number of days below −30 • C. That is, all of this warming could be explained by changes in winter temperatures.We can further confirm the general conclusion that winters show more variability now using our AWS data from 2002-2007, which shows a great correspondence with NNR in PDD in terms of magnitude and variability and a similarly large variation in NDD (Fig. 3a).However, the winter magnitudes themselves do not match as closely between datasets as summer ones, and our AWS record is not long enough to confirm the timing and magnitude of the winter warming indicated by the reanalysis data; the possibility remains that changes in input datasets to the NNR could cause this trend spuriously, though the most notable changes in input data occurred a decade earlier with the introduction of satellite-derived fields.In any case, we decided to use the range of 1961-1988 as representing what appears to be a relatively stable period of modern climate, a long-enough period to be considered climate, and seemingly unaffected by trends (real or spurious).The mean from 1961-1988 is shown in Fig. 4, along with the minimum and maximum daily temperatures.
When comparing to the "modern" period, it is also important to realize that this mean is likely about 1 • C warmer than actual (with biggest misfit in winter), with the bulk of the error likely caused by the model not capturing low-level inversions in winter.Given the recent warming, it appears that changes in NDD are driving changes in MAAT, but, as we show in the next section, it is really PDD that controls the existence of multiyear ice.So to the extent that modern conditions represent the past, it is important to note that a decrease in MAAT is likely primarily caused by a decrease in NDD and this may not influence multiyear lake-ice formation unless there is a corresponding decrease in PDD.In a companion paper in this issue, we explore this winter-warming trend further to determine that it is mostly occurring in subfreezing temperatures in spring and fall, with more summerlike weather patterns dominating in recent years during those seasons (Nolan et al., 2013).In this paper, of most relevance is that summer temperatures show little change in mean over the past 50 yr, that the NNR trend captures summer magnitude well, that summers are most important in controlling the initiation of multiyear lake-ice covers (as described in Sect.4), and that the warming of the past 20 yr is postmodern in the context of this paper.This mean is what we used for the "modern" era in comparisons with paleoconditions.

Precipitation comparisons
Precipitation between the two datasets did not compare as well as temperature in terms of magnitude, though the annual trends seem to be captured.We focused our precipitation analysis on rainfall because our local snowfall measurements were not suitable for long-term analysis, and because in this windy arctic environment, wind scour is the primary control on snow depth as it is in all other windy arctic environments (Liston and Sturm, 1998;Sturm and Liston, 2003;Sturm and Benson, 2004).Figure 2a compares the annual cumulative rainfall.Here we started the NNR cumulative calculation on day 150 to compare with the tipping bucket, as the NNR field does not distinguish between rain and snow; as noted in Nolan and Brigham-Grette, (2007), likely some of the summer precipitation captured by the tipping bucket was snow that landed in the funnel and subsequently melted, but we have not corrected for this here.As can be seen, the NNR rainfall can be four times that of the AWS (Table 1).Though NNR distinguishes high rainfall years from low rainfall years, it was not consistent in picking up intermediate years well nor was it consistent in picking up magnitudes of individual events.NNR was also not a good predictor of which days rainfall occurred.For example, in 2002, there were 126 days in summer when either the AWS or NNR indicated rainfall, but only 37 of those days were the same.Though it is certain our rain measurements were lower than actual, the likely range is only 10-50 % undercatch (Yang et al., 2005), not the 2-4 times difference seen here.Thus, use of the rainfall field of NNR as a proxy for local weather data is only done with errors on this level, as is the use of it to determine which weather patterns bring wetter or drier air.
Because the parameters in our models incorporate the influences of winter snowpacks, we also explored the use of NNR data to better define modern conditions against which paleo-precipitation could be compared (Melles et al., 2007,  2012).In addition to direct impacts on lake ice, snowpacks also relate to levels of sunlight reaching the lake water thereby inhibiting photosynthetic activity and lake water circulation, and rainfall places limits on whether the lake may have dried out due to winter sublimation exceeding inputs.The NNR data shows an increasing trend in precipitation totals over the past 20 yr (Fig. 2b), similar to the increasing air temperature trend.As can be seen, total precipitation rose about 50 %, from about 35 to about 55 cm a −1 .Unlike air temperatures, however, there is no clear-cut seasonal explanation -both summer and winter precipitation rose in roughly equal amounts with roughly the same interannual trends.Here we have defined the summer period as Day 150-275 based on our AWS tipping bucket, as previously described.The NNR data show that summer precipitation is slightly higher in annual percentage (57 %) than winter.However, based on our previous analysis we know that NNR is overestimating summer rain by 2-4 times measured, but winter NNR values are approximately what we have measured as end of winter snow water equivalent in 1998, 2003, and 2009, about 10-15 cm, and what we measured from lake pressure changes (15.0 cm) in 2002-2003 (Nolan and Brigham-Grette, 2007).So it is likely that winter and summer precipitation rates are roughly equal, or that winter may be a bit higher, but we do not have enough information to support this further.Also, our local measurements in 2001-2002 indicated a total annual precipitation of about 20 cm, which is roughly half of what NNR indicates, though we indicated that our measurements were probably low.In any case, reasonable values for "modern" precipitation are on the order of low decimeters and snowpack thickness seems to be controlled as much by wind scour and sublimation as it is by precipitation, since measured end-of-winter snowpacks range from 0 to 50 cm here but never higher except in drifts.Thus, it is likely good practice to consider that any analysis of precipitation data from paleomodels is not likely to be more accurate than our comparison of modern reanalysis to observed precipitation data and that these paleomodels are probably best used in terms of relative comparisons to modern models rather than actual magnitudes.

Model implementation: parameter selection and sensitivity testing
In this section we describe how we selected the model parameters based on our calibration data, and in the next section on results we describe the sensitivity tests we used and how our uncertainties would affect our results.The results are shown in Table 2.

Ice growth parameters
To model ice growth, we found setting α to 2.45 in Eq. ( 1) matched our calibration data bset (Fig. 5b).Our choice of α was optimized to meet the measured growth rates in 2008-2009 using the NNR air temperatures from that winter (since our AWS had failed that spring), with a resulting mean and RMS difference between the model and our 7 thickness measurements of 1.0 and 2.7 cm, respectively.As can be seen in Fig. 5b, the results compare well visually with field data.Unfortunately we did not measure maximum ice thickness in this year.Using this value of α applied to 1999-2000 yields a maximum ice thickness of 150 cm, or 30 cm less than predicted by our prior model (which has no thickness validation of its own).Similarly, these values predict a maximum ice thickness of 1.6 m in 2003, when we have measurements indicating 1.9 m.Setting α = 1.85 did match maximum observed 2003 thicknesses but created a large misfit in growth rates from 2008-2009, but these growth rates may not be applicable in 2002-2003 since we do not know when the ice formed or its rate of growth during that winter.No amount of parameter tuning could match the growth rates from 2008-2009 and simultaneously yield a maximum ice thickness of 1.8 m or more, but we do not have enough information to know if this is a misfit or not since we do not have enough field data from either year to determine this.Numerically the issue is the use of the square root in Eq. (1) as approximation of the Stefan solution likely begins to break down when NDD gets large, and most prior studies were not dealing with such large NDD.Therefore, we also experimented with varying the power of the root (0.5) in Eq. ( 1) and adding a multiplier to the NDD taken to that power, but these also failed to match both the growth rates and maximum thickness.A different  equation with more tuning parameters likely could represent these larger growth rates, but, as we show later, multiyear ice growth is much more sensitive to summer melt and the shortcomings of Eq. ( 1) do not affect our general conclusions, as using either value (2.45 or 2.85) does not change the ice-out date by more than 3 days since the ice melts at almost 7 cm day −1 .We selected 2.45 because it provides the best match to the available data.A variation of ±0.4 about 2.45 yields a difference in maximum thickness of ±25 cm, and given that this captures the range of our field data we believe it is a reasonable range of uncertainty, noting that use of 2.45 for modern era comparisons may be under-predicting maximum ice thicknesses by 20-30 cm.
Our choice of 340 degree days for the freeze-up delay parameter β was based on several things but relied most heavily on our remote sensing data.Observations from 1997-2000 indicate that 20 October (Day 293) is a typical date of the onset of freezing, even though the 2008 freeze-up occurred a month later.By mid-October the average air temperatures are usually sustained at −10 • C or below.At these temperatures ice is growing at over 5 cm per day once nucleated and wind can no longer mechanically break the pack completely after a day or two of such growth, as leads quickly skim over in the cold temperatures.The difference in freeze-up date between our 2001-2002 model curve (Day 255 in Fig. 5a) and reality is unknown, but the effect of wind is definitely to delay the onset of growth, as seen in 2008-2009 (Day 310 in Fig. 5b).The 2009 measurements were made in what was likely one of the latest freeze-ups in the modern era, likely caused by a storm in late-October to early-November that brought temperatures near freezing and strong winds (V.Neth, personal communication, 2011).By the time 650 negative degree days accumulate in all NNR years, air temperatures are sustained at −15 • C or colder, ice forms at rates well over 5 cm day −1 , and no significant warm spells occur after this date; thus this is likely the maximum of the possible range.Based on our remote-sensing date of 20 October, the mean NNR air temperature we use to define the modern era leads to a β = 340 • days.Given the large fetch of the lake, our modern data indicates that if paleo-wind was similar to today then air temperatures would need to be sustained below −10 • C and a similar amount of degree days accumulated before ice would form and be sustainable.Thus, the range of 0 to 650 captures the range of possibility without being able to predict it, and we believe using a median value of 340, which best fits our 4 yr of remote-sensing data, is probably the best choice in the absence of other information.
While snow thickness and timing exerts strong controls on ice growth, particularly in the beginning of winter, our lack of snowfall data likely does not negatively impact our modeling.Our remote sensing and field observations have shown that the north end of the lake is blown clear of snow every year in the modern environment, and based on our measurements we believe it likely that snow accumulation in the southern end does not get much deeper than 1 m due to wind scour.Therefore, it is an implicit assumption of all of our conclusions that the snowpack is limited by wind, not by precipitation, and thus there are no significant interannual differences in terms of the affect snowfall variations on ice growth.To the extent that the paleoenvironment was as windy as it is in the modern environment, the assumption will be similarly valid.

Ice melt parameters
To model ice melt, we found setting α to 0.58 matched our 1999-2000 melt curve best (Fig. 5a), described in Sect.3.2.Values for α vary in the literature from 0.2 to 2.0.Our prior work indicated average melt rates at 6.9 cm day −1 and this Clim.Past, 9, 1253-1269, 2013 www.clim-past.net/9/1253/2013/matched to within two days the observed ice-off dates from remote sensing.A value of 0.58 for α matches this average rate exactly using the same NNR driving temperatures used in our prior work.Varying this value by ±0.05 changes the rate from 6.2 to 7.4 cm day −1 and changes the ice-off date by ±2 days, respectively, and we consider this to be the likely range.This range is likely valid under differing climatic conditions, as it is largely dependent on latitude and topography which are essentially fixed over the time period of interest, and even major changes in cloudiness or orbital variations are unlikely to change rates beyond this range.By inspection of each year of NNR daily air temperatures, we found setting β to 10 PDD worked well for eliminating brief springtime warm spells from triggering the model to begin melt.

Results and discussion: modern baseline and paleoconditions
The results of primary interest to us in using these models are maximum ice thickness, length of open water season, and possibility of multiyear ice formation in both the modern and paleoenvironments.In the first section, we explore the sensitivity of our parameter choices on the modern and post-modern era using individual years of NNR data and in the next section we consider the paleoenvironment using the modern mean NNR air temperature as a baseline for comparison.

Modern ice dynamics
We modeled ice growth using the NNR 2 m air temperature field to explore the range of thickness variability and open water season from 1961-2009 and to provide further insights of the sensitivity of our choice of α , α , β , β , as well examine errors due to using NNR as opposed to AWS for driving data.The results from the modern era are summarized in Table 3. Figure 6a indicates that the mean ice thickness (blue lines) and mean ice melt potential (red lines) in the modern era was 163.9 cm ± 8.9 cm and 384.3 ± 70.8 cm, respectively.Because ice melt potential is so much larger than ice growth, there is no possibility of multiyear ice in the modern environment.The thick lines indicate the result for our choice of parameters (α = 2.45, α = 0.58, β = 340 degree days, β = 10 degree days).The thin lines indicate our sensitivity testing, with ±0.4 and ±0.05 for α and α , respectively, each resulting in ±25 cm uncertainty.Even with maximum parameter uncertainty multiyear ice is not a possibility, though 1965 and 1998 stand out as coming the closest.This analysis indicates that the initiation and presence of multiyear lake ice is much more sensitive to summer melt than winter growth.Growth increases as the square root of degree days whereas melt increases linearly with degree days.Thus, modeled melt is more sensitive to both actual summer temperatures and to the choice of α , and this likely reflects reality since actual growth rates decrease with time since it forms on the bottom of the ice where the thickening ice creating more insulation from the atmosphere and the actual decay rates stays fairly constant since it melts primarily at the air-ice surface.
From Fig. 6a it is clear that over the full modern period the range of extremes in ice thickness predicted, 145 to 179 cm (35 cm), is small despite the large magnitude changes apparent in MAAT and NDD, while the range of extremes in melt potential, 214 to 590 cm (76 cm), is much higher despite the magnitude of PDD variability being much lower than that of NDD (Fig. 4). Figure 6b indicates that the modern mean dates of the start of ice melt (red lines), ice-out (black lines), and start of freeze-up (blue lines), dates are days 157 ± 7 days, 197 ± 8 days, and 292 ± 7 days, corresponding to 6 June, 19 October, and 16 July.The dates of ice formation and ice-out match our remote sensing perfectly.We did not observe the start of ice melt in that work, but in 2003, lake-ice melt began on 8 June.Thus, we have strong confidence that this mean is representative of modern conditions for the purposes of paleo-ice modeling.
Thin red, blue, and black lines in Fig. 6b indicate our sensitivity testing for the parameters that most heavily control each result.The start of ice melt (red lines) is most heavily determined by β .Choosing a β value of zero produces substantial variation in starting data of ice melt as the first day above freezing is selected; as can be seen, this has the most effect from 2000 onwards due to warming of winters during this period.Though variations can be as much as two weeks, inspection of the daily data for each year quickly determines obvious errors on the early side caused by spikes of warm weather.However, we know that snowmelt begins in mid-May (Nolan et al., 2003), but it is likely a week or more before the snowpack is isothermal and probably another week before significant ice melt begins.Unfortunately, we do not have any field data on this to tune it further, but filtering out spikes by setting β to at least 10 is visually satisfying when looking at the daily data.The date of ice-out (black lines) turns out to be equally sensitive to the range we consider reasonable for both α and α , 0.4 and 0.05, respectively.Varying α causes the ice to be thicker or thinner, thus taking more or less time to melt; varying α directly controls   the rate of melt.In both cases the variation is ±5 to 7 days.The date of initial freeze-up is most strongly controlled by β , which directly delays this event until a specified number of degree days accumulate.We varied this from 10 to 670 degree days, as indicated by the thin blue lines.Setting β any lower causes any summer cold snap to initiate freezeup; while we have observed ice forming in early September, it is quickly destroyed by even minor wind.At 670 degree days, the full force of winter has set in and there is no possibility of preventing sustained ice formation at this point.Regardless of parameter choices, there seems no possibility of the ice-out date crossing the freeze-up date and leading to a multiyear ice pack.
This analysis allows for us to calculate the open water season in the modern era, which we show in Fig. 6c.Here, the mean open water season is 96 ± 10 days, with an extreme range of 65 to 114 days.Open water season is simply the time between ice-out and freeze-up.Ice-out dates are comparatively much more tightly constrained than freeze-up dates, as indicated in Fig. 6b.Ice-out dates are also comparatively much more stable in terms of the physics involved, since it is largely controlled by predictable solar inputs, whereas freeze-up dates involve wind and synopticscale weather patterns, which vary considerably year-to-year.Thus, the length of open water season is more sensitive to freeze-up date, which the model delays until β degree days have accumulated.
Note that the thin lines in Fig. 6 are not error bounds; they merely indicate how parameter choices affect the results and thus give an indication of uncertainty in choosing those parameters to fit a wide variety of natural variation.That is, our parameter choices were determined from comparisons with a limited set of field measurements that was not coincident with our local weather measurements, so we cannot test against other years' data not used in tuning to calculate errors.In our case, our parameters are tuned to match the field data almost exactly (such that the error bars in doing this are roughly the thickness of the lines on the graph), but these field data are not necessarily representative of all years and may have their own errors and uncertainties.Further, some of these parameters represent phenomena that have huge natural variability, such as the onset of freezeup, which is largely controlled by wind, so we estimate the uncertainty due to wind with the range of β .Another potential source of error here is our use of NNR model data compared to real measurements.We found that the mean differences of the six best years of AWS data from [2002][2003][2004][2005][2006][2007] and NNR led to model results within less than 10 % of each other, with the AWS-derived maximum ice thickness 10 cm thicker, ice melt potential 38 cm thicker, and open water season length of 11 days longer; that is, NNR leads to slightly thinner ice and thus more conservative towards estimating multiyear ice potential.These differences might be important for studying the dynamics of those particular years, but these differences are small compared to the interannual variability and uncertainties in our parameters and thus do not affect our conclusions based on using long-term NNR means.In any case, it may be important to reiterate that our parameter choices are not the last word in what to use.We have tried to explain what physical phenomena these parameters are parameterizing such that others can vary them to suit answering particular questions and encourage users to do so.

Paleo ice dynamics
Much of the existing paleoenvironmental analyses focused on Lake El'gygytgyn have emphasized changes relative to the modern environment (e.g., the "cold-wet", "cold-dry", etc. scenarios of Melles et al., 2007Melles et al., , 2012)).Consequently, we examined changes in ice dynamics caused by shifting the modern mean NNR air temperature record (Fig. 4, black line) toward warmer or colder conditions.We used −9 • C as a lower bound as this eliminates the thaw season completely and thus guarantees multiyear ice and up to +5 • C as this is predicted for the Arctic over the next century (e.g., ACIA, 2005).Shifting the mean in this way is crude and we know from the changes over the past 50 yr seen in Fig. 3 that summer and winter temperatures may vary independently of mean annual conditions.However, our purpose here is to demonstrate how to use the model and do so in a potentially useful way for constraining paleoclimate interpretations derived from the long-term sedimentary record.Even though this method is crude, melt and growth are treated separately by the model, with open water season length most heavily controlled by summer melt, so, using the figures presented here, a reader could estimate for themselves the impact that a colder summer could have if winters stayed the same or changed in a different way without having to run the model.
Figure 7a-c are qualitatively similar to Fig. 6a-c, using the same parameter choices, uncertainty variations, and colors, but substituting a change in temperature for year along the x-axis.Here we have varied the mean NNR air temperature record by uniformly shifting it by −9 • C to +5 • C and running it through the models.As can be seen, with warming temperatures the possibility for multiyear ice decreases substantially as maximum ice thickness changes very little compared to the huge changes in ice melt potential; at +5 • C nearly 7 m of ice could be melted if it were present.With cooling temperatures, the rapidly decreasing ice melt potential severely reduces the open water season until it eventually cannot melt the ice completely by end of summer.The model produces this result when mean annual temperature is 3.3 • C below the modern mean.Noting the thin lines representing parameter uncertainties, this could range from −2.4 to −4.2 • C; given that the degree-day model likely underpredicts maximum ice growth in cold winters, the warmer end of the ±0.9 • C range is more likely than the colder end if it were present, indicating that one exceptionally warm summer in the midst of millennia of colder weather could destroy a multi-year ice pack.The required summer degree days are shown on the upper x-axis in Fig. 7a, indicating that they would change from 608 to 322 degree days with a −3.3 • C shift.Even with such a change to our modern mean, the thaw season would still be about 2 months long with a month at greater than +4 • C, which can be visualized in Fig. 4 by shifting the 0 • C line down by 3.3 • C and examining the mean remaining above this reference line.
Whether a multiyear ice pack leads to anoxic water conditions seen in the sediment proxies is likely controlled by the aerial extent of the ice pack rather than its thickness.That is, just because the ice did not completely melt at the end of summer does not mean there was not significant water-atmospheric exchange, as moat formation and sub-ice mixing can accomplish this (Nolan et al., 2003;Nolan and Brigham-Grette, 2007).Our model predicts only ice thickness, not area, and our remote sensing validation (Nolan et al., 2003) indicates ice area, not thickness, so we have no direct guidance for this or any guidance from the literature that we could find.From our remote sensing we know that mechanical breakup due to wind shove does not begin until large moats form.On 14 June 2003, a week after the first moats formed a few meters wide, the ice was still 1.5 m thick; by 3 July it had been reduced to 50 cm with wind action doing considerable mechanical damage.By 12 July, a week before the ice pack disappeared completely, it was about 20 cm thick.So it seems that 1 to 1.5 m of ice thickness is required to maintain ice integrity in the real world.To limit oxygen exchange, likely at least 90 % of the area must be covered, and to eliminate it probably 100 % must remain.For lack of any further insights or field validation, we used Fig. 5a to estimate the required model ice thickness.On the day in 2000 that ice thickness reached 90 %, there was about 75 cm model ice thickness remaining and 100 % occurred with about 125 cm of model ice thickness or more.Note that this is a parameterized thickness and does not indicate an actual thickness of 75-125 cm, though these are probably close to real-world values.That is, varying this model thickness in essence varies the likelihood of oxygen exchange with the atmosphere, as it parameterizes springtime processes such as moat formation and the extent of sunlight entering the lake due to thinner, smoother ice (as we describe later).
Using these values of 75 to 125 cm in the multiyear ice curve in Fig. 7a (black line) indicates that the temperature would need to be −4.9 to −6.0 • C, or 1.6 to 2.6 • C colder than that required to simply have ice left at the end of summer.The uncertainty in multiyear ice thickness from the thin black lines is about ±1.0 • C. Taking the midpoint of the 75 to 125 cm range and adding the uncertainty yields the result that the "cold" condition of Melles et al. (2007) must be at least 5.5 • C ± 1.0 • C colder than the modern mean air temperature to create a multiyear ice pack that would limit oxygen exchange.

Validity of modern parameter values in the past
The models used in this study are empirical and explain nothing about the actual physics involved with lake-ice growth and decay and as such are tuned to give the best results given particular calibration data.Changes in the amount of solar radiation, cloudiness, and snow cover would require changes in the tuning parameters for maximum accuracy.However, our sensitivity testing and the results of others' use of these models at other lakes has shown that the models do a good job under a wide range of conditions, largely because air temperature is an excellent integrator of the other physical processes involved and because the parameters are more dependent for geography (e.g., latitude, basin shape) than for any particular temperature range.While we have spent some time in this paper trying to achieve the best parameter choices, use of other values within the ranges we have discussed yield comparatively minor changes as a result.Thus, an essential point of this work is not that our parameter choices are best, but that the degree-day approach works well and that future users should feel free to vary the parameters as they see fit to adapt them to paleoconditions, as suggested by comparison of modeled air temperatures to proxy data from sediment cores.For example, if snowfall was thought to be higher, then α could be decreased, or if wind was thought to be less, then B could be decreased.Ultimately, however, there are diminishing returns to tweaking parameter values to suit conditions in the past.The value of these models is really as a first pass to indicate how air temperatures were different in the past compared to today or to predict ice cover given an air temperature.To properly investigate the physical conditions required to explain particular core proxies, the physicallybased model that we used in the past (Nolan et al., 2003) or something similar should be implemented.

Implications for paleoclimate drivers
Our results indicate that there was little possibility for multiyear ice packs over the past 50 yr and strong support for claims that multi-year ice packs existed over much of the 3.6 million year record.To initiate such ice packs, much colder summers than today were necessary, regardless of what happened during winter.Given that glacial periods over the past 3 600 000 yr were often 6 • C colder than today (Petit et al., 1999;Lisiecki and Raymo, 2005;Melles et al., 2012), and lasted much longer than interglacials, we believe this paper lends strong support for the possibility of permanent ice cover persisting at Lake El'gygytgyn for as much as half of this time, perhaps more, since most of this time was characterized by temperatures at least 3.3 • C colder than the modern era.We address the synoptic drivers of air temperature in our companion paper (Nolan et al., 2013) along with some speculations as to which weather patterns likely dominated during the colder/warmer periods.
In any case, we can say with some confidence that it is possible to maintain a multiyear ice cover and still have a wide variety of Arctic summer scenarios.A 4 • C reduction in the NNR mean daily temperatures shows that there would still have been over 2 months of above freezing temperatures and 30 days during the peak with about a 5 • C mean (this can be visualized by shifting the temperature scale in Fig. 4 up by 4 • C and seeing what remains above the new 0 • C).Thus, many summertime biological processes could still have been active and contributing to the proxy record in the sediment core.Mean annual temperatures would have to be reduced by over 9 • C to completely eliminate PDD (assuming a uniform summer-winter change), and based on long-term reconstructions of air temperatures, such cold periods were rare, if ever observed, over the duration of the El'gygytgyn record (Petit et al., 1999;Lisiecki and Raymo, 2005).Therefore, just because the lake ice can persist through summer does not mean that air temperatures never rose above freezing; the question is more about the degree of summer conditions that existed, and we hope that proxies may be exploited that have temperature thresholds, which could further refine our crude modeling.It is again perhaps worth noting that if annual mean temperatures were lower in the past, this is unlikely to have been caused by a uniform decrease in temperature across all seasons.However, Fig. 7 can easily accommodate varying winter and summer differently by applying these shifts to the red and blue curves independently.
Similarly, given that a substantial summer might have existed with perennial ice cover, some warming and mixing of sub-ice lake water likely also occurred.We installed a thermistor string near the center of the lake and later recovered water temperature data from 2000-2003; these data are described in detail in Nolan and Brigham-Grette (2007), and show that warming and mixing began about the same time as snowmelt and prior to the likely start of ice melt, at least in the upper 30 m.We also found evidence for a lowflow, toroidal-shaped current which brought warm, dense water from the shallow shelves to the deepest part of the lake throughout winter, carrying both biota and nutrients (Nolan et al., 2003), and we have no reason to believe that such a flow would not also have existed beneath a perennial ice cover, at least as long as lake levels were high enough to cover the broad shelves where this warmer water was generated by heat transfer from the shallow sediments (Juschus et al., 2011).Thus, is it conceivable that full oxygenation of the lake water could have occurred despite a perennial ice cover, provided moats or cracks could supply part of the lake with that oxygen.

Multiyear ice growth limits
Our main goal in this study was to indicate conditions necessary to initiate multiyear ice growth rather than model ultimate thickness of multiyear ice, as our degree-day models would do a poor job of this without incorporating ice sublimation in some way.Also, as noted previously, the numerical approximations used in the growth model likely begin to underestimate growth after about 140 cm unless new tuning parameters are chosen.All else being equal, if there was a solid ice pack left once air temperatures dropped below freezing in fall, ice growth would begin immediately as wind disruption would be minimized, so the NDD delay factor would be eliminated.In this case, by the end of the following winter, the ice could have been an additional 50 cm thicker than it would have been had growth started in open water, plus the thickness of any remnant summer ice.Thus, after a few years of > 4 • C colder than modern, an ice pack of 3 m or more could have formed, which could have been maintained through subsequent years of only 1 • C colder temperatures based on where the red curve crosses 300 cm in Fig. 7a.However, to keep the ice pack thick enough to prevent oxygen exchange, a 3 • C shift is more likely to have been required.Multiyear ice growth rates would have continued to decrease as the thicker ice provided further insulation, and likely would have reached an equilibrium thickness with sublimation as it does in the Dry Valleys, and thus not continue to increase with time.Similarly, any number of scenarios can be imagined in which a multiyear ice cover can survive an occasional warm year, much in the way that multiyear sea ice reacts.That is, similar to sea ice, multiyear lake ice is harder to form than it is to maintain.However, it is beyond the scope of this paper to test these scenarios, and it is our hope that these simple formulas will be put to use by the paleoclimate community as a starting point for more elaborate scenario testing, leading to the use of a more sophisticated lake ice model to better elucidate the physical processes involved to the benefit of the paleoenvironmental record.

Fig. 1 .
Fig. 1.Automated weather station with wind-shoved lake ice in the background.This station ran essentially unattended from 2002 to 2009.Inset shows location of lake in northeastern Siberia, above the Arctic Circle.

Fig. 2 .
Fig. 2. Precipitation at Lake El'gygytgyn.In (a), cumulative rainfall from the AWS tipping bucket (thick lines) and NNR (thick lines).The NNR fails to capture daily events, but does better in capturing annual trends, though 2-3 times higher in magnitude.In (b), annual precipitation from NNR, showing that winter and summer magnitudes are roughly equal, and dry in general.

Figure 3 .Fig. 3 .
Figure 3. Air temperature at Lake El'gygytgyn.In A, mean annual air temperature and degree-days from NNR (thick lines) and AWS (thin lines), showing that NNR captures the AWS trends well and magnitude best in summer, which is most important to our study.Due to the pronounced warming trend starting in the late-1980s, we used the period 1961-1988 to define the modern era used in this study.In B, daily NNR is compared to local AWS showing that the trends and magnitude are captured well, though NNR fails to capture the coldest temperatures, causing about 1C of misfit on average.(b) Fig. 3. Air temperature at Lake El'gygytgyn.In (a), mean annual air temperature and degree days from NNR (thick lines) and AWS (thin lines), showing that NNR captures the AWS trends well and magnitude best in summer, which is most important to our study.Due to the pronounced warming trend starting in the late-1980s, we used the period 1961-1988 to define the modern era used in this study.In (b), daily NNR is compared to local AWS showing that the trends and magnitude are captured well, though NNR fails to capture the coldest temperatures, causing about 1 • C of misfit on average.

Figure 4 .
Figure 4. NNR daily minimum, maximum, and mean for 1961-1988.This mean is what we used for the 'modern' era in comparisons with paleoconditions.

Fig. 4 .
Fig. 4. NNR daily minimum, maximum, and mean for 1961-1988.This mean is what we used for the "modern" era in comparisons with paleoconditions.

Figure 5 .
Figure5.Modeled ice growth and decay calibration.In both plots, blue represents ice growth (Eq. 1) and red represents ice melt (Eq.2).Fall freeze-up temperatures are from the preceding year but presented like this to highlight summer dynamics.Dotted red lines indicate how much more ice would have been melted, if it were present at end of winter.In A, the black line shows our energybalance modeling results from 1999-2001 and black circles are percentage of ice-covered area on the lake from 2001 remote-sensing, using the ice thickness scale of 100 cm as 100%.Here we have fit the tuning parameters of Eq. 2 to match the melt rate of the previous model.In B, we have tuned parameters of Eq. 1 to match the growth rate measurements from 2009 (blue circles).

Figure 5 .Fig. 5 .
Figure5.Modeled ice growth and decay calibration.In both plots, blue represents ice growth (Eq. 1) and red represents ice melt (Eq.2).Fall freeze-up temperatures are from the preceding year but presented like this to highlight summer dynamics.Dotted red lines indicate how much more ice would have been melted, if it were present at end of winter.In A, the black line shows our energybalance modeling results from 1999-2001 and black circles are percentage of ice-covered area on the lake from 2001 remote-sensing, using the ice thickness scale of 100 cm as 100%.Here we have fit the tuning parameters of Eq. 2 to match the melt rate of the previous model.In B, we have tuned parameters of Eq. 1 to match the growth rate measurements from 2009 (blue circles).(b)Fig.5.Modeled ice growth and decay calibration.In both plots, blue represents ice growth (Eq. 1) and red represents ice melt (Eq.2).Fall freeze-up temperatures are from the preceding year but presented like this to highlight summer dynamics.Dotted red lines indicate how much more ice would have been melted if it were present at end of winter.In (a), the black line shows our energy balance modeling results from 1999-2001 and black circles are percentage of ice-covered area on the lake from 2001 remote-sensing, using the ice thickness scale of 100 cm as 100 %.Here we have fit the tuning parameters of Eq. (2) to match the melt rate of the previous model.In (b), we have tuned parameters of Eq. (1) to match the growth rate measurements from 2009 (blue circles).

Figure 6 .
Figure 6.Model output for the NNR record 1961-2009.In A, maximum ice-melt potential and ice thickness are shown (thick lines) with parameter uncertainties (thin lines) as described in the text.In B, dates of initial freezeup, ice melt, and ice out are shown (thick lines) with parameter uncertainties (thin lines) as described in the text.In C, the model open water season length in shown, indicating that there is usually at least 3 months of ice-free water in the modern environment.

Figure 6 .
Figure 6.Model output for the NNR record 1961-2009.In A, maximum ice-melt potential and ice thickness are shown (thick lines) with parameter uncertainties (thin lines) as described in the text.In B, dates of initial freezeup, ice melt, and ice out are shown (thick lines) with parameter uncertainties (thin lines) as described in the text.In C, the model open water season length in shown, indicating that there is usually at least 3 months of ice-free water in the modern environment.

Fig. 6 .
Fig. 6.Model output for the NNR record 1961-2009.In (a), maximum ice melt potential and ice thickness are shown (thick lines) with parameter uncertainties (thin lines) as described in the text.In (b), dates of initial freeze-up, ice melt, and ice-out are shown (thick lines) with parameter uncertainties (thin lines) as described in the text.In (c), the model open water season length in shown, indicating that there is usually at least 3 months of ice-free water in the modern environment.

Fig. 7 .
Figure 7. Mode the mean NNR uniformly frommaximum ice-m thickness, and thickness are s parameter unce described in the degree cooling, completely mel begins.In terms conditions, likel is required to le to prevent wate exchange.In B freeze-up, ice m shown (thick lin uncertainties (th in the text.In C water season le the value rapidl zero after 3C co modern temper

Table 1 .
Cumulative rainfall of the AWS and NCEP, shown in millimeters and as a percentage of 2006 total.

Table 2 .
Summary of parameter values, uncertainty ranges, and the resultant impact of these ranges on ice on/off dates.

Table 3 .
Models results for modern era 1961Models results for modern era  -2009.   .