Effects of Climate Change on Lake Thermal Structure and Biotic Response in Northern Wilderness Lakes

Mark B. Edlund 1,* ID , James E. Almendinger 1, Xing Fang 2 ID , Joy M. Ramstack Hobbs 1, David D. VanderMeulen 3, Rebecca L. Key 3 and Daniel R. Engstrom 1 ID 1 St. Croix Watershed Research Station, Science Museum of Minnesota, Marine on St. Croix, MN 55047, USA; dinger@smm.org (J.E.A.); jramstack@smm.org (J.M.R.H.); dengstrom@smm.org (D.R.E.) 2 Department of Civil Engineering, Auburn University, Auburn, AL 36849, USA; xzf0001@auburn.edu 3 Great Lakes Inventory and Monitoring Network, National Park Service, 2800 Lake Shore Drive East, Ashland, WI 54806, USA; david_vandermeulen@nps.gov (D.D.V.); rebecca_key@nps.gov (R.L.K.) * Correspondence: medlund@smm.org; Tel.: +1-651-433-5953


Introduction
Significant ecological changes are occurring in remote northern lakes [1][2][3].These changes do not fit the traditional paradigm of nutrient loading and eutrophication due to obvious land use and may be the result of recent climate warming.Climate change has the potential to severely disrupt aquatic ecosystems directly though changes in temperature, wind, and precipitation, indirectly through watershed effects, and in concert with man-made stressors such as nutrient pollution, invasive species, and land-use change.We currently have only a limited scientific grasp of how these forces will interact or how lakes will respond, yet predicting these effects will be critical to both resource management and public understanding of changes that have already begun.
Recent studies have documented a series of possible climate-induced changes in boreal-region lakes, including a longer ice-free season [4][5][6][7], stronger or weaker thermal stratification [8,9], increased inputs of dissolved organic carbon and increased carbon burial [2,10], shifts in algal communities [1], and an increased frequency of cyanobacterial blooms [3].These changes have been noted in remote lakes far removed from direct human disturbance, with the strongest evidence coming from the

Study Sites
Eight lakes located in two National Park units were identified for thermal modeling and paleoecological reconstruction (Figures 1-3).Four lakes were in VOYA: Cruiser, Ek, Little Trout, and Peary lakes (Figure 2).In ISRO, study sites were Ahmik, Harvey, Richie, and Siskiwit lakes (Figure 3).These lakes were chosen to capture a gradient of lake types across the boreal ecoregion ranging from small and shallow to large and deep, and they also representing a spectrum of lake productivity (Table 1).Notes: 1 Latitude (Lat) and longitude (Lon) in decimal degrees for lake center estimated from Google Earth.Altitude for Voyageurs National Park (VOYA) lakes from National Park Service (NPS) bathymetry data; altitude for Isle Royale National Park (ISRO) lakes estimated from Google Earth.Lake area, watershed area (Aws), inlets, and outlets from Ulf Gafvert, NPS-GLKN, personal communication, 2012.Maximum lake depth (Zmax), Secchi-disk depth (Zseci), TP (total phosphorus), Chl a (chlorophyll-a), and dissolved organic carbon (DOC) from Joan Elias, NPS, personal communication, 2012, based partly on Great Lakes Network (GLKN) lake monitoring.VOYA lakes were sampled 2006-2011, and ISRO lakes 2007-2011.TP and Chl a were sampled three times each year (ice-free season); DOC was sampled three times in the first year and once per year (mid-summer) thereafter.Chlorophyll-a values are based on fluorometric readings.Earlier values from spectrophotometer analyses were adjusted to be consistent with fluorometric readings.For each variable, values in the table above are the multi-year average of each year's average value.

Model Mechanics
The MINLAKE models solve the one-dimensional, unsteady differential equation describing heat transfer with depth in a lake, i.e., the change in temperature over time for each depth in a lake [24].The equation is solved numerically with a finite-difference scheme, thereby describing change in temperature over discrete time intervals (e.g., one day) for discrete layers in a lake (e.g., one meter increments).Heat transfer related to atmospheric exchanges (radiation, evaporation, conduction) is applied to the top-most layer, whereas short-wave solar radiation is applied as a heat source to all layers down to the effective depth of light penetration.MINLAKE2012 has submodels handling the details of ice formation [28] and heat transfer across the sediment/water interface [22,29].The upper meter of the lake is finely divided into at least eight layers (depths of 2. 5,5,10,20,40,60,80, and 100 cm) to allow the model to simulate the dynamics of ice formation and melting [28].Including the sediment heat capacity adds an important storage component that can smooth seasonal water temperatures, especially for shallow lakes [24].The MINLAKE2012 model also solves the one-dimensional, unsteady transport equation for changes in dissolved oxygen (DO) concentration with depth [24].We did not use the DO capabilities of MINLAKE2012 in this project; however, we recognize that changes in DO may be a primary driver of ecological change in lakes with altered thermal structure due to climate warming.

Model Input Variables
MINLAKE requires geometric, optical, and thermal properties for each lake, as well as a time series of weather data to drive the model [24].Lake and watershed areas, average values for Secchi disk depth, total phosphorus, chlorophyll-a, and dissolved organic carbon were measured from 2006 to 2011 (Table 1).In accordance with Stefan et al. [30] and Heiskary and Wilson [31], lakes were categorized by area and depth (Table 2) and range from small and shallow (Ahmik in ISRO) to large and deep (Siskiwit in ISRO), with the other six lakes being medium in area, depth, or both (Tables 2 and 3).Area and depth influence light interception, wind fetch, and lake volume, and are thus critical in determining the likelihood of thermal stratification.Based on geometric ratio, a measure of the likelihood of stratification [32], the study lakes range from shallow and small lakes with a weak likelihood of stratification to deep and medium/large lakes with a strong likelihood (Table 3).Trophic status was inferred from the Secchi depth, total phosphorus (TP), and chlorophyll-a measurements [31,33].About half the study lakes were oligotrophic, and about half mesotrophic, depending on the parameter (Table 3).Richie Lake (VOYA) reached the eutrophic category, based on Secchi depth.Note: 1 Modified from Table 2 in Reference [30] and on Reference [31].

Park
The optical, physical, chemical, and biological characteristics of each lake required by the model were inferred from the data discussed above [24].Based on the total light attenuation determined by Secchi depth and light attenuation due to chlorophyll content, light attenuation due to the water matrix was determined by taking the difference of the two (Table 4).Water-matrix attenuation is generally proportional to dissolved organic carbon (DOC), as expected (compare with Table 1).The wind-sheltering coefficient was estimated from the lake area.Multipliers of diffusion coefficients for the metalimnion and hypolimnion were inferred from the lake size based on empirical relations [24].
Weather data include daily values for air temperature, dew point temperature, wind speed, solar radiation, and precipitation as rain or snow [24].Most of these data are available from National Weather Service (NWS) first class weather stations, of which there are six in Minnesota, including Duluth and International Falls, for which we compiled climate data for the period 1961-2011.Weather data for International Falls were used for modeling the VOYA lakes, and data for Duluth were used for modeling the ISRO lakes.While lakes within VOYA are relatively close (about 30-50 km) to International Falls, lakes in ISRO are about 240-320 km from Duluth.All values based on Reference [24] and references therein for total light attenuation, chlorophyll-a (chl a) light attenuation, and wind sheltering coefficient (Wstr).Sb20, sediment oxygen demand (SOD) at 20

