Trends and abrupt changes in 104 years of ice cover and water temperature in a dimictic lake in response to air temperature , wind speed , and water clarity drivers

The one-dimensional hydrodynamic ice model, DYRESM-WQ-I, was modified to simulate ice cover and thermal structure of dimictic Lake Mendota, Wisconsin, USA, over a continuous 104-year period (1911–2014). The model results were then used to examine the drivers of changes in ice cover and water temperature, focusing on the responses to shifts in air temperature, wind speed, and water clarity at multiyear timescales. Observations of the drivers include a change in the trend of warming air temperatures from 0.081 C per decade before 1981 to 0.334 C per decade thereafter, as well as a shift in mean wind speed from 4.44 m s before 1994 to 3.74 m s thereafter. Observations show that Lake Mendota has experienced significant changes in ice cover: later ice-on date(9.0 days later per century), earlier ice-off date (12.3 days per century), decreasing ice cover duration (21.3 days per century), while model simulations indicate a change in maximum ice thickness (12.7 cm decrease per century). Model simulations also show changes in the lake thermal regime of earlier stratification onset (12.3 days per century), later fall turnover (14.6 days per century), longer stratification duration (26.8 days per century), and decreasing summer hypolimnetic temperatures (−1.4 C per century). Correlation analysis of lake variables and driving variables revealed ice cover variables, stratification onset, epilimnetic temperature, and hypolimnetic temperature were most closely correlated with air temperature, whereas freeze-over water temperature, hypolimnetic heating, and fall turnover date were more closely correlated with wind speed. Each lake variable (i.e., ice-on and ice-off dates, ice cover duration, maximum ice thickness, freeze-over water temperature, stratification onset, fall turnover date, stratification duration, epilimnion temperature, hypolimnion temperature, and hypolimnetic heating) was averaged for the three periods (1911–1980, 1981–1993, and 1994–2014) delineated by abrupt changes in air temperature and wind speed. Average summer hypolimnetic temperature and fall turnover date exhibit significant differences between the third period and the first two periods. Changes in ice cover (ice-on and iceoff dates, ice cover duration, and maximum ice thickness) exhibit an abrupt change after 1994, which was related in part to the warm El Niño winter of 1997–1998. Under-ice water temperature, freeze-over water temperature, hypolimnetic temperature, fall turnover date, and stratification duration demonstrate a significant difference in the third period (1994–2014), when air temperature was warmest and wind speeds decreased rather abruptly. The trends in ice cover and water temperature demonstrate responses to both long-term and abrupt changes in meteorological conditions that can be complemented with numerical modeling to better understand how these variables will respond in a future climate. Published by Copernicus Publications on behalf of the European Geosciences Union. 1682 M. R. Magee et al.: Trends and abrupt changes in 104 years of ice cover and water temperature

depth and cooler deep waters (Hocking and Straskraba, 1999; Tanentzap et al., 2008). 1 Previous studies have also shown that increased water clarity was correlated with deeper 2 mixed-layer depth (Mazumder and Taylor 1994;Fee et al. 1996;Schindler et al. 1996) and 3 increased hypolimnetic heating rate (Yan 1983). Water clarity has important consequences for 4 photosynthesis and vertical distribution of biota. Accompanied by a warming climate, 5 increased evaporation and changes in precipitation patterns can alter the inputs of nutrients 6 and dissolved organic carbon (DOC) into lakes, resulting in changes in water clarity 7 (Schindler et al. 1990; Schindler et al. 1996). 8 Abrupt shifts, rather than linear changes, can occur in both climate and lake variables over 9 long time scales. In the past, the majority of referenced studies used some form of linear 10 regression analysis to look for changes, which assumes that lakes undergo monotonic changes 11 over time (Van Cleave et al., 2014). This assumption can mask the occurrence of step changes 12 (Liu et al., 2013;North et al., 2013) and can hide the underlying mechanisms that are 13 responsible for the changes (North et al., 2014). As noted in other studies (Mueller et al., 14 2009; Van Cleave, 2014), long-term changes in lake thermodynamic variables can be 15 associated with a pronounced, nonlinear step change that indicates a rapid switch between 16 stable states or regimes (Scheffer et al., 2001;Rodionov, 2004;North et al., 2013). Identifying 17 this type abrupt change in climate drivers and corresponding abrupt or gradual change in lake 18 variables may shed light on the role of changing climate on the corresponding ecosystem. 19 The purpose of this study is to investigate how long-term changes in air temperature, wind 20 speed, and water clarity affect the ice cover and thermal structure of a dimictic Lake Mendota, 21 Wisconsin, USA during the past century using a long-term 104-year simulation model. We 22 hypothesize that changes in lake ice cover (ice-on and ice-off dates, ice cover duration, and 23 maximum ice thickness) and thermal structure (stratification onset, fall turnover date, 24 stratification duration, water temperatures, hypolimnetic heating) variables may be 25 characterized by periods of abrupt change rather than gradual trends, based on observations of 26 rapid change in the climate drivers of air temperature and wind speed. To address this, a one-27 dimensional hydrodynamic model with ice cover is first validated using a long-term (104-28 year) observational dataset. The model is then employed to simulate long-term (1911-2014) 29 ice cover and water temperature in the lake. With the knowledge of lake responses to these 30 past conditions, we aim to reveal how lakes might respond to future changes in air 31 temperature, wind speed, and water clarity. 32 Hydrol. Earth Syst. Sci. Discuss., doi: 10.5194/hess-2015-488, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 19 January 2016 c Author(s) 2016. CC-BY 3.0 License. located in Madison (MSN) Dane County Regional Airport (Truax Field), the same site as that 1 used in 1988. Since Robertson (1989) adjusted all historical data to that collected in 1988, no 2 adjustments are applied to the recent data except for wind. In 1996, a discontinuity in the 3

