Response of water temperatures and stratification to changing climate in three lakes with different morphometry

Water temperatures and stratification are important drivers for ecological and water quality processes within lake systems, and changes in these with increases in air temperature and changes to wind speeds may have significant ecological consequences. To properly manage these systems under changing climate, it is important to understand the effects of increasing air temperatures and wind speed changes in lakes of different depths and surface areas. In this study, we simulate three lakes that vary in depth and surface area to elucidate the effects of increasing air temperatures and decreasing wind speed on lake 10 thermal variables (water temperature, stratification dates, strength of stratification, and surface heat fluxes) over a century (1911-2014). For all three lakes, epilimnetic temperatures increased, hypolimnetic temperatures decreased, the length of the stratified season increased due to earlier stratification onset and later fall overturn, stability increased, and longwave and sensible heat fluxes at the surface increased. Overall, lake depth influences presence of stratification, Schmidt stability, and differences in surface heat flux, while lake surface area influences differences in hypolimnion temperature, hypolimnetic 15 heating, variability of Schmidt stability, and stratification onset and fall overturn dates. Larger surface area lakes have greater wind mixing due to increased surface momentum. Climate perturbations indicate that larger lakes have more variability in temperature and stratification variables than smaller lakes, and this variability increases with larger wind speeds. For all study lakes, Pearson correlations and climate perturbation scenarios indicate that wind speed plays a large role on temperature and stratification variables, sometimes greater than changes in air temperature, and wind can act to either amplify or mitigate the 20 effect of warmer air temperatures on lake thermal structure depending on the direction of local wind speed changes.