Model Runs and Output Variables
MINLAKE was run for a 51-year period, from mid-April 1961 through 2011.The model was run from mid-April, when lakes are approximately isothermal, to provide uniform initial conditions for all lakes; to reduce concerns about thermal inertia, temperature output from the first partial year (1961) was discarded to allow for model equilibration, and all further post-processing was based on the 50 full years of output from 1962 to 2011.Even though MINLAKE provides useful data summaries, significant post-processing was needed to build the annual output variables selected to be compared with the paleoecological proxies.Scripts were written within the R statistical framework [34] to process the raw MINLAKE output and produce the summary plots, tables, and statistical analyses of the selected annual output variables.
We selected output variables that seemed representative of factors that could influence biological activity in lakes.To allow comparison with the paleoecological records, these output variables were summarized as time series of annual values.Water temperatures during the growing season are obviously important in affecting not only biological activity in the upper waters but also microbiological activity at depth.Hence, average summertime (June-July-August: JJA) shallow and bottom water temperatures were calculated for each year of simulation.The shallow water temperature was taken at the 1 m depth in the model, and the bottom water temperature was taken at a depth corresponding to 0.9*Zmax, but no closer than 1 m to the bottom to avoid any gradients near the sediment surface.Temperatures at these depths should be representative of epilimnetic and hypolimnetic temperatures for those lakes deep enough to stratify.The time-series plots of these annual mean temperatures suggested that their year-to-year variability was increasing over time, and we quantified this variability by calculating a nine-year running standard deviation.We chose a nine-year window to give an approximately decadal measure of variability that could be centered on the middle year of the window.Windows of seven to 11 years resulted in very similar patterns of variability over time.
The formation of thermal gradients can affect lake mixing and thus the distribution of key nutrients and escalation of redox-sensitive microbiological activity.The depth of the steepest gradient and whether it forms earlier, dissipates later, and has greater continuity through the year are potential mechanisms by which climate warming may impact lake biology.We calculated thermal gradients with the central-difference method (i.e., negative T 2 − T 1 divided by Z 2 − Z 1 , applied at the midpoint (average) depth between Z 1 and Z 2 ).The thermocline depth was defined as the depth of the steepest gradient each day [35].Output was tested for a variety of gradients, from 1 to 4 • C per meter.While a gradient of 1 • C per meter is commonly invoked as defining a thermocline [35], the model simulations indicated that such gradients formed and dissipated frequently in the model during the spring and fall and thus did not seem representative of features that could strongly isolate top and bottom waters.Hence, we focused on steeper gradients in order to identify less-easily formed but more persistent, stable features with a greater potential to isolate top from bottom waters.The steeper the selected gradient, the fewer days each year attained that gradient.Gradients of 3 to 4 • C per meter or larger were too infrequent to be considered, and so we settled on gradients of 2 • C per meter.For each year of the model run, we determined the date the gradient first formed, the total number of days and the longest continuous period the gradient existed, and the date the gradient dissipated.In addition, for those days that achieved the selected gradient or larger, we also calculated the average summertime depth of the steepest daily gradient (i.e., the thermocline depth).We treated all model output equally without regard to lake depth, recording all gradients as they occurred, even in shallow lakes.However, we note that such gradients have short duration and probably limited impact in shallow lakes that are generally well-mixed.
The ecological impact of thermocline formation can be profoundly magnified when hypolimnetic anoxia releases phosphorus from the sediments and spurs noxious algal blooms [36].Whether a hypolimnion goes anoxic depends on several factors.The dissolved oxygen (DO) depletion in the hypolimnion is driven by the respiration of benthic bacteria and other organisms, and hence is directly proportional to the lake-bottom area in contact with the hypolimnion; we here used the thermocline area as a proxy for this hypolimnion/sediment-contact area.Hypolimnetic anoxia should also be directly proportional to the duration of time that a thermocline has isolated the bottom waters.In contrast, anoxia should be inversely related to hypolimnetic volume simply because of the larger amount of DO to be consumed.For simplicity, we calculated hypolimnetic volume as the lake volume below the thermocline (the plane of the steepest temperature gradient), rather than defining a metalimnion of finite thickness.Combining these factors, we created a simple anoxia susceptibility index of thermocline duration (days), times thermocline area (m 2 ), and divided this by hypolimnetic volume (m 3 ).Dividing volume by area produces a mean hypolimnetic thickness, and hence the final dimensions of the index are days of thermocline duration per meter of hypolimnion thickness.Whether hypolimnetic anoxia results in the release of substantial phosphorus from the sediments depends on the form and concentration of the sediment-bound phosphorus [37], and whether the algal community responds to such a phosphorus release depends on its timing [38].
Ice cover is perhaps the most obvious thermally driven feature of a lake in cold regions that could affect lake biology.Although photosynthesis can take place below ice if snow cover is not too deep, algal biomass and productivity are generally small under the ice during winter [35], indicating that the growing season for aquatic vegetation is largely determined by the length of the ice-free season.All other factors remaining equal, the length of the growing season can directly alter the annual net primary productivity and rates of nutrient cycling.Length of ice cover further affects the eventual redox state of the lake water over the course of winter, which in turn can affect the fish community and nutrient release from the sediments.Consequently, MINLAKE output was processed to extract the dates of first open water (ice-out) in the spring, last ice in the spring, first ice in winter, and last open water in winter to search for trends in how early ice-out occurred in the spring or how late ice-on occurred in the winter.In addition, the durations of continuous ice cover in the winter and open water in the summer were calculated for each year of the simulations.

Model Statistics
Model goodness-of-fit compared modeled water column thermal profiles with observed values.The observed data comprised temperature profiles collected several times per year on each lake, from 2006 to 2010 for VOYA lakes and from 2007 to 2010 for ISRO lakes.Multiple field measurements taken within a 0.1-m depth on the same day were considered duplicates and replaced with their averaged depth and temperature.The modeled data set comprised simulated temperatures for the same day and depth as observed data.Modeled temperatures at the observed depth were calculated by linear interpolation of temperatures at the depth midpoints of the adjacent model layers.Model goodness-of-fit was qualitatively assessed by visual comparison of observed temperature profiles with simulated profiles (Supplementary Materials).In addition to visual inspection, quantitative assessment of model goodness-of-fit used the coefficient of determination, or R-squared (R 2 ; range 0 to 1), to summarize the deviations between each observed value and its corresponding simulated value.The Nash-Sutcliffe coefficient of efficiency (NSE; range negative infinity to 1) was also calculated to estimate model goodness-of-fit.It is similar to R 2 in that it attempts to quantify the fraction of variance explained by the model.If a model is perfect, the predicted model value (x-axis) and the observed value (y-axis) should fall exactly on the 1:1 line.The difference is that in linear regression, the regression line is optimally placed to minimize residuals and thus maximize R 2 , whereas in the NSE calculation there is no such optimization.R 2 is the upper limit to NSE, and the closer the regression intercept and slope are to 0 and 1, respectively, then the closer NSE will be to R 2 .
Change over time for weather input and model output variables was assessed with the simple but robust method of dividing the 50 years of data into two equal time periods (1962-1986 and 1987-2011) and determining the significance of their differences with a Mann-Whitney U Test [34].The significance value does not identify the direction of change, only that the two halves of the data output are different.Direction of change may be inferred from the loess curves that were plotted over time for selected output variables with default lowess settings and a span of 0.75 [34].