Light extinction 18
Seasonal Secchi depths can be used to determine the light extinction coefficient for 19 DYRESM-WQ-I. Lathrop et al. (1996) compiled Secchi depth data for Lake Mendota 20 between 1900 and 1993 (1,701 daily Secchi depth readings from 70 calendar years), and 21 summarized the data for six seasonal periods: winter (ice-on to ice-out), spring turnover (ice-22 out to 10 May), early stratification (11 May to 29 June), summer (30 June to 2 September), 23 destratification (3 September to 12 October), and fall turnover (13 October to ice-on). After 24 1993, Secchi depths were obtained from the North Temperate Lakes -Long Term Ecological 25 Research (NTL-LTER) (https://portal.lternet.edu/nis/home.jsp#). For years with no Secchi 26 data, the long-term mean seasonal Secchi depth was used to estimate light extinction. Light

River inflow and outflow 1
Daily inflow measurements have been made on selected Lake Mendota tributaries since 1974, 2 and daily outflow from the lake has been measured since 1975. These measurements were 3 used to calculate total daily inflow and outflow. Total daily inflow and outflow from 1930 to 4 1974 were estimated from measurements made at the gauging station downstream of the lake 5 using the drainage-area ratio method (Maidment 1993). Prior to 1930, total inflow and 6 outflow were calculated using a water budget approach, i.e., the balance of inflow/outflow, 7 precipitation, evaporation, and lake level changes. The inflow/outflow is the residual 8 unknown term of the water balance where other terms include evaporation, rainfall, and water 9 level. The calculation is performed at the interval of water level measurements and the 10 residual term is then distributed evenly across the number of days between water level 11 measurements. Water level of Lake Mendota has been recorded since 1916 at the Yahara 12 River inlet in the north of the lake. Prior to 1916, the long-term mean lake level was assumed 13 in water budget calculations. River water temperature has been recorded at the Yahara River 14 Long-term water temperature records were obtained from Robertson (1989) and the NTL-3 LTER dataset (https://portal.lternet.edu/nis/home.jsp#). The frequency of water temperature 4 data varied widely from only one or two profiles per year to several profiles for a given day, 5 while there were no data collected in some years between 1931 and 1970. The vertical 6 resolution of the water profiles varied from 0.5 m to 2 m and sometimes 5 m when the water 7 column was weakly stratified. (15% reduction) based on the sequential t-test STARS method (Rodionov, 2004). Figure 1c  17 shows Secchi depth for Lake Mendota; however, no statistical trend was obtained due to the 18 incompleteness of the data. 19 Combining the statistically significant breakpoint in air temperature trend that occurred in

Model validation and long-term simulations 27
Here we provide results of model validity comparisons followed by results of long-term 28 changes in both observed (where available) and simulated conditions for: snow and ice 29  Figure 2 shows the comparison of the long-term (from 1911 to 2014) simulation and 4 measured data for blue ice, white ice, and total ice thickness and snow depth. Mean absolute 5 differences between simulated and measured thickness of total ice, blue ice, snow ice, and 6 snow cover are 7.8 cm (n=251; RMSE=7.9 cm), 5.5 cm (n=21; RMSE=3.3 cm), 1.9 cm 7 (n=21; RMSE= 2.3 cm), 4.1 cm (n=49; RMSE= 3.7 cm), respectively. Some discrepancy 8 between the model results and measurements may be due to the 1-D structure of the model. 9 The DYRESM-WQ-I model provides only the lakewide average ice and snow depths, hence 10 no horizontal variation. For a medium-to-large lake, such as Lake Mendota, spatial variations 11 in ice and snow depths should be expected (Bengtsson 1986  Manitoba produced similar errors to Lake Mendota between modeled and observed ice 25 thickness and snow cover. In general DYRESM-WQ-I predicts the ice and snow depths fairly 26 accurately. Figure 3a shows the evolution of simulated ice and snow thickness with multiple 27 ice and snow thickness measurements taken during the winter of 2009-2010. Once ice 28 develops, water beneath the ice continues to freeze, adding to the ice thickness, as heat is 29 conducted from water to the air. In this year, lake ice grew fast initially and then slowed 30 substantially. Eventually, the ice reaches its maximum thickness and stops growing (Ashton 31 1986). Once air temperatures begin to warm up, ice quickly thins until ice out. Overall, the 32 model is able to capture the growing ice cover as well as the temporally variable snow cover 1 found on the lake.  , 1911-1980, 1981-1993, and 1994-2010, (see Table 1 and   bottom sediment although we cannot discount that cold inflows could also penetrate into the 1 lake and affect temperature at the selected depths (20-25 m).After ice break-up and just prior 2 to seasonal stratification, temperature in near-bottom waters is usually identical to surface 3 waters and the whole water column undergoes a period of sustained increase in temperature. 4 Hypolimnetic temperatures then stay relatively constant during the stratified period with 5 limited heat exchange associated with strong temperature gradients across the thermocline.   23 We characterize summer stratification by the dates of the onset of stratification and fall 24 turnover, and the total duration of stratification. Accurately describing these conditions 25 requires high frequency observations, which is challenging for long-term datasets. As a result, 26 we used modeled data to determine these variables. The dates of stratification onset/fall 27 turnover are defined as the day when the surface-to-bottom temperature difference is 28 greater/less than 2 o C (Robertson and Ragotzkie 1990). Figure 11  has become earlier by 12.3 days (p<0.05), and fall turnover has become later by 14.6 days 1 (p<0.05), resulting in the stratification period increasing by 26.8 days (p<0.05). Stratification 2 onset date shows no significant difference among the three regimes (see Table 1); however, 3 fall turnover date is significantly (p<0.05) different between periods 2 and 3 and between 4 period 1 and 3, and the duration of stratification is significantly (p<0.05) different between 5 periods 1 and 3, as shown in Table 1.

Significance of lake drivers 13
Changes in long-term simulated ice cover and thermal variables in Lake Mendota from 14 1911−2014 appear to occur over short time frames within this period, synchronous with rapid 15 changes in measured drivers (i.e., air temperature, wind speed, and water clarity). The 16 variables examined included ice on/off dates, maximum ice thickness, freeze-over water 17 temperature, mid-summer epilimnetic and hypolimnetic temperatures, summer hypolimnetic 18 heating, and dates of stratification onset and fall turnover. To identify the the relationship 19 between the drivers and each of these variables, we employed Pearson correlation analysis on 20 the detrended driver and simulation data. While using this method does not allow us to 21 directly determine causality, we may identify possible related drivers to changes in lake 22 variables and the relative importance of relationship of the three drivers (i.e. air temperature, 23 wind speed, and water clarity) to the lake variables. To calculate correlation coefficients, each 24 lake driver was averaged over a fixed period (e.g. April-May or November-December) and 25 then paired with each of the simulated lake variables from the same period. The averaging 26 period for air temperature and wind speed was chosen based on the fixed period that yielded 27 the best correlation. We chose to employ fixed periods because thorough testing indicated that 28 using a dynamic period introduced seasonal meteorological variations into the analysis. For 29 instance, the correlation coefficient between onset date and the mean air temperature averaged 30 over the prior month is expected to be high simply due to the seasonal variations over the 31 course of the year, i.e., the later in the spring, the warmer air temperature. Seasonal Secchi 1 depths were compiled for six seasonal periods: winter, spring turnover, early stratification, 2 summer, destratification, and fall turnover (Lathrop et al. 1996), and the period closest to the 3 averaging period used for air temperature and wind speed was used for correlation analysis. 4 The averaging periods used for calculating correlation coefficients of each pair of variables 5 are listed in Table 2. We discuss below the correlations for each pair of variables, shown in 6  Fall turnover occurs when the lake water cools, driven largely by colder air temperature and 12 wind-induced mixing. The date of fall turnover is significantly negatively correlated with 13 September wind speed (r = -0.43) so that higher wind speed is related to earlier fall turnover. 14 Air temperature (r = 0.11) and Secchi depth (r = -0.15) are not significantly correlated with 15 the timing of fall mixing. Interestingly, the mid-summer epilimnion-hypolimnion temperature 16 difference has a higher correlation with turnover date (r = 0.62) than any of the three lake 17 drivers. This finding suggests that the influences of spring conditions may be transmitted to 18 the following summer and fall seasons. In other words, the cooler hypolimnion and the greater 19 epilimnion-hypolimnion temperature difference after a warmer and less windy spring may be 20 related to later fall turnover and a longer stratification period. 21 22

Abrubt changes in lake variables 23
To investigate the effects of abrupt changes in air temperature and wind speed on lake ice 24 cover and water temperatures, we used the hydrodynamic model DYRESM-WQ-I to describe 25 changes in several lake variables during 1911−2014. Simulation results were used to examine 26 differences in mean values of these lake variables between specific periods (see Sect. 4.1). 27 For each period of the selected periods, mean lake variables were calculated and the 28 differences between periods were analysed with t-tests to determine if they were significantly 29 different. Table 1  variables shows a shift to warmer air temperature; period 2 to period 3 (1994-2014) represents 1 an abrupt change to lower wind speed; and period 1 to period 3 represents a shift in to warmer 2 air temperature combined with an abrupt change to lower wind speeds. 3

Ice cover
4 Three simulated ice cover variables (maximum ice thickness, ice-on date, and ice-off date) 5 show no significant difference in means between period 1 and 2. In other words, the abrupt 6 change in air temperature trend does not result in a different ice regime even though the ice 7 cover variables are all highly correlated with air temperature (r >0.70). This may be because 8 the change in air temperatures was not of sufficient magnitude to cause a particularly large 9 change in ice cover or it may signify that other drivers are contributing to changes in ice cover 10 variables. Additionally, no significant difference is observed between period 2 and period 3 11 for the ice cover variables since the wind speed and ice cover variables are only weakly 12 correlated. The ice variables do show a statistically significant difference in mean values 13 between periods 1 and 3, indicating that a significant shift in these variables occurs only after 14 a sufficiently large increase in air temperature and an abrupt shift in wind speed within the 15 time between periods 1 and 3. In other words, air temperature need to increase sufficiently to 16 observe a statistically significant difference in ice cover. Ice cover duration, however, shows a 17 significant difference in the mean between all regimes (1−2, 2−3, and 1−3), indicating that 18 distinct differences in ice cover duration can be affected by both trends in air temperature 19 (i.e., there was a large enough change in air temperature between each period) and an abrupt 20 shift in the wind speed. The combined effects of slightly later ice-on dates and earlier ice-off 21 dates during each of the three periods resulted in statistically significant difference in mean 22 ice cover duration values between each of the three periods. 23 Analysis of simulated maximum ice thickness, ice-on, ice-off, and ice cover duration using 24 the method of Rodionov (2004) shows that the most statistically significant timing of the shift 25 in these ice cover variables occurs in the winter of 1997-1998, but a major shift in the air 26 temperature or wind speed data was not observed at that time. indicating that water clarity may act to inhibit heating regardless of changes in wind speed, or 28 it may be acting to filter or mitigate the effects of the wind speed shift. Finally, mean onset 29 date of stratification and mid-summer epilimnetic temperature exhibit no difference among 30 the three periods. This may be due to two processes: (i) the climate signal is being filtered out 31 the lake; or (ii) the external perturbation of the system is not yet strong enough to trigger a 1 major shift in the system's internal dynamics. 2

Ecological significance of long-term changes in lake variables 3
Changes in lake ice cover and thermal structure are of great ecological significance. Summer 4 stratification inhibits the vertical transport of oxygen and nutrients. The vertical temperature 5 gradient in Lake Mendota increased (Fig. 9b) periods, while ice-on, ice-off, and maximum ice thickness only show a significant difference 1 between period one and three, indicating that only with a large change in air temperature and 2 an abrupt shift in wind speeds are change in the ice cover variables statistically different. Mid-3 summer hypolimnetic temperature and fall turnover date both reveal significant (p <0.05) 4 differences in the mean value in 1994, corresponding with the abrupt shift toward lower wind 5 speeds. Some lake variables (under-ice water temperature, freeze-over water temperature, 6 epilimnion-hypolimnion temperature difference, and stratification duration) may not be driven 7 by either the change in air temperature trend or the abrupt shift in wind speed alone, but a 8 shift in the mean of the lake variables does occur in 1994 when both the air temperatures are 9 warmest and the wind speed experienced an abrupt shift. The exact timing of shifts may be 10 difficult to define because of extreme changes in weather in specific years and it may mask 11 the longer term changes in meteorological conditions (i.e. abrupt shifts). 12 Linear trends from the extensive data collection may be misleading as they masks these 13 sudden shifts in the air temperature and wind speed drivers and lake ice cover and thermal 14 structure variables. Examining the mean differences in lake variables in response to shifts in   (1) and (2) (2) and (3) (1) and (3) Lake   1910 1920 1930 1940 1950 1960 1970 1980 1990 1910 1920 1930 1940 1950 1960 1970 1980 1990