Introduction
The past century has experienced global changes in air temperature and wind speed. Land and ocean surface temperature anomalies increased from 1850 to 2012 (IPCC, 2013). Mean temperature anomaly across the continental United States has increased (Hansen et al., 2010), and studies suggest that more intense and longer lasting heat waves will continue in the future 25 (Meehl and Tebaldi, 2004). Additionally, global change in wind speed has been heterogenous. For example, wintertime wind 2 energy increased in Northern Europe (Pryor et al., 2005), but other parts of Europe have experienced decreases in wind speed in part due to increased surface roughness (Vautard et al., 2010), while modest declines in mean wind speeds were observed in the United States (Breslow and Sailor, 2002). Similarly, on regional scales,  showed a decrease in Madison, Wisconsin wind speeds after 1994, but Austin and Colman (2007) found increased wind speeds in Lake Superior, North America. Significant changes to air temperature and wind speed observed in the contemporary and historical periods 5 are likely to continue to change in the future. Lake water temperature is closely related to air temperature and wind speed. Previous studies have mainly focused on warming air temperatures, showing increased epilimnetic water temperatures (Dobiesz and Lester, 2009;Ficker et al., 2017;O'Reilly et al., 2015;Shimoda et al., 2011), increased strength of stratification (Hadley et al., 2014;Rempfer et al., 2010), prolonged 10 stratified period (Ficker et al., 2017;Livingstone, 2003;Woolway et al., 2017a), and altered thermocline depth (Schindler et al., 1990). However, hypolimnetic temperatures have undergone warming, cooling, and no temperature increase (Butcher et al., 2015;Ficker et al., 2017;Shimoda et al., 2011). Wind speed also strongly affects lake mixing (Boehrer and Schultze, 2008), lake heat transfer (Boehrer and Schultze, 2008;Read et al., 2012), and temperature structure (Desai et al., 2009;Schindler et al., 1990). Stefan et al. (1996) found that decreasing wind speeds resulted in increased stratification and 15 increased epilimnetic temperatures in inland lakes, and Woolway et al. (2017a) similarly found that decreasing wind speeds increased days of stratification for a polymictic lake in Europe. In Lake Superior, observations show that the complex nonlinear interactions among air temperature, ice cover, and water temperature result in water temperature increases (Austin and Allen, 2011), contrary to the expected decreases in water temperature from increased wind speeds (Desai et al., 2009). In recent years, our understanding of the effects of air temperature and wind speed on changes in water temperature and stratification has 20 improved (Kerimoglu and Rinke, 2013;Woolway et al., 2017a), but research into the response of lakes to isolated and combined changes in air temperature and wind speed has been limited.
Changes in lake water temperature influence lake ecosystem dynamics (MacKay et al., 2009). For example, increasing water temperatures may change plankton community composition and abundance (Rice et al., 2015), alter fish populations (Lynch 25 et al., 2015), and enhance the dominance of cyanobacteria (Jöhnk et al., 2008). Such changes affect the biodiversity of freshwater ecosystems (Mantyka-Pringle et al., 2014). Furthermore, increased thermal stratification of lakes can intensify lake anoxia (Ficker et al., 2017;Palmer et al., 2014), increase bloom-forming cyanobacteria (Paerl and Paul, 2012), and change internal nutrient loading and lake productivity (Ficker et al., 2017;Verburg and Hecky, 2009). Variations in water temperature impact the distribution, behaviour, community composition, reproduction, and evolutionary adaptations of organisms ( Thomas 30 3 et al., 2004). Improved understanding of response of lake water temperatures and ecosystem response to air temperature and wind speed can better prepare management, adaptation, and mitigation efforts for a range of lakes.
Lake morphometry complicates the response of lake water temperatures to air temperature and wind speed changes because it alters physical processes of wind mixing, water circulation, and heat storage (Adrian et al., 2009). Mean depth, surface area, 5 and volume strongly affect lake stratification (Butcher et al., 2015;Kraemer et al., 2015). Large surface areas increase the effects of vertical wind mixing, an important mechanism for transferring heat to the lake bottom (Rueda and Schladow, 2009), and changes in thermocline depth from warming air temperatures may be dampened in large lakes where thermocline depth is constrained by lake fetch (Boehrer and Schultze, 2008;MacIntyre and Melack, 2010). Lake size has been demonstrated to influence the relative contribution of wind and convective mixing to gas transfer (Read et al., 2012), and lake size can influence 10 the magnitude of diurnal heating and cooling in lakes (Woolway et al., 2016), which both have implications for calculating metabolism and carbon emissions in inland waters (Holgerson et al., 2017). Winslow et al. (2015) showed that differences in wind-driven mixing may explain the inconsistent response of hypolimnetic temperatures between small and large lakes.
Previous research efforts have investigated the response of individual lakes (Voutilainen et al., 2014) and the bulk response of lakes in a geographic region to changing climate (Kirillin, 2010;Magnuson et al., 1990), but few studies have focused on 15 elucidating the effects of morphometry, specifically lake depth and surface area, on changes in lake water temperature in response to long-term changes in air temperature and wind speed.
The purpose of this paper is to investigate the response of water temperatures and stratification in lakes with different morphometry (water depth and surface area) to changing air temperature and wind speed. To do this, we employ an existing 20 one-dimensional hydrodynamic lake-ice model to hindcast water temperatures for three lakes with different morphometry.
These lakes vary in surface area and depth and are nearby (<30 km distance) to experience similar daily climate conditions (air temperature, wind speed, solar radiation, cloud cover, precipitation) over the period 1911-2014. Long-term changes in water temperature, stratification, heat fluxes, and stability were used to investigate how lake depth and surface area alter the response of thermal structure to air temperature and wind speed changes for the three study lakes. 25 4 2 Methods

Study sites
Three morphometrically different lakes, Lake Mendota, Fish Lake, and Lake Wingra, located near Madison, Wisconsin, United States of America (USA), were selected for this study. These lakes are chosen for (i) their morphometry differences, (ii) their proximity to one another, and (iii) the availability of long-term data for model input and calibration. 5 Lake Mendota (43°6' N; 89°24'W; Fig. 1a; Table 1), is a dimictic, eutrophic, drainage lake in an urbanizing agricultural watershed (Carpenter and Lathrop, 2008). The lake stratifies during the summer, and typical stratification periods lasts from May to September. Summer (1 June -31 August) mean surface water temperature is 22.4 °C, and hypolimnetic temperatures vary between 11°C to 15 °C. Normal Secchi depth during the summer is 3.0 meters (Lathrop et al., 1996). Fish Lake (43°17'N; 10 89°39'W; Fig. 1b; Table 1) is a dimictic, eutrophic, shallow seepage lake located in northwestern Dane County. From 1966 to 2001, lake level rose by 2.75 meters due to increased groundwater flow from higher than normal regional groundwater recharge (Krohelski et al., 2002). Krohelski et al. (2002) hypothesized that the increase in recharge may be the result of increased infiltration from snowmelt after increased snowfall and less frost-covered soil. Summer stratification lasts from the beginning of May to mid-September. Mean surface water temperature 23.9°C and hypolimnetic temperatures are normally near 8°C 15 during summer months; however, some years reach temperatures of only 5-6 °C in the hypolimnion due to shortened spring mixing durations. Average Secchi depth during the summer months is 2.4 m. Lake Wingra (43°3' N; 89°26' W; Fig. 1c; Table   1) is a shallow, eutrophic, drainage lake. It stratifies on short timescales of hours to weeks (Kimura et al., 2016), but does not experienced sustained thermal stratification. Summer mean water temperature is 23.9°C, and mean Secchi depth is 0.7 meters.
All three lakes have ice cover during winter months, and a description of ice on the lakes can be found in Magee and Wu 20 (2016)

Model description
To hindcast water temperature and stratification in the three study lakes we use the vertical heat transfer model, DYRESM-WQ (DYnamic REservoir Simulation Model-Water Quality; Hamilton and Schladow, 1997), which employs discrete horizontal Lagrangian layers to simulate vertical water temperature, salinity, and density with input including inflows, 25 outflows, and mixing (Imberger et al., 1978). The model has been previously used on a variety of lake types and is accepted as a standard for hydrodynamic lake modelling (Gal et al., 2003;Hetherington et al., 2015;Imberger and Patterson, 1981;Kara et al., 2012;Tanentzap et al., 2007). DYRESM-WQ adopts a one-dimensional layer structure based on the importance of vertical density stratification over horizontal density variations. A one-dimensional assumption is based on observations that 5 the density stratification found in lakes inhibits vertical motions while horizontal variations in density relax due to horizontal advection and convection (Antenucci and Imerito, 2003;Imerito, 2010). Surface exchanges include heating due to shortwave radiation penetration into the lake and surface fluxes of evaporation, sensible heat, long wave radiation, and wind stress (Imerito, 2010). Surface layer mixing is based on potential energy required for mixing, and introduction of turbulent kinetic energy through convective mixing, wind stirring, and shear mixing (Imerito, 2010;Yeates and Imberger, 2003). Layer mixing 5 occurs when the turbulent kinetic energy, stored in the topmost layers, exceeds a potential energy threshold (Yeates and Imberger, 2003). Yeates and Imberger (2003) improved performance of the surface mixed layer routine within the model by including an effective surface area algorithm (see Eq 32 in Yeates and Imberger, 2003) that reduced surface mixing in smaller, more sheltered lakes. Details of the surface mixed layer algorithm are not reproduced here, but can be found in Eq 27-34 of Yeates and Imberger (2003). Hypolimnetic mixing is parameterized through a vertical eddy diffusion coefficient, which 10 accounts for turbulence created by the damping of basin-scale internal waves on the bottom boundary and lake interior (Yeates and Imberger, 2003). Detailed equations on the simulation of water temperature and mixing can be found in Imberger and Patterson (1981), Imerito (2010), andYeates andImberger (2003).
Sediment heat flux is included as a source/sink term for each model layer. A diffusion relation from Rogers et al. (1995) is 15 used to estimate qsed, heat transfer from the sediments to the water column.

= K
(1) where Ksed represents the sediment conductivity with a value of 1.2 Wm -1 °C -1 , and dT/dz is estimated as: where dT/dz is the temperature gradient across the sediment-water interface, Tw is the water temperature adjacent to the 20 sediment boundary, zsed is the distance beneath the water-sediment interface at which the sediment temperature becomes relatively invariant, and is taken to be 5 m (Birge et al., 1927). Ts derived from Birge et al. (1927) and seasonally variant as follows: where D is the number of days from the start of the year and TD is the total number of days within a year. 25 The ice component of the model, DYRESM-WQ-I, is based on the three-component MLI model of Rogers et al., (1995), with the additions of two-way coupling of the hydrodynamic and ice models and time-dependent sediment heat flux for all horizontal layers. The model assumes that the time scale for heat conduction through the ice is short relative to the time scale of meteorological forcing (Patterson and Hamblin, 1988;Rogers et al., 1995), an assumption which is valid with a Stefan 6 number less than 0.1 (Hill and Kucera, 1983). The three-component ice model simulates blue ice, white ice, and snow thickness (see Eq. 1 and Fig. 5 of Rogers et al., 1995). Further description of the ice model can be found in  and Hamilton et al. (in review). Details on ice cover simulations in response to changing climate for the three lakes can be found in  5 Model inputs include lake hypsography, initial vertical profiles for water temperature and salinity, Secchi depth, meteorological variables, and inflows/outflows. The model calculates the surface heat fluxes using meteorological variables: total daily shortwave radiation, daily cloud cover, air vapor pressure, daily average wind speed, air temperature, and precipitation. Water temperature, water budget, and ice thickness is calculated at 1 hr timesteps. Snow ice compaction, snowfall and rainfall components are updated at a daily time step, corresponding to the frequency of meteorological data input. Cloud 10 cover, air pressure, wind speed, and temperature are assumed constant throughout the day, and precipitation is assumed uniformly distributed. Shortwave radiation distribution throughout the day is computed based on lake latitude and the Julian day. Parameters relevant to the open water period are provided in Table 2

Lake morphometry
Height (m), area (m2), and volume (m3) which describe the hyposgraphic curves for each lake were calculated using bathymetric maps of each lake from the Wisconsin DNR. 20

Initial conditions
Initial conditions for each lake include a temperature and salinity profiles for the first days of the simulations. For Lake Mendota, initial conditions were obtained from the NTL-LTER database on the first day of simulation [NEED CITATION HERE]. For Fish Lake and Lake Wingra, initial conditions after ice off were unavailable for 1911, and were assumed to be the average of all available initial conditions for the lake from ±7 days of the Julian start date for all years with available data. 25

Light extinction coefficient
Seasonal Secchi depths were used to determine the light extinction coefficients. Lathrop et al. (1996) compiled Secchi depth data for Lake Mendota between 1900 and 1993 (1701 daily Secchi depth readings from 70 calendar years), and summarized 7 the data for six seasonal periods: winter (ice-on to ice-out), spring turnover (ice-out to 10 May), early stratification (11 May to 29 June), summer (30 June to 2 September), destratification (3 September to 12 October), and fall turnover (13 October to ice-on). After 1993, Secchi depths are obtained from the North Temperate Lake Long Term Ecological Research (NTL-LTER) program (https://portal.lternet.edu/nis/home.jsp#). Open water and under-ice Secchi depths were collected for various longterm ecological research studies, including the NTL-LTER study, and used here to better characterize temperature profiles 5 throughout the year including under ice cover. Secchi depth data for Fish Lake and Lake Wingra were available only from 1995 to the present and collected from the NTL-LTER program. For years with no Secchi data, the long-term mean seasonal Secchi depths were used. Light extinction coefficients were estimated from Secchi depth using the equation from Williams et al. (1980): where kd is the light extinction coefficient and zs is the Secchi depth (m).

Meteorological data
Meteorological data used to hindcast the historical period consisted of daily solar radiation, air temperature, vapor pressure, wind speed, cloud cover, rainfall, and snowfall over a period of 104 years from 1911 to 2014. Air temperature, wind speed, vapor pressure, and cloud cover were computed as an average of the whole day, while solar radiation, rainfall, and snowfall 15 were the daily totals. Meteorological data was gathered from Robertson (1989), who compiled a continuous daily meteorological dataset for Madison Wisconsin from 1884 to 1988 by adjusting for changes in site location. Appended to this dataset is data from the National Climate Data Center weather station at the Dane County Regional Airport. All data other than solar radiation can be obtained from http://www.ncdc.noaa.gov/, for Madison (MSN), and solar radiation can be obtained from http://www.sws.uiuc.edu/warm/weather/. Adjustments to wind speed were made based on changes in observational techniques 20 occurring in 1996techniques 20 occurring in (McKee et al., 2000 by comparing data from Dane County Airport with that collected from the Atmospheric

Inflow and outflow data 25
Daily inflow and outflow data for Lake Mendota was obtained and described in detail by . Details of data collection and gap-filling can be found there and are not reproduced for brevity. Inflow and outflow data for Fish Lake and Lake Wingra follow a similar process. Inflow and outflow were estimated as the residual unknown term of the water budget 8 balancing precipitation, evaporation, and lake level. USGS water level data from  (http://waterdata.usgs.gov/wi/nwis/dv/?site_no=05406050&agency_cd=USGS&amp;referred_module=sw) was used to estimate inflow and outflow from surface runoff and groundwater inflow. For early years of simulation, where lake level information was not available, the long-term mean lake level was assumed for calculations. Krohelski et al. (2002) determined that surface runoff accounted for two-thirds of inflowing water while groundwater inflow accounted for one-third of total 5 inflow over the period 1990-1991. Using these values, we attributed two-thirds of the inflowing water as surface runoff using air temperatures to estimate the runoff temperature similar to the method for Lake Mendota in  and onethird of inflowing water as groundwater inflow using an average of groundwater temperature measurements (Hennings and Connelly, 2008). For Lake Wingra, water level was recorded sporadically during the period of interest, and was assumed to be the long-term mean lake level for water budget calculations. As in Fish Lake, Lake Wingra has no surface inflow streams, with 10 inflow values attributed equally to direct precipitation, surface runoff, and groundwater inflow (Kniffin, 2011). Groundwater inflow temperatures were estimated using an average of measurements (Hennings and Connelly, 2008), and surface and direct precipitation were estimated as air temperature.

Observation data
Observation data used for model calibration came from a variety of sources. For Lake Mendota, long term water temperature 15 records were collected from Robertson (1989) (1965); and the NTL-LTER program (2012a). Frequency of temperature data varied from one or two profiles per year to several profiles for a given week. Additionally, the vertical resolution of the water profiles varied greatly. For Fish Lake and Lake Wingra, water temperature data were collected from NTL-LTER only from 1996-2014 (2012b). 20

Model calibration and evaluation
Model calibration consisted of two processes: (1) closing the water balance to match simulated and observed water levels and (2) adjusting the minimum water level thickness to match simulated and observed water temperatures for each lake. Water balance for all three lakes was closed using the method described in Section 2.3.5 to match measured water levels to known 25 values and to long term average water levels when elevation information was unknown. Model evaporation rates were not validated; we assume that evaporative water flux and heat flux were properly parameterized by the model. Model parameters were derived from literature values (Table 2). To calibrate water temperature, minimum layer thickness was varied from 0.05 to 0.5 m in intervals of 0.025 m for the period 1995-2000 for all three lakes, similar to the method in Tanentzap et al. (2007) and Weinburger and Vetter (2012). One minimum layer thickness was chosen for all three lakes, and the final thickness was chosen to be 0.125 m as it minimized the overall deviation between simulated and observed temperature values for the three lakes.
Three statistical measures were used to evaluate model output against observational data (Table 3): absolute mean error (AME), 5 root mean square error (RMSE), and Nash-Suttcliffe efficiencies (NS) were used to compare simulated and observed temperature values for volumetrically-averaged epilimnion temperature, volumetrically-averaged hypolimnion temperature, and all individual water temperature measurements for unique depth and sampling time combinations. Simulated and observed values are compared directly, except for aggregation of water temperature measurements to daily intervals where sub-daily intervals are available. Water temperatures were evaluated for the full range of available data on each lake. 10

Analysis
Modelling results were analysed using linear regression, sequential t-test, and Pearson correlation coefficient. Linear regression was used to determine the trend of long-term changes in lake variables. Breakpoints in variables were determined using a piecewise linear regression Ying et al., 2015). A sequential t-test (Rodionov, 2004;Rodionov and 15 Overland, 2005) was used to detect abrupt changes in the mean value of lake variables. The variables were tested on data with trends removed using a threshold significance level of p = 0.05, a Huber weight parameter of h = 2, and a cut-off length L = 10 years. Coherence of lake variables (Magnuson et al., 1990) for each lake and between lake pairs was determined with a Pearson correlation coefficient (Baron and Caine, 2000). The three lakes were paired to compare coherence of lake variables with surface area difference (Mendota/Fish pair), depth differences (Fish/Wingra pair), and both surface area and depth 20 differences (Mendota/Wingra). Additionally, temperature, stratification, and heat flux variables for all three lakes are correlated to air temperature and wind speed drivers, ice date and durations, and to temperature, stratification, and heat flux variables within each lake.
To determine the sensitivity of lake water temperature and stratification in response to air temperature and wind speed, we 25 perturbed these drivers across the range of -10°C to +10°C in 1°C temperature increments and 70% to 130% of the historical value in 5% increments, respectively. For each scenario, meteorological inputs remained the same as for the original simulation and snowfall (rainfall) conversion if the air temperature scenarios increase (decrease) above 0°C. Similarly, the water balance is maintained so that the long-term water levels in both lakes matches the historical record. Inflow temperatures are recalculated for each lake to account for increases or decreases in temperature because of air temperature changes. 30

Model evaluation
Simulated temperatures agreed well with observations for all three lakes (Fig. 3, Table 3). The model was validated with all available data for all three lakes during the period 1911−2014. AME and RMSE for all variables were low and less than standard deviations. NS efficiencies were high (>0.85) and most above 0.90, indicating high model accuracy. 15

Summer water temperatures
Lake Mendota and Lake Wingra had similar increasing epilimnetic water temperature trends, while Fish Lake had a larger increase ( Lake Wingra has an abrupt change after 1930 from 23.13°C to 24.02°C. These increases in epilimnetic temperatures are similar to those found for European lakes (Woolway et al., 2017b) in response to regional climate changes, although Woolway et al. (2017b) demonstrated a substantial increase in annually averaged lake surface water temperatures in the lake 1980s in response to an abrupt shift in the climate, which is not apparent in epilimnion water temperatures for our study lakes. Additionally, Van Cleave et al. (2014) showed a regime shift in July -September Lake Superior surface water temperatures after 1997, driven 25 by El Niño in 1997-1998; however we do not find a similar regime shift in our study lakes, which may be due to geographical differences in meteorology or morphometric differences from the larger Lake Superior.

11
Lake Mendota and Fish Lake hypolimnions were defined as 20−25 m and 13−20 m, respectively, based on the long-term bottom depth of the metalimnion. Lake Mendota has a larger decrease in hypolimnetic temperature than Fish Lake (Table 4), and neither has an abrupt change in temperature nor a significant breakpoint in linear trend during the study period (Fig. 4).
Change in summer (15 July -15 August) hypolimnetic heating was an order of magnitude larger for Mendota than for Fish Lake (Table 4). 5

Stratification and stability
We characterize summer stratification by stratification onset, fall overturn, and duration of stratification (Fig. 5). Onset of stratification and fall turnover were defined as the day when the surface-to-bottom temperature difference was greater than (for stratification) or less than (for overturn) 2°C (Robertson and Ragotzkie, 1990). Lake Wingra experienced only short-term stratification (timescale of days-weeks) and is excluded from this analysis. 10 Lake Mendota has larger trend in earlier stratification onset, fall overturn, and stratification duration than Fish Lake (Table 4), with most of the difference in stratification duration caused by larger change in stratification onset date for Lake Mendota. For both lakes, a significant (p<0.01) shift in onset date occurred at similar times, with shift of 13.3 days earlier for Lake Mendota after 1994 and 15.1 days earlier for Fish Lake after 1993. No change in trend occurred for stratification onset or overturn, but 15 stratification duration shifted from +0.067 days decade -1 to +4.5 days decade -1 after 1940 for Lake Mendota and from -0.19 days decade -1 to +9.6 days decade -1 after 1981 for Fish Lake (Fig. 5).
We quantify resistance to mechanical mixing with a Schmidt number (Idso, 1973). Lake Mendota showed greater stability in general than Fish Lake (Fig. 6) and had a larger trend of change than Fish Lake (Table 4), possibly due to a larger change in 20 stratification and hypolimnion temperature, increasing stability. There was no significant abrupt shift or change in trend for any of the three lakes during the study period.

Surface heat fluxes
Modelled surface heat fluxes included net shortwave, net longwave, sensible heat, latent heat, and total heat fluxes (Fig. 7).
Magnitude of shortwave, longwave, and sensible heat fluxes are similar for all three lakes, but Lake Wingra has a larger 25 magnitude of both latent and net heat fluxes. Net longwave is negative for all three lakes and increased in magnitude (Table   4), and sensible heat flux decreased in magnitude (became less negative; Table 4). There is no significant trend in other surface heat flux variables. Lake Wingra has a much smaller change in trend for longwave radiation than Mendota or Fish, but a larger change in trend for sensible heat flux, indicating that depth likely influences the response of those heat fluxes to air temperature and wind speed changes.

Coherence between lake pairs
Pearson correlations for all variables and lake pairs are significant (Table 5). Shortwave, longwave, sensible, and latent heat fluxes show high correlation for lake pairs, suggesting that morphometry has little impact on variability among lakes. Similarly, 5 epilimnion temperatures have high temporal coherence. However, Fish Lake pairs have lower correlations, which may be a result of changes to lake depth (Krohelski et al. 2002) compared to stable water levels in Mendota and Wingra. Low coherence between the Mendota/Fish pair for hypolimnion temperature and stratification dates suggest that fetch differences impact variability. Stability, however, is lower for pairs with Lake Wingra, indicating that lake depth plays a role in temporal coherence of stability. Similarly, Lake Wingra pairs have lower coherence of net heat flux although the coherence of heat flux 10 components is relatively high. Depth may be influencing a non-linear response of net heat flux that is not present in the components of the flux.

Correlations between lake variables
Generally, direction and magnitude of Pearson correlation between lake variables are similar for each of the three lakes, however, there are some notable exceptions (Fig. 8). Ice off dates are significantly correlated with stratification onset dates 15 and hypolimnetic temperature on Fish Lake, but those correlations do not exist for Lake Mendota. Stratification onset is significantly correlated with hypolimnetic temperature and stability in Lake Mendota, but not significantly correlated on Fish Lake. Summer air temperatures are more highly correlated with stability than summer wind speed for Lake Mendota and Fish Lake, but the opposite is true for Lake Wingra, where summer air temperature is not significantly correlated. Additionally, hypolimnion temperature is more highly correlated with stability in Lake Mendota, whereas epilimnion temperature is more 20 highly correlated with stability in Fish Lake.

Sensitivity to changes in air temperature and wind speed
Response of stratification onset, fall overturn, and hypolimnetic temperature to air temperature and wind speed perturbation scenarios for Lake Mendota and Fish Lake are discussed in the following. Other variables are omitted for brevity and Lake Wingra did not experience prolonged stratification under any sensitivity scenarios, so are excluded from the analysis. In our 25 analysis, we refer to the "base case" as the meteorological values that represent the historical period from 1911-2014, and refer to perturbations in air temperature and wind speed from that base case. the response of median onset dates to changes in air temperature is the same for both increases and decreases from the base case (-2.0 days °C -1 ) for Lake Mendota, but for Fish Lake, the change is magnitude of change is larger for air temperature decreases (-1.5 days °C -1 for temperature increases and +2.7 days °C -1 for temperature decreases). Variability in Lake Mendota onset remains consistent across scenarios, but decreases for Fish Lake as air temperatures increase. This may be from 5 interaction between ice cover and stratification onset on Fish Lake but not on Lake Mendota. For Lake Mendota stratification onset dates have a larger change for wind speed increases than decreases (-3.4 days (m s -1 ) -1 for decreases and +10.5 days (m s -1 ) -1 for wind speed increases), as does Fish Lake(-3.6 days (m s -1 ) -1 for decreases and +8.1 days (m s -1 ) -1 for wind speed increases). Variability in onset dates decreases with lower wind speeds and increases with higher wind speeds.

10
Fall overturn typically occurs slightly early on Lake Mendota than Fish Lake for all scenarios (Fig. 10). For Lake Mendota, stratification overturn dates change at a rate of +0.68 days °C -1 with positive and negative perturbations in temperature, while Fish Lake has a larger change for air temperature increases (+1.81 days °C -1 for temperature increases and -0.77 days °C -1 for temperature decreases) from the historical condition. Standard deviation in overturn dates decreased slightly for Lake Mendota as air temperature increase, but remains consistent for Fish Lake. For stratification overturn dates on Lake Mendota, the change 15 is +13.9 days (m s -1 ) -1 for wind speed decreases and -17.1 days (m s -1 ) -1 for wind speed increases. For Fish Lake, the change is +16.4 days (m s -1 ) -1 for wind speed decreases and -8.5 days (m s -1 ) -1 for wind speed increases. Like onset dates, variability in overturn dates decreases with lower wind speeds and increases with higher wind speeds.
For both lakes, increases in air temperature increase hypolimnetic temperatures, while decreases in wind speed decrease 20 temperatures (Fig. 11). Simulations show that the response of median hypolimnetic temperatures to changes in air temperatures is consisten for Lake Mendota (+0.18°Chypolimnion Cair temperature -1 ) for both increases and decreases from the base case, but not so for Fish Lake (+0.25°Chypolimnion Cair temperature -1 for air temperature increases and -0.18 °Chypolimnion Cair temperature -1 for air temperature decreases). Standard deviations under varying air temperature scenarios remain consistent for both lakes.
Hypolimnion temperatures change inconsistently with increases and decreases in wind speed for both lakes. For Lake Mendota, 25 the change is -1.1°C (m s -1 ) -1 for decreases and +1.8°C (m s -1 ) -1 for wind speed increases. For Fish Lake, the change is -1.2°C (m s -1 ) -1 for decreases and +0.8°C (m s -1 ) -1 for wind speed increases. Variability decreases for lower wind speeds in Lake Mendota, but remains constant for Fish Lake.

Model performance and comparison
The DYRESM-WQ-I model reliably simulated water temperatures over long-term (1911-2014) simulations ( Figure 3, Table   4). Generally, simulated temperatures were lower than observed values. Some may be attributed to timing of observations, which in most instances occur during midday, when water temperatures may be slightly higher than daily averages, as output 5 from the model. Slight deviation is also expected due to averaging of air temperature and wind speeds. In general, thermocline depths were within 1 m of observed values, but some years differ by as much as 2.5 m, contributing additional error in water temperature comparison for depths near the thermocline. for Freshwater Lake (FLake) model (Golosov et al., 2007;Kirillin et al., 2012). Similar to this study, errors in the upper layers were lower than those in the bottom of the water column (Perroud et al., 2009). Fang and Stefan (1996) gave standard errors of water temperature of 1.37°C for the open water season and 1.07°C for the total simulation period for Thrush Lake, MN, 15 similar to those here. Nash-Sutcliff efficiency coefficients for all 3 study lakes were within the ranges found in Yao et al. (2014) for the Simple Lake Model (SIM; Jöhnk et al., 2008), Hostetler, Minlake (Fang and Stefan, 1996), and General Lake Model (GLM; Hipsey et al., 2014) for Harp Lake, Ontario, Canada water temperatures.
Model parameters used to characterize the lake hydrodynamics were taken from literature values. These values may be 20 expected to have small variability between lakes; however, previous studies have shown that many of the hydrodynamic parameters are insensitive to changes of ±10% (Tanentzap et al., 2007). Here the model was validated against an independent dataset for each lake to determine if the model fits measured data and functions adequately, with errors within the range of those from other studies. Adjustments were made to limit uncertainty and errors associated with changes in location and techniques of meteorological measurements. Inflow and outflow measurements were assessed by the USGS for quality 25 assurance and control, but uncertainty for both quantity and water temperature is unknown. The effects of these uncertainties may not be large as the inflow and outflow are small in comparison to lake volume. The combination of uncertainties in parameters and observed data may be high; however, as all parameters and observational methods were kept consistent among the three lakes, the validity of the model in predicting differences among the three lake types is adequate.

15
The main limitation in the model and resulting simulations is the assumption of one-dimensionality in both the model and field data. Quantifying the uncertainty from this limitation can be difficult (Gal et al., 2014;Tebaldi et al., 2005) Small, stratified lakes generally lack large horizontal temperature gradients (Imberger and Patterson, 1981), allowing the assumption of onedimensionality to be appropriate. However, short-term deviations in water temperature and thermocline depth may exist due to internal wave activity, especially in larger lakes (Tanentzap et al., 2007), and spatial variations in wind stress can produce 5 horizontal variations in temperature profiles (Imberger and Parker, 1985). To address the role of internal wave activity and benthic boundary layer mixing, the pseudo two-dimensional deep mixing model by Yeates and Imberger (2003) is employed here. This mixing model has been shown to accurately characterize deep mixing that distributes heat from the epilimnion into the hypolimnion, thus weakening stratification, and the rapid distribution of heat entering the top of the hypolimnion from benthic boundary layer mixing, which strengthens stratification (Yeates and Imberger, 2003). 10 Additionally, light extinction significantly impacts thermal stratification (Hocking and Straškraba, 1999) and light extinction estimated from Secchi depths can have a large degree of measurement uncertainty (Smith and Hoover, 2000), which may result in uncertainty in water temperatures. To address this uncertainty, where available, we use measured Secchi depth values, which has been shown to improve estimates of the euphotic zone over fixed coefficients (Luhtala and Tolvanen, 2013). Secchi depths 15 were unavailable for portions of the simulation period, and average values for the season were used. Analysis comparing using the method of known Secchi depths to both seasonally-varying average Secchi depths and constant Secchi depths for the lakes indicates that seasonally-varying averages do not significantly decrease model reliability when compared to year-specific values, but do show improvement over constant Secchi depths.

Importance of wind speed and other variables 20
While many have addressed the importance of changing air temperatures on water temperatures and water quality (e.g. Adrian et al., 2009;Arhonditsis et al., 2004;O'Reilly et al., 2015;Shimoda et al., 2011), fewer have investigated wind speed as a specific driver of changes to lakes Snortheim et al., 2017). However, results here show that correlations between wind speeds and lake temperature variables are as high as, or higher than, correlations air temperature and lake temperature variables (Fig. 8), highlighting the importance of considering wind speeds as drivers of lake temperature and 25 stratification changes. For many variables (e.g. stratification dates, epilimnetic temperatures, stability), correlation is opposite for air temperature and wind speed variables, indicating that wind speed increases may offset the effect of air temperature increases, while locations with decreasing wind speeds may experience a greater impact on water temperature and stratification than with air temperature increases alone. This is further supported through sensitivity analysis on stratification onset and overturn ( Fig. 9 and 10), which show that for Madison-area lakes, increasing air temperatures and decreasing wind speeds 30 16 have a cumulative effect toward earlier stratification onset and later overturn. This is similar to results of Woolway et al. (2017a), which found that for polymictic Lake Võrtsjärv, wind speed is the key influence on the number of stratified days and that the influence of air temperature increases was minimal; results of Kerimoglu and Rinke (2013), which found that a 30% increase in wind speed can compensate up to a 5.5 K increase in air temperature; and Hadley et al. (2014), which suggest that the combination of increased air temperature and decreased wind are the primary drivers of enhanced stability in Harp Lake 5 since 1979, although no significant change in the timing of onset, breakdown, or duration of stratification was observed.
However, for hypolimnetic temperatures, correlations and sensitivity indicate that decreasing wind speeds may cool hypolimnetic temperatures, while increasing air temperatures warm hypolimnetic temperatures. Arvola (2009) showed that hypolimnion temperatures were primarily determined by the conditions that pertained during the previous spring turnover, which is consistent with our results showing significant (p<0.01) correlation between hypolimnion temperatures and wind 10 speed (Fig. 8), but no significant correlation with air temperature or summer conditions. This could explain the conflicting results of previous research showing both warming and cooling trends in different lakes (Gerten and Adrian, 2001). Hindcasted hypolimnion temperatures (Fig. 4) show decreasing trends for Lake Mendota and Fish Lake. Combining the effects of air temperature and wind speed, it suggests that wind speed decreases are a larger driver to hypolimnetic water temperature changes than increasing air temperatures for both lakes. This is particularly notable as current research into changes in lake 15 water temperature and stratification have been dominated by studying air temperature effects to the neglect of the role of wind speed changes.
Ultimately, lake warming or cooling may depend on the magnitudes and directions of changes of air temperature, wind speed, and other variables as climatic variables humidity, cloud cover, and solar radiation and water clarity variables are important in 20 determining lake water temperatures. Schmid and Köster (2016) demonstrated that 40% of surface water warming in Lake Zurich was caused by increased solar radiation, and Wilhelm et al. (2006) showed that daily extrema of surface equilibrium temperature responded to shifts in wind speed, relative humidity, and cloud cover in addition to changes in air temperature.
However, neither study looked at lakes with seasonal ice cover, which may not account for changes in ice sheet formation and the resulting influence on lake water temperatures (Austin and Colman, 2007). Other studies have demonstrated that ice cover 25 changes do not directly influence summer surface water temperatures (Zhong et al., 2016), in agreement with our modelling results (Fig. 8). Changes in underwater light conditions from increased dissolved organic carbon concentrations combined with reduction in surface wind speeds can result in cooling whole-lake average temperatures despite substantial air temperature increases, as was the case for Clearwater Lake, Canada (Tanentzap et al., 2008). Water clarity has seen both increases and decreases since the early 1990's (Rose et al., 2017), with precipitation playing a critical role in year-to-year variability (Rose 30 et al., 2017). Further investigation into the combined effects of these climatic and lake-specific variables is warranted.

Lake depth
Lake depth plays a key role in determining thermal structure and stratification of the three lakes in this study. Even under the extreme increases in air temperature, Lake Wingra remained polymictic and did not become dimictic like Lake Mendota or Fish Lake. Additionally, Schmidt stability exhibited no trend on the shallow lake, unlike for the deeper two (Table 4). Due to 5 lower heat capacity, shallow lakes respond more directly to short-term variations in the weather (Arvola et al., 2009), and heat can be transferred throughout the water column by wind mixing (Nõges et al., 2011). This was particularly evident in correlations among drivers and lake variables, where air temperature did not have a significant correlation with stability for Lake Wingra, but wind speed was highly correlated. For shallow lakes, wind speed may be a larger driver to temperature structure and stability, with the importance of air temperature increasing with lake depth. Deep lakes have a higher heat 10 capacity so that greater wind speeds are required to completely mix the lake during the summer months, resulting in more temperature stability and higher Schmidt stability values for deeper Lake Mendota and Fish Lake. Our study is consistent with previous research showing mean lake depths can explain the most variation in stratification trends and lakes with greater mean depths have larger changes in their stability (Kraemer et al., 2015). Overall, Lake Wingra had a larger magnitude of latent and net heat fluxes than the deeper lakes. Diurnal variability in surface temperatures is larger for shallow lakes, promoting increased 15 latent heat fluxes in these lakes (Woo, 2007). This increased response may also explain the larger change in trend for sensible heat flux since Lake Wingra responds more quickly to changes in air temperature, thus, have a larger change in sensible heat flux during each day. Interestingly, net heat flux of Lake Wingra is less coherent with the deeper lakes than the deep lakes are with each other. This may be due to the combination of more extreme temperature variability, increasing sensible and latent heat fluxes during the open water season and the lower sensitivity of ice cover duration in Lake Wingra compared to the deeper 20 lakes . Ice cover significantly reduces heat fluxes at the surface (Jakkila et al., 2009;Leppäranta et al., 2016;Woo, 2007), and larger changes in ice cover duration for Lake Mendota and Fish Lake compared to Lake Wingra would reduce synchrony of heat fluxes among the three lakes.

Surface area
Lake surface area impacts the effects of climate changes on water temperatures and stratification. Air temperature is 25 significantly correlated (p<0.01) with epilimnion temperature for all three lakes, as is wind speed (p<0.05). Increasing air temperatures are well documented to increase epilimnetic water temperatures (Livingstone, 2003;Robertson and Ragotzkie, 1990), since air temperature drives heat transfer between the atmosphere and lake (Boehrer and Schultze, 2008;Palmer et al., 2014). However, wind mixing can act as a mechanism of heat transfer (Nõges et al., 2011), and cool the epilimnion through 18 increased surface mixed-layer deepening. Decreasing wind speeds may increase epilimnion temperatures above that from air temperature increases alone (Fig. 8). Surface area plays a role in lake-wide average vertical heat fluxes from boundary processes (Wüest and Lorke, 2003), and the model accounts for this by including an effective surface area algorithm to scale transfer of momentum from surface stress based on lake surface area (Yeates and Imberger, 2003). This increases transfer momentum from surface stress and vertical heat transfer for lakes with larger fetch. Accounting for this larger fetch increases 5 mixing and vertical transfer of heat to bottom waters, reducing epilimnion water temperatures (Boehrer and Schultze, 2008) and increasing the rate of lake cooling (Nõges et al., 2011). For this reason, Lake Mendota with the large fetch experiences a smaller increase in epilimnetic water temperature compared to Fish Lake (Table 5). Additionally, momentum from surface stress scales linearly with lake area and non-linearly with wind speed (Yeates and Imberger, 2003, see Eq. 31 and 33), making momentum from surface stress, and thus, mixing, stratification, and hypolimnion temperatures more variable for lakes with 10 larger fetch and even more variable when wind speed is increased (see Fig. 9-11). Greater variability in momentum and mixing corresponds to larger variability of Schmidt stability for Lake Mendota, with the larger surface area. Greater transfer of momentum in Lake Mendota results in the slightly deeper thermocline for the larger surface area lake (~10 m in Lake Mendota and ~6 m in Fish Lake), which may play a role in filtering the climate signals into hypolimnion temperatures. Low hypolimnetic temperature coherence between Mendota and Fish suggest that lake morphometry plays a role. This result is consistent with 15 other studies that show lake morphometry parameter affects the way temperature is stored in the lake system (Thompson et al., 2005). Increased momentum on Lake Mendota from the larger surface area may also limit the impact of ice off dates on stratification onset and hypolimnetic temperatures because the lake has ample momentum to sustain mixing events regardless of ice off dates, while Fish Lake's small surface area limits mixing making ice-off dates and stratification more highly correlated. 20

Conclusion
The combination of increasing air temperatures and decreasing wind speeds in Madison-area lakes resulted in warmer epilimnion temperatures, cooler hypolimnion temperatures, longer stratification, increased stability, and greater longwave and sensible heat fluxes. Increased stratification durations and stability may have lasting impacts on fish populations (Gunn, 2002;Jiang et al., 2012;Sharma et al., 2011) and warmer epilimnion temperatures affects the phytoplankton community (Francis et 25 al., 2014;Rice et al., 2015). Shallow lakes respond more directly to changes in climate, which drives differences in surface heat flux compared to deeper lakes, and wind speed may be a larger driver to temperature structure than air temperatures, with importance of air temperatures increasing as lake depth increases. Larger surface area lakes have greater wind mixing, which influences differences in temperatures, stratification, and stability. Climate perturbations indicate that larger lakes have more variability in temperature and stratification variables than smaller lakes, and this variability increases with greater wind speeds.
Most significantly, for all three lakes, wind speed plays as large as, or a larger role in temperature and stratification variables than does air temperatures. This reveals that air temperature increases are not the only climate variable that managers should plan for when planning mitigation and adaptation techniques. Previous research has shown uncertainty in the changes in hypolimnion water temperatures for dimictic lakes, however the perturbation scenarios indicate that while increasing air 5 temperature always increases hypolimnion temperature, wind speed is a larger driving force, and the ultimate hypolimnion temperature response may be primarily determined by whether the lake experiences an increase or decrease in wind speeds.

5
R is the range of onset, overturn, or duration, M is the mean date for onset, overturn, or duration length, and SD is the standard deviation in dates for the study period.