Core Collection
Sediment cores were recovered from central depositional basins from the eight lakes as part of other projects [8,13,39,40].All cores were retrieved with a 6.5-cm inner diameter polycarbonate tube fitted with a piston and operated from the lake surface from an anchored boat.Rigid drive rods were used to push the tube into the lake sediments while the piston was held in place by a cable [41].Sediment cores were extruded vertically from the coring tube in the field at 1-cm increments to 60-85 cm core depth and placed in polypropylene jars; the remainder of the core was capped.The material was transported to 4 • C laboratory storage for further processing.Sediment samples were homogenized and subsampled for loss-on-ignition and diatom analyses.Remaining sediments were freeze-dried within the sample jars for 210 Pb dating and geochemical analysis and are archived at the St. Croix Watershed Research Station (Science Museum of Minnesota).

Isotopic Dating and Geochemistry
Sediment cores were analyzed for 210 Pb activity to determine age and sediment accumulation rates for the past 150 to 200 years.Lead-210 activity was measured from its daughter product, 210 Po, which is considered to be in secular equilibrium with the parent isotope.Aliquots of freeze-dried sediment were spiked with a known quantity of 209 Po as an internal yield tracer and the isotopes were distilled at 550 • C after treatment with concentrated HCl.Polonium isotopes were then directly plated onto silver planchets from a 0.5 N HCl solution.Activity was measured for 1-3 × 10 5 s using an Ortec alpha spectrometry system.Supported 210 Pb was estimated by mean activity in the lowest core samples and subtracted from upcore activity to calculate unsupported 210 Pb.Core dates and sedimentation rates were calculated using the constant rate of supply model [42].Dating and sedimentation errors represented first-order propagation of counting uncertainty [43].
Dry-density (dry mass per volume of fresh sediment), water content, organic content, and carbonate content of sediments were determined by standard loss-on-ignition techniques [44].Organic matter content and accumulation were converted to whole basin dry mass organic carbon (OC) burial using the sediment focusing factor for each lake and an organic matter-OC conversion factor [45,46].

Diatom Analysis
Diatoms are a group of microscopic algae that are common in all aquatic habitats and characterized by their ornamented cell wall made of biogenic silica.Diatom cell walls accumulate in lake sediments to provide a record of historical diatom productivity and community structure through time.Sediment was prepared for diatom analysis following Ramstack et al. [39].Coverglasses with diatoms were permanently attached to microscope slides using Zrax mounting medium (MicrAP, Philadelphia, PA, USA), and all slides and material are archived at the St. Croix Watershed Research Station, Marine on Saint Croix, MN, USA.

Statistical Approach
To summarize the timing and degree of change in diatom communities within cores from ca. 1850 to 2010, detrended correspondence analysis (DCA) was used.DCA is an indirect gradient analysis technique used to identify major axes of variation in species data.Axis 1 DCA scores are scaled in standard deviation (SD) units and represent units of species turnover based on beta-diversity [59].Diatom species occurring at 1% or greater abundance in more than one core section were included in the analysis, with downweighting of rare taxa, detrending by segments, and nonlinear rescaling applied using the software R [34,39].The DCA axis 1 scores were plotted stratigraphically to visualize the amount of turnover in SD units; the standardized units of DCA allow for direct comparison among cores.Axis scores were plotted for each lake together with the downcore relative abundance of predominant diatom species (two or more occurrences at >5% relative abundance).The presentation of downcore data is truncated at 1960-2010 for consistency with the thermal model output.

Thermal Modeling
Results of all model runs are presented in Supplemental Materials (Figures S1-S8), arranged alphabetically by lake name.For each lake, figures show 50-year time series of mean summer (JJA) epilimnion temperature, hypolimnion temperature, epilimnion temperature variability (nine-year running standard deviations), hypolimnion temperature variability, calendar days with a 2 • C/m thermocline, average summertime depth of thermocline, days of continuous thermocline duration, and days of continuous thermocline duration per m of hypolimnion thickness.

Climate Drivers
Annual mean air temperatures increased significantly at both the Duluth and International Falls weather station sites over the 50 years of model runs (Figures 4 and 5; Table 5).However, much of this warming was due to an increase in winter temperatures; a summer temperature increase was significant (p < 0.05) only for the Duluth station and not for International Falls.Although an increase in winter temperatures may seem removed from summer water temperatures and algal productivity, warmer winters would contribute to thinner ice and thus allow more of the spring energy inputs to warm the water rather than being consumed by the latent heat needed to melt ice.Annual solar radiation was significantly different between the two time periods (Table 5), but appeared to decrease rather than increase (Figures 4 and 5), and any change was non-significant during summer for both stations (Table 5).Hence, changes in radiation were unlikely to be the driver of warmer water temperatures.Wind speed changed significantly for all seasons and the annual average (Table 5); in this case, the change was a decrease (Figures 4 and 5), which has implications for thermocline formation, stability, and depth.Snowfall showed some changes during fall and winter (Table 5), but rainfall did not change significantly in any season over the 50-year period.
Water 2017, x, x 12 of 34 rather than increase (Figures 4 and 5), and any change was non-significant during summer for both stations (Table 5).Hence, changes in radiation were unlikely to be the driver of warmer water temperatures.Wind speed changed significantly for all seasons and the annual average (Table 5); in this case, the change was a decrease (Figures 4 and 5), which has implications for thermocline formation, stability, and depth.Snowfall showed some changes during fall and winter (Table 5), but rainfall did not change significantly in any season over the 50-year period.Annual means

Model Goodness-of-fit
The model fit for each lake may be assessed qualitatively by scanning the depth profiles of simulated versus observed temperatures (Supplementary Materials; Figures S9-S18).In general, the model fits were quite good for the deeper lakes (Cruiser and Little Trout in VOYA and Richie and Siskiwit in ISRO), with a wide range of temperatures and thermoclines situated at about the correct depth.While the simulated temperatures at any depth might be off by several degrees, the overall shapes of the profiles were similar between the simulated and observed values.The model had more trouble with the shallower lakes, particularly for the ISRO lakes, Ahmik and Harvey.The profiles for the shallower lakes in VOYA, Ek and Peary, appeared reasonable but not quite as good as those for the deeper lakes.
The quantitative measures of model goodness-of-fit support these qualitative observations (Table 6).For lakes within both national parks, the NSE increased as maximum depth (Zmax) and range of observed temperatures both increased.NSE values were generally large and just a few points below their effective maximum (R 2 ) in most cases.Commonly, an NSE value of 0.5 or greater is considered to be acceptable [60].The only lake falling below this threshold was Ahmik Lake in ISRO, whose negative NSE indicated that a simple mean was a better predictor of temperature than was the International Falls, MN Annual means

Model Goodness-of-Fit
The model fit for each lake may be assessed qualitatively by scanning the depth profiles of simulated versus observed temperatures (Supplementary Materials; Figures S9-S18).In general, the model fits were quite good for the deeper lakes (Cruiser and Little Trout in VOYA and Richie and Siskiwit in ISRO), with a wide range of temperatures and thermoclines situated at about the correct depth.While the simulated temperatures at any depth might be off by several degrees, the overall shapes of the profiles were similar between the simulated and observed values.The model had more trouble with the shallower lakes, particularly for the ISRO lakes, Ahmik and Harvey.The profiles for the shallower lakes in VOYA, Ek and Peary, appeared reasonable but not quite as good as those for the deeper lakes.
The quantitative measures of model goodness-of-fit support these qualitative observations (Table 6).For lakes within both national parks, the NSE increased as maximum depth (Zmax) and range of observed temperatures both increased.NSE values were generally large and just a few points below their effective maximum (R 2 ) in most cases.Commonly, an NSE value of 0.5 or greater is considered to be acceptable [60].The only lake falling below this threshold was Ahmik Lake in ISRO, whose negative NSE indicated that a simple mean was a better predictor of temperature than was the model.The mean error was positive for all lakes (Table 6), indicating that the model systematically overestimated temperatures, with Ahmik again having not only the largest mean error but also the largest root mean squared error (RMSE; see also Figures 6 and 7).The clouds of points for the deeper lakes had wide ranges of values, and the regression lines (red) were nearly collinear with the blue 1:1 line.In contrast, the shallower lakes had smaller clouds of points with less elongation, and hence the regression line was less constrained to be similar to the 1:1 line.Most points fell just to the right of the blue 1:1 line, demonstrating that the model overestimated temperatures slightly but consistently.

Model Runs
As noted earlier, model output was distilled into time series of annual values from 50 years of MINLAKE model run, from 1962 to 2011 (Figures S1-S8).For each lake, figures show 50-year time series of output variables.Three output variables were selected for comparison with paleolimnological markers of lake change: epilimnion temperature, epilimnion temperature variability, and longest duration of thermocline (Figures 8-15).

Model Runs
As noted earlier, model output was distilled into time series of annual values from 50 years of MINLAKE model run, from 1962-2011 (Figures S1-S8).For each lake, figures show 50-year time series of output variables.Three output variables were selected for comparison with paleolimnological markers of lake change: epilimnion temperature, epilimnion temperature variability, and longest duration of thermocline (Figures 8-15).
Ahmik Lake Epilimnion Temperature     Cruiser Lake Epilimnion Temperature    Ek Lake Epilimnion Temperature     Harvey Lake      Epilimnion Temperature   Peary Lake Epilimnion Temperature   Richie Lake Epilimnion Temperature   Not surprisingly, the most direct relation found in the lake output to climate warming was in shallow-water temperatures during summer (JJA).Every lake showed an apparent warming at a 1m depth over time (Table 7).The loess curves indicated that shallow temperatures rose from about Siskiwit Lake Epilimnion Temperature  Not surprisingly, the most direct relation found in the lake output to climate warming was in shallow-water temperatures during summer (JJA).Every lake showed an apparent warming at a 1-m depth over time (Table 7).The loess curves indicated that shallow temperatures rose from about 1962 to 1980, remained steady from 1980 to 1992, and then continued upward from 1992 to 2011.In addition, the increase in year-to-year variability in shallow temperatures was strongly significant for all lakes (Table 7), indicating a more uncertain habitat that may favor opportunistic species.These water-temperature patterns were apparently driven by increases in air temperatures and decreased wind speed at International Falls and Duluth, the two weather stations used in the MINLAKE runs (Figures 4 and 5).Trends in JJA bottom-water temperatures were not as uniform.For well-mixed shallow lakes, bottom water temperatures should follow similar trends as those for shallow water.While this generalization held true for the shallower ISRO lakes (Ahmik and Harvey; Figures S1 and S4), it did not for the shallower VOYA lakes (Ek and Peary; Figures S3 and S6) because these two lakes were in fact deep enough to stratify for substantial periods of time, as confirmed by observed profiles (Supplementary Materials) and their geometric ratios.The deep lakes in both VOYA and ISRO generally showed no significant bottom-water temperature trends, with the exception of Little Trout (Figure S5), whose temperature appeared to decline.This counterintuitive result may have been caused by a slightly shallower thermocline during the second period, thereby reducing heat transfer to the hypolimnion.

Little Trout Lake
Trends were assessed for the formation and duration of thermoclines, here equated to temperature gradients equaling or exceeding 2 • C per meter depth.Climate warming might be expected to result in a higher frequency (total number of days) and continuous duration of these steep gradients, and this appeared to be generally true for most of the deep lakes (Little Trout, Richie, and Siskiwit, with Cruiser being the exception) (Table 7).In general, these results imply an earlier spring formation and later fall dissipation of these gradients-even though, individually, these dates did not commonly show significant trends (Figures S2, S5, S7, and S8).Climate warming, especially if accompanied by stronger winds, might also be expected to cause an increase in thermocline depth.Diatom data from Siskiwit Lake is in fact consistent with this hypothesis, suggesting a large increase (a deepening) in thermocline depth [8].In direct contrast, the model results indicated that the thermocline depths decreased slightly but significantly for three of the deep lakes (Table 7 and Supplemental Figures), which was supported by the decrease in wind velocity from 1962 to 2011 noted above.

Paleolimnological Records
Although cores were analyzed back to sediment levels that pre-date Euroamerican settlement (ca.1850-1870; [8,13,40]), here we report data in a comparable fashion to the thermal model output using downcore data from 1960 to 2010 (Figures 8-15).In most cores, sampling intervals from 1960 to 2010 represent an approximate decadal sampling resolution.The relative abundance predominant species in the cores captures biological changes in community makeup; species autecology is used in interpreting these data.The DCA axis 1 scores are reported as a measure of community compositional turnover [59].Although thermal model output treats only the time period 1960-2010, DCA scores were calculated for the entire core record (~1800-2010; [13]) but are only reported for 1960-2010.Carbon accumulation rates are also plotted for each core as a measure of historical in-lake productivity and terrestrial C sources to each lake.

Diatom Species Response
In addition to nutrient levels, pH, and conductivity, the depth of lakes and lake morphometry are strong determinants of diatom communities primarily through physicochemical structuring and habitat availability.As such, the diatom communities and changes between 1960 and 2010 in the eight lakes are well differentiated between the shallower and deeper lakes, but also differ between the two park units.The downcore distribution of select diatoms that were present at >5% abundance in two or more core depths provides an overview of biological changes that have occurred during the last 50 years in these wilderness lakes.The stratigraphic records of additional taxa mentioned in the results are available in Reference [13].
Shallower lakes of VOYA (Ek, Peary)-Sediment records from the two shallower lakes of VOYA were dominated by similar diatom species, notably the acidophilic taxon Eunotia zasuminensis and Tabellaria flocculosa IIIp (Figures 10 and 13).These taxa behaved similarly in both lakes between 1960 and 2010 with E. zasuminensis decreasing and T. flocculosa IIIp simultaneously increasing in abundance.Additionally, Staurosirella pinnata decreased in abundance in Ek Lake [13] and Discostella stelligera decreased slightly in Peary Lake.
Deeper lakes of VOYA (Cruiser, Little Trout)-The plankton of Cruiser Lake was dominated by few species, including Discostella stelligera and Asterionella formosa [13] and two morphs of Cyclotella ocellata (Figure 9).The two former taxa showed little change in abundance between 1960 and 2010, whereas the two morphs of C. ocellata showed opposite trends with the nominate form increasing in abundance in recent decades.In Little Trout Lake, the plankton comprised a more diverse diatom community of several Cyclotella species, Asterionella formosa, Aulacoseira subarctica (Figure 12), Tabellaria flocculosa IV and D. stelligera [13].Cyclotella ocellata and C. michiganiana increased in abundance since the 1980s, whereas D. stelligera had variable abundance [13].In contrast, Asterionella formosa and Aulacoseira subarctica decreased in abundance in recent decades.
Shallower lakes of ISRO (Ahmik, Harvey)-In contrast to the shallower VOYA lakes, ISRO's Ahmik and Harvey lakes were dominated by what are known as small fragilarioid diatoms (Figures 8 and 11).These species are considered ecological generalists and can live on benthic surfaces and be entrained in the water column during mixing events, such as tychoplankton.In Ahmik Lake, Staurosirella pinnata and Staurosira construens have recently increased in abundance, whereas Pseudostaurosira brevistriata decreased in abundance.Harvey Lake saw the increase of Staurosirella pinnata and an undescribed morphotype called Staurosira construens SWMN2, simultaneous with a decreasing abundance of S. construens [13] and S. elliptica.
Deeper lakes of ISRO (Siskiwit, Richie)-The two deeper lakes at ISRO have strong differences in productivity, watershed area, and morphometry.Richie Lake is mesotrophic to borderline eutrophic, and its deep waters are limited to a south basin.Much of the lake remainder is relatively shallow and, in recent years, Richie has been experiencing mid-to late-summer cyanobacterial blooms.Siskiwit Lake is much larger and deeper than Richie, is oligotrophic, and has a coldwater lake trout and cisco fishery.The plankton diatom flora of Siskiwit is dominated by cyclotelloid species [13] such as Cyclotella planktonica, C. atomus, C. michiganiana, the C. comensis-group, and Discostella stelligera (Figure 15).Patterns of change over the last 50 years include an increasing abundance of C. planktonica and C. atomus [13], and a declining abundance of D. stelligera and C. michiganiana.In contrast to Siskiwit, the plankton flora preserved in the core from Richie Lake comprises a meso-to eutrophic flora of Aulacoseira species [13], Fragilaria crotonensis, Asterionella formosa, and Tabellaria flocculosa IIIp (Figure 14).The last 50 years have seen decreased abundances of Aulacoseira subarctica, A. granulata, and T. flocculosa IIIp [13].The period of time since the 1980s has been characterized by a rapid increase in abundance of F. crotonensis and Asterionella formosa.

Diatom Community Turnover
Community turnover among diatoms based on DCA axis 1 scores varied among the eight lakes.In general, shallower lakes had lower species turnover compared to deeper lakes.For example, Ek and Peary lakes at VOYA showed only 0.22 and 0.31 SD units of turnover, respectively, between 1960 and 2010 (Figures 10 and 13).Cruiser and Little Trout lakes had DCA axis scores that ranged over 0.60 and 1.04 SD units during the same time period (Figures 9 and 12).Of note is the rapid turnover in diatom communities in Little Trout Lake between the mid-1970s and mid-1990s (Figure 12).At ISRO, a similar pattern emerged with a somewhat lower turnover in the shallower lakes, Harvey and Ahmik, with a community turnover of 0.48 and 0.025 units, respectively (Figures 8 and 11).The larger and deeper lakes of ISRO, Siskiwit and Richie, had turnovers of 0.46 and 1.14 SD units (Figures 14 and 15).Lake Richie in particular had the highest community turnover rates among the eight lakes (Figure 14); this wilderness lake has been noted to have periodic cyanobacteria blooms since 2007 [3].

Carbon Burial
Several patterns of OC burial are evident in the lakes.First, the shallower lakes of VOYA (Ek, Peary) and ISRO (Ahmik, Harvey) tend to have greater OC burial rates than the deeper lakes in the parks ( Figures 8,10,11 and 13).Second, almost all lakes show trends of increasing OC burial between 1960 and 2010 (Figures 8-15).Little Trout and Peary at VOYA are exceptions, showing little change in OC burial between 1960 and 2010 (Figures 12 and 13).All lakes show longer-term increases when records from mid-1800s to 2010 are considered [13], with the increases strongly continuing from 1960-2010 (Harvey, Siskiwit, Ek, and Cruiser; Figures 9-11 and Figure 15), plateauing before 1960 (Little Trout and Peary; Figures 12 and 13), or flattening out between 1960 and 2010 (Richie and Ahmik; Figures 8 and 14).

Discussion
Climate change impacts lakes through direct and indirect effects by shifting energy and mass inputs to lakes [61].Direct effects include the transfer of energy to lakes by heat, radiant energy, and wind, as well as inputs of mass (e.g., precipitation, nutrients) directly to the lake.Indirect effects of climate are mediated through the watershed, for example by changes in the transfer of energy to the lake by shading or by inputs of mass from the watershed (nutrients, DOC, and mineral matter, both particulate and dissolved).In many lake-rich regions, human activities in the lake or watershed may mask changes driven by climate [62][63][64].We focused on changes in wilderness lakes where direct human impact is minimal and watershed disturbances are primarily limited to natural forces.To explore potential climate impacts on wilderness boreal lakes, we paired the paleolimnological analysis of sediment cores with the thermal modeling of eight lakes in ISRO and VOYA to link historical lake thermal structure with biological response.

Model Performance
The reason the thermal structure of shallow lakes was not simulated as well as that of deep lakes may be related to less well-constrained factors in the model that have a disproportionate effect on smaller and shallower lakes.Heat loss to sediments and to outflowing water (both surficial outlets and groundwater recharge) would have a greater effect on lakes with small volumes.Lakes shallow enough to have sunlight reach their bottoms would be susceptible to small errors in the solar absorbance of the sediment surface.Lakes with small surface areas are more influenced by the wind sheltering effects of surrounding topography and vegetation, especially the effects of winds from different directions.MINLAKE allows the user to set a constant wind-sheltering coefficient to aid calibration (Table 4), but the coefficient is a general term meant to account for average or typical conditions, not for details of different wind speeds from different directions.
The proximity of the weather station used to drive the model may have influenced our results.The model runs for the VOYA lakes used data from the nearby International Falls weather station, and model fits for these four lakes were all quite good, with even the shallowest lake (Peary) having a good NSE value of 0.58.In contrast, the nearest weather station with model-ready data to the ISRO lakes was Duluth, 240 to 320 km away.Large lakes like Siskiwit may have enough thermal inertia to respond to a smoothed average of weather conditions and so may be reasonably modeled based on a somewhat distant weather station.Smaller lakes with faster responses would be more prone to mismatches due to wind and radiation differences between ISRO and the weather station in Duluth.Small and shallow lakes would have relatively large variations in water temperature and stratification characteristics during the course of each day from early morning to late afternoon.The simulated temperature using MINLAKE2012 with a daily time step is closer to the temperature at late afternoon because daily solar radiation is used to heat the lake, whereas observed profiles in these lakes may have been taken earlier or later on the day of sampling.The mismatch of times for simulated and measured temperature profiles would create larger errors in shallow lakes than in deep lakes.Nonetheless, good model results for the ISRO lakes indicate that weather data from Duluth are adequate for modeling the thermal structure, suggesting that the principal drivers (temperature, solar radiation, and wind) are reasonably similar between ISRO and Duluth.
Although other thermal models exist for lakes (e.g., [9,65]), we feel that MINLAKE2012 is appropriate and well-suited for these wilderness lakes and could be applied to other lakes and for forecasting future lake responses.We did not use the DO capabilities of MINLAKE2012 in this project.However, we recognize that changes in DO may be a primary driver of ecological change in lakes with altered thermal structure due to climate warming [66].

Climate-Lake Response
Specific MINLAKE2012 model results were extracted and analyzed to determine trends in shallow and deep water temperatures, interannual water temperature variability, and the behavior and timing of water column gradients over two time periods .The most common significant trend was increased summer shallow-water temperatures for all eight lakes; this trend was also accompanied by a significant increase in year-to-year variability in shallow water temperatures.The next significant trend was, for the deep lakes, an increase in the frequency and duration of thermal gradients equaling or exceeding 2 • C per meter.Surprisingly, there were no significant trends identified among the lakes in timing or duration of ice cover.
The modeled water-temperature and stratification patterns were coincident with increases in mean annual air temperatures and decreased wind speed at the two weather stations.Seasonally, increased winter air temperature was significant at both weather stations, but increased summer air temperature was only significant at the Duluth station.Wind speed decreased significantly for all seasons and the annual average at both meteorological stations.These meteorological trends are likely key for understanding how changes in thermal structure may have impacted lake structure and biology during the simulation period, and are consistent with other regional meteorological syntheses [2].
The most common significant trend we modeled was an increase in summer shallow-water temperatures for all lakes.Throughout the world, air temperature increases are correlated with surface water temperature increases [62,67,68] and linked to shorter ice cover, decreased cloud cover, increased summer air temperatures.Based on monitoring data for 26 lakes worldwide from 1970 to 2010, lake surface water temperatures, whole lake temperatures, and bottom water temperatures all increased [66].Other trends included higher latitude lakes warming more than low latitude lakes, and surface water temperatures warming more than deep water temperatures [66].In Europe, surface water temperatures increased up to 1.4 • C between 1960 and 2009, with greater temperature increases in larger (surface area) lakes [69].In Wisconsin (USA) lakes, observed trends included shallow water warming in all lakes, but deep water warming only in large lakes, potentially due to the enhanced wind sheltering of small lakes [70].In northeast North America, Richardson et al. [62] synthesized broad regional lake records to show that surface water temperatures are rising more rapidly than air temperatures and more rapidly than other regions; rising temperatures are accompanied by increased lake stratification; deep water temperatures are not showing consistent warming patterns among lakes; clear lakes are warming most rapidly; and lakes farther from ocean influence are warming more quickly.Broadly similar warming trends, especially in shallow water temperatures, were evident in our modeled data and are corroborated by worldwide trends.Other model output, especially warming trends in deep water temperatures, were not common to all of our model hindcasts.
Other trends in our modeled results include significantly greater inter-annual variability in summer shallow and deep water temperatures.This variability appears to be producing greater inter-annual extremes in thermal structure, which may be critical in generating regime shift lake responses to direct and indirect climate drivers [71].Unfortunately, the temporal resolution of paleolimnological proxies (diatom abundance, OC burial) are limited to five to 10 years in this study, making direct connection to inter-annual modeling results difficult.However, the lake with the greatest diatom community turnover, Lake Richie, has experienced unprecedented change likely linked to inter-annual extremes.In 2007, the lake had modeled high hypolimnetic temperatures and shallow thermocline depth.For the first time in park history, the lake experienced its first recorded noxious cyanobacterial bloom that required park managers to post warning for wilderness campers.This wilderness lake may be experiencing a threshold ecological response [71] as diatom, cyanobacterial, and fish (loss of coregonids) communities are all undergoing dramatic restructuring [3,13,72].
An increase in the frequency and duration of thermal gradients, i.e. stratification, was the other significant modeled trend in our deep lakes, except VOYA's Cruiser Lake, a trend that appears to be common among modeled and monitored lakes.Hadley et al. [9] used the heat transfer model DYRESM to model daily temperature-depth profiles from 1979 to 2009 in Ontario's Harp Lake (Canada), a small, deep (37.5 m), oligotrophic, forested Canadian Shield lake.Similar to our results, Harp Lake showed a significant increase in modeled thermal stability (Schmidt's stability) in response to increased air temperatures and decreased winds; however, in contrast to our results, increased stability in Harp Lake was not accompanied by any significant change in the timing of onset, duration, or breakdown of stratification.Winder and Schindler [73] documented increased stability in Lake Washington (King County, WA, USA) between 1962 and 2002 that was accompanied by earlier onset and longer duration of stratification.Similar trends were also shown for Lake Simcoe (Barrie, ON, Canada), where the lake was experiencing six additional days of stratification between 1980 and 2008 [74].Increased stability and duration of stratification is not limited to northern lakes, as stratification patterns have been more affected globally in deep lakes [66].

Diatom Community Shifts
Biological changes in the lakes between 1960 and 2010 differed among the shallow and deep lakes and also between the two park units.The latter difference most likely reflects differences in lake type and water chemistry between the parks [75].
Shallow lakes showed a combination of primary shifts in benthic or tychoplanktonic taxa and secondary shifts in planktonic species.In the shallow lakes of ISRO, the small benthic and tychoplanktonic fragilarioid species were the main diatoms to respond, with Staurosirella pinnata and Staurosira construens increasing in abundance and Pseudostaurosira brevistriata decreasing in abundance in Ahmik Lake, and several small Staurosirella and Staurosira taxa shifting dominance in Harvey Lake.Changes in abundance of these small fragilarioid forms are poorly constrained ecologically.Most paleoecologists treat them as an ecological group ("small fragilarioids"; [76]) and consider their overall community response.Species-level shifts likely represent changes in benthic microhabitat availability, shallow lake mixing events, and may reflect the overall trend of shallow water warming we show with thermal modeling.In VOYA shallow lakes, changes in E. zasuminensis and T. flocculosa IIIp reflect subtle increases in lake pH and nutrients as noted by Edlund et al. [13] and Hobbs et al. [3], and the decrease of Discostella stelligera in Peary Lake possibly reflect the weakening of a thermocline [8] or more rapid spring warming [77], although these trends were not identified as significant in our thermal modeling output.
Most diatom shifts in the deep lakes appear to track changes in thermal structure, as evidenced by shifts in the cyclotelloid genera Cyclotella and Discostella.The ecologies of these deep lake species have not been completely resolved; however, we generally consider the small Cyclotella and Discostella species to be members of the deep chlorophyll layer, i.e., associated with the thermocline [8,76,78].In many northern lakes, their increased abundance is considered indicative of longer and more stable periods of stratification [1,76,79].Others have also shown that the size structure within populations of these species may also reflect climate drivers [80].Their ecologies may also reflect light regime, seasonality, and nutrient supply [77,78,81].Other diatom species (Asterionella formosa Fragilaria crotonensis, and Aulacoseira subarctica) shown to change in the deep lakes are more likely responding to the length of mixing and nutrient supplies [3,82].A better understanding of species-level ecology achieved through modern sampling and experimentation (see [8,78]) would provide a more mechanistic understanding of how diatom communities in boreal lakes respond to direct and indirect climate drivers.
For the deep lakes at VOYA, changes in the abundance of the deep chlorophyll layer taxa support modeled changes in the thermal structure.The increased abundance of small Cyclotella species (C.ocellata in both Cruiser and Little Trout and C. michiganiana in Little Trout) support the thermal modeling results that show Cruiser Lake to be trending toward longer periods of stability and deepening of the thermocline.Little Trout Lake has exhibited a slight shallowing of the thermocline but an extended period of stability in recent decades.At ISRO, the two deep lakes, Siskiwit and Richie, differ in their productivity, basin morphometry, and ecology [75].Siskiwit is oligotrophic, supports a cold-water fishery, and is much deeper.Richie Lake is meso-to borderline eutrophic and has only a single deep basin with much of the lake being very shallow.In Siskiwit, the increased abundance of C. planktonica and decreased abundance of Discostella stelligera has been used in models to suggest that the thermocline in Siskiwit Lake has deepened significantly in recent decades [8], possibly in response to stronger winds around Lake Superior [83], a conclusion that was not well supported by the thermal modeling.While the trend analysis indicates a significant change in thermocline depth, the graphical presentation suggests only a subtle change that appears to be a shallowing, and not a deepening, of the thermocline.We are not able to readily resolve these conflicting results.It may be a reflection of the meteorological data used for modeling ISRO lakes being from the Duluth (Minnesota) weather station and not accurately capturing conditions on an island in Lake Superior that is over 280 km away.Nonetheless, the meteorological data were good enough to result in a very good fit of the model to observed values: of all the lakes, Siskiwit had the lowest RMSE and the second highest NSE and R 2 values (Table 5).Alternatively, the mixing depth model used by Saros et al. [8] is based on observations and experimental data from Rocky Mountain lakes and may be limited in its applicability to Siskiwit Lake.Detailed species-level ecological data on the dominant diatom plankters are needed to fully resolve the conflicting interpretations.
In contrast to Siskiwit Lake, Richie Lake (ISRO) showed changes in many diatom species that suggest increased nutrient availability and changes in spring mixing.Its meso-to eutrophic diatom flora of Aulacoseira species, Fragilaria crotonensis, Asterionella formosa, and Tabellaria flocculosa IIIp showed decreased abundances of Aulacoseira subarctica, A. granulata, and T. flocculosa IIIp during the modeled period, coincident with the increased abundance of F. crotonensis and Asterionella formosa.The decreased abundance of Aulacoseira species is a widespread indicator of shorter spring mixing [1,77]; these heavily silicified taxa are dependent on strong mixing to remain in the water column.The ecologies of F. crotonensis and Asterionella formosa, which have recently increased in abundance in Richie Lake, show strong seasonality in relation to temperature and mixing [77,84], but they are also strong indicators of both nutrient enrichment and warming [3,85].In particular, the two species have shown broad regional response to nitrogen (N) addition [82].Two changes in Richie Lake may be driving shifting nutrient availability.As noted below, increased export of dissolved inorganic and organic nitrogen (DIN and DON, respectively) from the watersheds at ISRO is occurring as a response to climate (especially winter) warming [86].MINLAKE modeling results also show greater inter-annual variation in many physical parameters (ice duration, duration of stratification, depth of thermocline) that would affect internal nutrient cycling in Richie Lake.The increased abundance of F. crotonensis and A. formosa has continued in recent sediment samples taken since the initial cores were collected and support the notion that a combination of changes in the mixing and nutrient environment (greater availability) of Lake Richie are at play in driving recent unprecedented change in this lake, including periodic noxious cyanobacteria blooms since 2007 [3,87].
Diatom community turnover is accentuated in the deeper lakes compared to the shallow lake sediment records of VOYA and ISRO.The greatest turnover is noted in ISRO's Richie Lake, where noxious cyanobacterial blooms have been commonplace since the mid-2000s [3].The less productive deep lakes (Cruiser, Siskiwit, and Little Trout) also show higher community turnover rates than most of the shallower lakes.Because of the physical thermal environment (stratification, seasonal vs. polymixis) associated with deeper lakes, their diatom communities have been identified as the strongest to respond to direct climate drivers, whereas shallow lakes show more nuanced changes in community response, likely controlled by indirect climate drivers [1,78,88].

Carbon Burial
Changes in burial rates of organic carbon (OC) in lake sediments reflect shifts in both in-lake productivity and external or watershed inputs of C [2].Among our modeled lakes, paleolimnological records show that shallow lakes have greater OC burial rates than deep lakes, and most lakes show increasing OC burial during 1960-2010, a trend that has continued since the mid-1800s [13].Subtle shifts in accumulation rates of biogenic silica, a proxy for diatom productivity, hint that siliceous algae may be becoming more productive in some lakes in the last few decades (Cruiser, Ek, Richie, Harvey; [13]).
Heathcote et al. [2] showed that rates of OC burial have been increasing in the last ~100 years in northern lakes around the world.Many of these lakes are similar to our study lakes, i.e., they are not undergoing dramatic land-use change or development [26,27].In similar lakes where land use is changing, the effects of climate drivers are often masked [62,89].Warming temperatures could not fully explain the increased burial rates; N deposition and cycling were suggested to be important drivers of recent OC burial trends.Hobbs et al. [3] summarized monitored N deposition trends against δ 15 N sediment records from our same model lakes.At VOYA, N deposition has declined since the 1990s, and at ISRO there is no trend in deposition; however, there are decreasing precipitation trends (especially in winter), suggesting that the N concentration in precipitation is increasing [86].However, at ISRO, the long-term N deposition monitoring that has been in place since 1985 is coupled with the monitoring of watershed output of DIN, DON, and DOC.Increased export of DIN, DON, and DOC from the watershed is occurring in spite of the lack of trends in atmospheric deposition as warmer winters and less snow cover increase mineralization and export of watershed N [86].

Conclusions
Information from this modeling exercise indicated that lakes and in-lake habitats in the northern lakes region may be responsive to climate drivers that include increased mean annual temperature and decreased wind speed.Responses of our study lakes based on the thermal model MINLAKE2012 from 1960 to 2010 included the warming of shallow waters in all lake types and increased frequency and duration of thermal gradients equaling or exceeding 2 • C per meter in deeper lakes.Biological and biogeochemical responses preserved in paleolimnological records including diatom species abundances, diatom community turnover, and organic carbon burial were compared to model hindcasts of lake temperatures, stratification, and mixing regimes.Increased rates of carbon burial were common among lakes and continued a trend that began in the late 1800s for most lakes.Diatom community turnover was greater for deep lakes compared to shallow lakes, likely reflecting stronger stratification, less intense spring mixing, and possible nutrient shifts from internal and watershed loading.Diatom species that are responding at timescales similar to modeled thermal trends have autecologies that support changes in nutrient availability and modeled trends in mixing, stratification, and lake warming.

Figure 4 .
Figure 4. Annual and summer weather patterns at Duluth, Minnesota, 1962-2011.Left panels are the annual means of air temperature (°C), wind speed (m/s), and solar radiation (Langley).Right panels are the summer means of the same parameters.

Figure 4 .
Figure 4. Annual and summer weather patterns at Duluth, Minnesota, 1962-2011.Left panels are the annual means of air temperature ( • C), wind speed (m/s), and solar radiation (Langley).Right panels are the summer means of the same parameters.

Figure 5 .
Figure 5. Annual and summer weather patterns at International Falls, Minnesota, 1962-2011.Left panels are the annual means of air temperature (°C), wind speed (m/s), and solar radiation (Langley).Right panels are the summer means of the same parameters.

Figure 5 .
Figure 5. Annual and summer weather patterns at International Falls, Minnesota, 1962-2011.Left panels are the annual means of air temperature ( • C), wind speed (m/s), and solar radiation (Langley).Right panels are the summer means of the same parameters.

Figure 6 .
Figure 6.Observed versus simulated temperatures (T), showing regression line (red) through the cloud of observed values and the 1:1 line (blue) for four lakes in Voyageurs National Park.(a) Shallow Ek Lake; (b) shallow Peary Lake; (c) deep Cruiser Lake; (d) deep Little Trout Lake.

Figure 7 .
Figure 7. Observed versus simulated temperatures (T), showing regression line (red) through the cloud of observed values and the 1:1 line (blue) for four lakes in Isle Royale National Park.(a) Shallow Ahmik Lake; (b) shallow Harvey Lake; (c) deep Richie Lake; (d) deep Siskiwit Lake.

Figure 6 .
Figure 6.Observed versus simulated temperatures (T), showing regression line (red) through the cloud of observed values and the 1:1 line (blue) for four lakes in Voyageurs National Park.(a) Shallow Ek Lake; (b) shallow Peary Lake; (c) deep Cruiser Lake; (d) deep Little Trout Lake.

Figure 6 .
Figure 6.Observed versus simulated temperatures (T), showing regression line (red) through the cloud of observed values and the 1:1 line (blue) for four lakes in Voyageurs National Park.(a) Shallow Ek Lake; (b) shallow Peary Lake; (c) deep Cruiser Lake; (d) deep Little Trout Lake.

Figure 7 .
Figure 7. Observed versus simulated temperatures (T), showing regression line (red) through the cloud of observed values and the 1:1 line (blue) for four lakes in Isle Royale National Park.(a) Shallow Ahmik Lake; (b) shallow Harvey Lake; (c) deep Richie Lake; (d) deep Siskiwit Lake.

Figure 7 .
Figure 7. Observed versus simulated temperatures (T), showing regression line (red) through the cloud of observed values and the 1:1 line (blue) for four lakes in Isle Royale National Park.(a) Shallow Ahmik Lake; (b) shallow Harvey Lake; (c) deep Richie Lake; (d) deep Siskiwit Lake.

Water 2017, x, x 15 of 34 Figure 7 .
Figure 7. Observed versus simulated temperatures (T), showing regression line (red) through the cloud of observed values and the 1:1 line (blue) for four lakes in Isle Royale National Park.(a) Shallow Ahmik Lake; (b) shallow Harvey Lake; (c) deep Richie Lake; (d) deep Siskiwit Lake.

Figure 8 .
figures and legends toge

Figure 9 .
Figure 9. Model results and paleolimnological data for Cruiser Lake, Voyageurs National Park, 1960-2011.(a) Epilimnion temperature (°C); (b) nine-year running average of epilimnion temperature (SD, figures and le Comment [M and coauthor used in Figs 8 been cut off i to be aware o

Table 1 .
Study lake locations and assembled background data 1 .

Table 2 .
Definitions of general lake categories based on maximum depth, area, and Secchi-disk depth 1 .

Table 4 .
MINLAKE-specific parameters for study lakes1.Lake parameters for Voyageurs National Park (VOYA) and Isle Royale National Park (ISRO) lakes.

Table 5 .
Significance of change over time from the first halfto the second halfof the data set for selected climate variables 1 .

Table 5 .
Significance Mann-Whitney U Test used to determine if the early period is significantly different from the later period.Dark cells are significant at the 0.05 level or less; shaded cells are significant at the 0.1 level.Abbreviations: VOYA, Voyageurs National Park; ISRO, Isle Royale National Park; °C, degrees Celsius; m/s, meters per second; mm, millimeters; NA, not applicable.
of change over time from the first halfto the second halfof the data set for selected climate variables1.Notes:1Mann-Whitney U Test used to determine if the early period is significantly different from the later period.Dark cells are significant at the 0.05 level or less; shaded cells are significant at the 0.1 level.Abbreviations: VOYA, Voyageurs National Park; ISRO, Isle Royale National Park; • C, degrees Celsius; m/s, meters per second; mm, millimeters; NA, not applicable. 1

Table 6 .
Goodness-of-fit statistics relative to lake depths and range of observed temperatures 1 .
Notes: 1 NSE, Nash-Sutcliffe coefficient of efficiency; R 2 , coefficient of determination, the fraction of variance explained by linear regression; RMSE, root mean square error, essentially the standard deviation of the errors; Zmax, maximum lake depth; N, number of observations compared against simulated values; Min, minimum; Max, maximum; VOYA, Voyageurs National Park; ISRO, Isle Royale National Park.

Table 6 .
Goodness-of-fit statistics relative to lake depths and range of observed temperatures 1 .
Notes: 1 NSE, Nash-Sutcliffe coefficient of efficiency; R 2 , coefficient of determination, the fraction of variance explained by linear regression; RMSE, root mean square error, essentially the standard deviation of the errors; Zmax, maximum lake depth; N, number of observations compared against simulated values; Min, minimum; Max, maximum; VOYA, Voyageurs National Park; ISRO, Isle Royale National Park.

Table 6 .
Goodness-of-fit statistics relative to lake depths and range of observed temperatures 1 .
Notes: 1 NSE, Nash-Sutcliffe coefficient of efficiency; R 2 , coefficient of determination, the fraction of variance explained by linear regression; RMSE, root mean square error, essentially the standard deviation of the errors; Zmax, maximum lake depth; N, number of observations compared against simulated values; Min, minimum; Max, maximum; VOYA, Voyageurs National Park; ISRO, Isle Royale National Park.

Table 7 .
Significance of change over time from the first halfto the second halfof MINLAKE runs for selected lake-thermal variables1.Mann-Whitney U Test used to determine if the early period is significantly different from the later period.Dark bold cells are significant at the 0.05 level; shaded cells are significant at the 0.1 level.Summary data for each lake are given at the top of the table, to help search for patterns.Abbreviations: VOYA, Voyageurs National Park; ISRO, Isle Royale National Park; m, meters; ha, hectares; µg/L, micrograms per liter; mg/L, milligrams per liter.