Transparency , Geomorphology and Mixing Regime Explain Variability in Trends in Lake Temperature and Stratification across Northeastern North America ( 1975 – 2014 )

Lake surface water temperatures are warming worldwide, raising concerns about the future integrity of valuable lake ecosystem services. In contrast to surface water temperatures, we know far less about what is happening to water temperature beneath the surface, where most organisms live. Moreover, we know little about which characteristics make lakes more or less sensitive to climate change and other environmental stressors. We examined changes in lake thermal structure for 231 lakes across northeastern North America (NENA), a region with an exceptionally high density of lakes. We determined how lake thermal structure has changed in recent decades (1975–2012) and assessed which lake characteristics are related to changes in lake thermal structure. In general, NENA lakes had increasing near-surface temperatures and thermal stratification strength. On average, changes in Water 2017, 9, 442; doi:10.3390/w9060442 www.mdpi.com/journal/water Water 2017, 9, 442 2 of 22 deepwater temperatures for the 231 lakes were not significantly different than zero, but individually, half of the lakes experienced warming and half cooling deepwater temperature through time. More transparent lakes (Secchi transparency >5 m) tended to have higher near-surface warming and greater increases in strength of thermal stratification than less transparent lakes. Whole-lake warming was greatest in polymictic lakes, where frequent summer mixing distributed heat throughout the water column. Lakes often function as important sentinels of climate change, but lake characteristics within and across regions modify the magnitude of the signal with important implications for lake biology, ecology and chemistry.


Introduction
Northeastern North America (NENA) has an exceptionally high density of lakes in comparison to the rest of the world [1].These lakes provide fundamental ecosystem services to the people and biota of the region and contribute substantially to Earth's aquatic biodiversity (e.g., [2]).Yet, much like other regions on the planet, NENA is experiencing a changing climate that results in warmer and, in NENA's case, wetter conditions relative to other locations in the USA [3].Across the NENA region between 1895 and 2011, we have seen air temperature rising by 0.09 • C decade −1 with acceleration of warming in more recent decades, 10 mm decade −1 increases in average annual precipitation and a 70% increase in extreme storm events since 1958 [3,4].Climate change is a particular threat to lakes because these systems are the lowest points on the landscape, effectively integrating climate and environmental forcing patterns from within their water-and air-sheds [5][6][7].NENA lakes may be changing faster than those in other mid-latitude regions of the world due to amplified climate forcing, i.e., warmer and wetter conditions [8].
Despite the significance of lakes in the NENA region and the potential impacts of warming, Torbick et al. [9] noted that relatively few systematic and long-term in situ lake temperature trend datasets exist for the area.Torbick et al. [9] used satellite-derived data to study NENA lake skin temperatures from 1984 to 2014, finding that lake surface temperatures increased more rapidly than air temperatures during this period.These findings are consistent with other regional-and global-scale studies that reveal that lake surface temperatures are warming at a rate that is faster than expected given surrounding air temperature warming [10][11][12][13][14].In addition to understanding surface water temperature trends, however, we need to determine if there is evidence of long-term changes in lake thermal structure and deepwater temperatures in the NENA region because these thermal habitats are critical to aquatic biota.
Globally, surface water temperatures are warming in response to climate drivers, primarily air temperature warming [5,12,15], and this has been documented across wide geographic regions over the past several decades [12,13,16].However, individual lake surface water warming trends (LSWT) can be quite variable, with some lakes showing stronger trends than others, and a few showing cooling trends in surface water temperature [12].Geographic heterogeneity in surface warming suggests that local or regional characteristics contribute to variability in lake thermal structure trends [12,13,16].
In contrast to surface water temperature trends, we know far less about how deepwater (i.e., hypolimnion) temperatures and stratification patterns have changed through time and what drives these changes.In general, deepwater temperatures are related to spring weather and ice-off phenology in dimictic lakes [17].Few studies have provided insight into regional patterns and variability of deepwater temperature changes (see [16,18]).Winslow et al. [18] hypothesized that wind sheltering in small lakes may reduce mixing and account for reduced warming in deeper waters compared to large lakes.The changes in thermal stratification patterns over time are likely to be highly variable across lakes.For example, effects of climate change on lake stratification varied with lake depth, surface area and water clarity in thermal simulation models [19] and varied with morphometry and average lake temperature in long-term monitoring records [16].Changes in lake stratification patterns are also likely region-specific, since the strength of stratification is largely related to average lake temperature [16], which will vary with latitude.
This study is the first systematic study within NENA aimed at understanding the combined effects of lake characteristics and climate change on long-term deepwater temperatures and lake stratification.We use a macrosystems approach [20] to study spatio-temporal lake thermal structure across this region in order to understand whether differential warming is being transferred below the lake surface.We assembled data on 231 study lakes spread geographically across NENA (Figure 1), with a range of mixing regimes (i.e., dimictic to polymictic); they also vary in transparency, trophic status and geomorphometry.Using this macroecological approach across a region so rich in lakes allows us to determine which types of lake (i.e., across a gradient in mixing regimes, sizes, water clarity and geomorphometry) may be most vulnerable to rapid rates of warming.Furthermore, we augment the recent NENA skin temperature study [9] by using a regionally-dense dataset of in situ long-term temperature profile records to examine surface temperature, deepwater temperature and thermal stratification changes.
Water 2017, 9, 442 3 of 22 stratification patterns are also likely region-specific, since the strength of stratification is largely related to average lake temperature [16], which will vary with latitude.This study is the first systematic study within NENA aimed at understanding the combined effects of lake characteristics and climate change on long-term deepwater temperatures and lake stratification.We use a macrosystems approach [20] to study spatio-temporal lake thermal structure across this region in order to understand whether differential warming is being transferred below the lake surface.We assembled data on 231 study lakes spread geographically across NENA (Figure 1), with a range of mixing regimes (i.e., dimictic to polymictic); they also vary in transparency, trophic status and geomorphometry.Using this macroecological approach across a region so rich in lakes allows us to determine which types of lake (i.e., across a gradient in mixing regimes, sizes, water clarity and geomorphometry) may be most vulnerable to rapid rates of warming.Furthermore, we augment the recent NENA skin temperature study [9] by using a regionally-dense dataset of in situ long-term temperature profile records to examine surface temperature, deepwater temperature and thermal stratification changes.We posed two overarching questions.First, to what extent has the thermal structure of lakes changed in NENA in recent decades?Second, which lake types, based on geographic, morphologic, trophic state or mixing regimes, are more vulnerable to rapid rates of climate warming?We hypothesized that clear, deep, large, dimictic lakes would respond to climate change with the most rapid rates of surface and deepwater warming [21,22].Further, given that lakes are sentinels of their surroundings, we expected that inland lakes may be more subject to extreme surface and deepwater temperature changes than lakes closer to the Atlantic coast due to the moderating influence of the ocean on local climate.We posed two overarching questions.First, to what extent has the thermal structure of lakes changed in NENA in recent decades?Second, which lake types, based on geographic, morphologic, trophic state or mixing regimes, are more vulnerable to rapid rates of climate warming?We hypothesized that clear, deep, large, dimictic lakes would respond to climate change with the most rapid rates of surface and deepwater warming [21,22].Further, given that lakes are sentinels of their surroundings, we expected that inland lakes may be more subject to extreme surface and deepwater temperature changes than lakes closer to the Atlantic coast due to the moderating influence of the ocean on local climate.

Data Acquisition and Quality Control
We obtained lake thermal profile data from various sources (e.g., state and provincial agencies, research stations) throughout the temperate deciduous forest region of northeastern North America, from Maine to Pennsylvania, USA, and southern Ontario, CA (Figure 1).For each of the 231 lakes, we collected the following descriptor data, if available: latitude (lat), longitude (Long), elevation (Elev), mean and maximum lake depth (Depth, MaxD), lake surface area (Sarea), Secchi disk depth (Secc), total phosphorus (TP), total nitrogen (TN), dissolved organic carbon concentration (DOC) and chlorophyll-a (Chl a).When multiple values existed for a single lake, variables were averaged over the entire dataset to produce one descriptor value for each lake.Additionally, we calculated two more descriptor variables: distance to Atlantic coast (Dist.)and fetch (Fetch).Dist.was measured as the minimum straight-line Euclidean distance (km) from each lake's centroid to the eastern coastline of North America, and fetch measures the maximum length of open water for wind to act on the surface of the lake (i.e., the longest unobstructed distance of a lake).We performed these calculations in ArcGIS with an equidistant conic projection ( [23], functions DIST and BOUNDING).
For all lakes, we carefully examined each temperature profile from all dates available for anomalous values.All extreme values that were deemed to be errors (>40 • C or <−5 • C) were removed, and all profiles with fewer than two depths of temperature readings were removed.To identify physically impossible profiles and data recording or entry errors, we examined each profile for monotonically-increasing water densities [24] from the surface to the deepest measurement and iteratively removed non-monotonically increasing densities.Following this automated process, we visually examined profiles using heat maps for any additional anomalous values, which were removed.Unless otherwise noted, all data analyses and data management were conducted in the R statistical environment [25].

Annual Profile Selection
Most of the temperature profiles were collected manually at a relatively low frequency, ~1-2-times per month.For each lake and each year, we typically had monthly profiles during the open-water period; however, there was a range from a minimum of one profile per year up to daily profiles.Two lakes had high-frequency temperature sensors that recorded every 10 min, but only for a recent subset of the years during our study period.Given the variability of profile frequency across lakes, we decided to select a single profile for each lake during peak summer stratification.This provided the most relevant period for assessing trends in lake stratification and allowed us to maximize the size of the dataset (data using this method were available for more than 200 lakes).This method is generally robust in long-term analyses of lake temperature trends, but does not give insight into short-term responses to weather patterns or changes in specific phenological events that would be possible with higher resolution data.
To determine the timing of peak summer stratification, we calculated the maximum Brunt-Väisälä buoyancy frequency at the seasonal thermocline (max.buoyancy frequency, with units of cycles per hour, cph) for each profile using the following equation: where g is gravitational acceleration (9.81 m s −1 ), ρ is water density (kg m −3 ) and z is depth (m), using the R software package "rLakeAnalyzer" [26].For each year, we found the ordinal day with the max.buoyancy frequency; this approximates timing of the strongest stratification.From those summer peak stratification days, we used the median to identify the ordinal day of maximum thermal stratification for each lake.For this region, peak stratification was typically in late July to early August.
For each lake, we selected one profile that was closest to this day within ±30 days.If there were no profiles within ±30 days for a given year, no profile was selected for that year.For each profile, we down-sampled the data to 1-m depth increments by averaging all temperature measurements ±0.5 m around the designated meter interval (e.g., <0.5 for 0 m deep and 0.5-1.5 for 1 m deep).This accounted for some profiling or depth measurements that resulted in temperature measurements at 0.9 m or 1.1 m.Finally, we truncated the bottom depth to be the deepest profile depth for each lake with at least 80% of the years' profiles having a temperature measurement at that depth ("cut-off depth"; Figure S1).

Lake Thermal Metrics
We calculated six metrics of lake thermal structure from each selected profile for each lake, which included three metrics of lake temperature and three metrics of the strength of thermal stratification.Near-surface temperature was defined as the average temperature readings between 1 m and 2 m (inclusive).This depth range was our attempt to exclude temporary thermoclines that can form at the surface in response to diurnal warming rather than as a response to longer-term changes.Deepwater temperature was defined as the average temperature from 1 m to 2 m above the lake-specific cut-off depth.Mean water column temperature was defined as the depth-weighted mean temperature of the 1-m resolution temperature profile from surface to the cut-off depth.
We calculated three metrics to capture the strength of stratification for each temperature profile.Temperature difference was defined as the near-surface minus deepwater temperature.We calculated the max.buoyancy frequency for each profile as above (Equation ( 1)).The density gradient was defined as the absolute value of the density of the near-surface water temperature minus the density of the deepwater temperature.We used the absolute value so the values of the density gradient were consistent with those for temperature difference and max.buoyancy frequency (i.e., larger positive values mean stronger thermal stratification).We note that most lakes have 1-m temperature measurement resolution, but there are differences in the vertical resolution of temperature measurements between the lakes.While the differences in resolution will affect the density gradient, the buoyancy frequency takes into account different depths between temperature measurements.
Near-surface temperature was calculated for all lakes, regardless of cut-off depth.Very shallow lakes were potentially problematic with overlap for near-surface and deepwater temperature calculations and for metrics of thermal stratification; therefore, for lakes 2 m or shallower, only near-surface temperature was calculated.For lakes with a cut-off depth between 3 m and 4 m, only near-surface temperature and mean water column temperature were calculated.For lakes with a cut-off depth greater than 4 m, all six metrics of thermal structure were calculated.For most lakes, the cut-off depth was >70% of maximum depth (Figure S1); for lakes with low cut-off depths relative to maximum depth, we considered only near-surface measurements in trend analysis.

Data Cohorts
Because of the variation across lakes in the span of data, we chose to analyze trends across two time periods, one beginning in 1975, and the other beginning in 1985.To be included in the cohort, each lake's thermal metrics were required to have at least one data point occurring between 1975 and 1980 (the 1975 cohort) or between 1985 and 1990 (the 1985 cohort), and at least one data point occurring between 2009 and 2014 (both cohorts).Within these required start and end years, a minimum of 50% of years spanning at least 15 years of data were required.In both cohorts, >50% of lakes had 80% or more years with data; the average was ~80% of years with data across all metrics.Lakes were included if at least one of the six response metrics met the completion criteria described above.Eight lakes had multiple data sources or multiple basins.For those, we selected (1) the primary/main basin, or (2) the most complete record, or (3) the deepest source/site; in all eight cases, all three were consistent.A total of 85 lakes met our data criteria for the 1975 cohort and 226 for the 1985 cohort.All lakes in the 1975 cohort were also included in the 1985 cohort, with the exception of 5 lakes.

Temporal Trends in Thermal Metrics
After calculating the thermal metrics for each summer temperature profile for the qualifying lakes, we assessed temporal trends using the Theil-Sen estimator (Sen's slope) for each lake thermal metric using the R package "wq" [27].We applied one-sample t-tests to each cohort's thermal metric slopes compared to a mean slope of zero, using a significance level of α = 0.05, to determine if the thermal metrics were, on average, changing over time.We used two-sample t-tests to compare the distribution of the slopes of each thermal metric between the two cohorts, also using a significance level of α = 0.05.All data were either normally distributed (Shapiro-Wilk test) or had a sufficiently large sample size (n > 187).

Air Temperature Data
Air temperature data were collected from the Climate Research Unit (CRU), University of East Anglia (TS Version 3.23).This dataset interpolated monthly mean air temperature values from land-based meteorological stations across the world spanning 1901-2014 at 0.5 • × 0.5 • globally-gridded resolution.For each cohort, we created a minimum convex polygon around all lakes included in the respective cohort.We calculated an annually-averaged mean air temperature for each gridded centroid from the CRU dataset that was within this polygon, spanning from 1975 to 2014 (1975 cohort) and 1985 to 2014 (1985 cohort).For each centroid's time series, we calculated the temporal rate of change as Sen's slope and took the average and standard error across all centroids as the regional air temperature rate of change.

Comparison with Air and Lake Near-Surface Temperature Trends
To compare our near-surface water temperature trends with regional air temperature trends, we applied two-sample t-tests comparing the near-surface temperature trends from each cohort to the respective air temperature slopes (one for the 1975 cohort and one for the 1985 cohort).We acknowledge that this approach is limited because regional air temperature trends were based on interpolated data; hence, variance in mean air temperatures will be underestimated, but differences between air and water temperature trends would still be evident with accurate uncertainty estimates [28].To compare to trends reported in other global studies, we used two-sample t-tests to compare to global studies where all individual LSWT data were available in the published study and one-sample t-tests where trend data were not available.

Median Buoyancy Frequency
From the selected annual profiles used for this analysis, we calculated median buoyancy frequency (MBF) as a metric for the mixing regime of each lake.Data providers indicated the mixing regime for a subset of lakes; the majority were dimictic (n = 87) with several polymictic lake (n = 8) in the 1985 cohort.Using a two-sample t-test, we compared the MBF between lakes identified as polymictic or dimictic to see if MBF was an appropriate proxy for mixing regime.Following that analysis, we ran a series of correlation analyses between the MBF and all 6 thermal metric slopes for each cohort.We used a family-wise α = 0.05 with Šidák correction for individual comparisons across each cohort.

Boosted Regression Tree Analysis
We built boosted regression tree (BRT) statistical models to define relationships between environmental explanatory variables (Table 1) and each of the 6 lake thermal metrics in our study.The BRT method is a combination of regression decision trees and a boosting algorithm through machine learning that can account for non-linear relationships between explanatory and response variables and can accommodate both categorical and numeric variables [29].We combined the two cohorts into a single dataset where every lake had only one unique value.We picked the longer time series in cases where lakes were represented in both cohorts resulting in a dataset with 231 unique lakes."Cohort" was added as a categorical explanatory variable in the BRT model.For the other explanatory variables, correlations existed between multiple sets of explanatory variables.For example, variables that reflect lake transparency and productivity were correlated.Total phosphorus and chlorophyll a concentrations were each negatively correlated with Secchi depth (r = −0.63,r = −0.57,respectively, p < 0.01 for 1975 cohort) and positively correlated with each other (r = 0.77, p < 0.01).To avoid multicollinearity in our model, we included latitude, elevation, distance to Atlantic coast, fetch and Secchi depth as explanatory variables in each model while excluding other variables (Table 1; see Figure S2 for principal component analysis of driver variables).Secchi depth was positively correlated with mean depth (r = 0.57, p < 0.001).Therefore, we included a depth index (mean depth 1.5 /Secchi depth) that was independent of Secchi depth [19].This depth index was not correlated with Secchi depth (r = 0.03, p = 0.62) and represents the effect of light penetration into water and resistance to mixing related to depth and independent of water clarity.BRT analysis was completed using the "dismo" [30] and "gbm" [31] packages in R. All boosted regression trees were fit with a model complexity value of 4 and a learning rate of 0.001.We used the "gbm" package to automatically assess the optimal number of trees included in each BRT model.Final fitted models were assessed using root mean squared error (RMSE) and the Pearson's correlation coefficient between observed trends in each of the 6 thermal metrics and those modeled by the BRT.We estimated the relative influence of each explanatory variable using the "gbm" package [32].The relative influence is based on the number of times a variable was selected for splits in the regression trees weighted by the squared improvement to the model resulting from the split averaged over all trees.The relative influence value for each variable is divided by the sum of relative influence values across all explanatory variables such that higher numbers indicate stronger influence on the response variable.
We present plots of the relationships between each response variable and explanatory variable with the largest relative influence from the boosted regression tree analysis (e.g., [33]).The plots include modeled slope estimates from the BRT ("model fits") where observed values of all explanatory variables were used to make estimations.The plots also include the modeled BRT slope estimates for each observed explanatory variable separately, where all other explanatory variables were held at their median value from the dataset ("single predictor fits").

Meta-Analysis
We surveyed primary literature databases, including Google Scholar and Web of Science, for lake temperature change rates using key words "lake*" and "climate change", as well as backward and forward literature searching in papers with trends.From the literature identified, we selected papers with lake surface water temperature change trends.We recorded trends, timing of the study, if the survey was for many lakes across the globe, a set of regional lakes or a single lake and air temperature trends, if reported.

Results
In the NENA lakes, transparency and trophic status ranged from clear-water, low nutrient lakes to less transparent, highly productive lakes.The range of lake types included dimictic and polymictic lakes and large (>200 km 2 ) to very small lakes (<0.02 km 2 ) (Table 1; Figure S1; see Table S1 for a full list of lakes and characteristics).Lake characteristics are similar for both the 1975 cohort (n = 85) and the 1985 cohort of lakes (n = 226, Table 1).
Across all lakes in both the 1975 and 1985 cohorts of lakes, near-surface and mean water column temperatures warmed, and the strength of thermal stratification increased (Figure 2).For all thermal metric trends, there were no significant differences between the 1975 and 1985 cohorts (two-sample t-test, all p > 0.05; Figure 2).Near-surface temperatures warmed in both the 1975 and 1985 cohorts (one-sample t-tests: both cohorts p < 0.001) at mean rates of 0.054 • C y −1 and 0.052 • C y −1 , respectively.However, deepwater temperatures as a whole did not change significantly in either cohort, with mean rates of change of −0.008 • C y −1 and −0.002 • C y −1 (p = 0.18, p = 0.62, respectively).The majority of lakes in the 1975 and 1985 cohorts had some degree of near-surface warming trends (89% and 81%, respectively).However, near-surface warming trends were not always accompanied by corresponding deepwater warming, with only 48% of lakes exhibiting both near-surface and deepwater warming, similar in both cohorts.A significant positive, but weak relationship (16% and 6% of the variance explained for 1975 and 1985 cohorts, respectively) between near-surface and deepwater temperature trends was evident: generally, lakes with the most rapid warming trends in near-surface waters also had the most rapid warming in deepwater temperatures (Figure S3).In both cohorts, the near-surface and deepwater temperatures cooled in about ~10% of lakes.A small percentage of lakes (<5%) had near-surface cooling and deepwater warming.
if the survey was for many lakes across the globe, a set of regional lakes or a single lake and air temperature trends, if reported.

Results
In the NENA lakes, transparency and trophic status ranged from clear-water, low nutrient lakes to less transparent, highly productive lakes.The range of lake types included dimictic and polymictic lakes and large (>200 km 2 ) to very small lakes (<0.02 km 2 ) (Table 1; Figure S1; see Table S1 for a full list of lakes and characteristics).Lake characteristics are similar for both the 1975 cohort (n = 85) and the 1985 cohort of lakes (n = 226, Table 1).
Across all lakes in both the 1975 and 1985 cohorts of lakes, near-surface and mean water column temperatures warmed, and the strength of thermal stratification increased (Figure 2).For all thermal metric trends, there were no significant differences between the 1975 and 1985 cohorts (two-sample t-test, all p > 0.05; Figure 2).Near-surface temperatures warmed in both the 1975 and 1985 cohorts (one-sample t-tests: both cohorts p < 0.001) at mean rates of 0.054 °C y −1 and 0.052 °C y −1 , respectively.However, deepwater temperatures as a whole did not change significantly in either cohort, with mean rates of change of −0.008 °C y −1 and −0.002 °C y −1 (p = 0.18, p = 0.62, respectively).The majority of lakes in the 1975 and 1985 cohorts had some degree of near-surface warming trends (89% and 81%, respectively).However, near-surface warming trends were not always accompanied by corresponding deepwater warming, with only 48% of lakes exhibiting both near-surface and deepwater warming, similar in both cohorts.A significant positive, but weak relationship (16% and 6% of the variance explained for 1975 and 1985 cohorts, respectively) between near-surface and deepwater temperature trends was evident: generally, lakes with the most rapid warming trends in near-surface waters also had the most rapid warming in deepwater temperatures (Figure S3).In both cohorts, the near-surface and deepwater temperatures cooled in about ~10% of lakes.A small percentage of lakes (<5%) had near-surface cooling and deepwater warming.All three metrics of thermal stratification increased for both the 1975 and 1985 cohorts (p < 0.001 for all), including: temperature difference (0.062 and 0.048 • C y −1 ); max.buoyancy frequency (1.61 and 1.35 cph y −1 ); and density gradient (0.002 and 0.001 kg m −3 y −1 ).The majority of lakes had increasing stratification indicated by temperature difference, max.buoyancy frequency and density gradient (87-88% for the 1975 cohort; 72-80% for the 1985 cohort).Mean water column temperature significantly increased in both cohorts, but to a lesser degree than near-surface temperatures (0.008 • C y −1 for both cohorts; p = 0.049, p = 0.039, respectively).The majority of lakes had warming mean water column temperatures, 60% and 55% for 1975 and 1985 cohorts.

Comparison of Air and Lake Surface Temperature Trends
Near-surface temperatures of lakes in this study warmed significantly faster than corresponding mean annual air temperatures for both cohorts (two-sample t-tests: t = −4.2 and −3.7, df = 85 and 250, respectively, both p < 0.001; Figure 3).Relative to annual air temperature trends, near-surface waters were warming an average of 1.7-times faster than annual air temperature.However, on an individual lake basis, the range was highly variable (−5.4-8.4 times faster, Figure 3).

Comparison of Air and Lake Surface Temperature Trends
Near-surface temperatures of lakes in this study warmed significantly faster than corresponding mean annual air temperatures for both cohorts (two-sample t-tests: t = −4.2 and −3.7, df = 85 and 250, respectively, both p < 0.001; Figure 3).Relative to annual air temperature trends, near-surface waters were warming an average of 1.7-times faster than annual air temperature.However, on an individual lake basis, the range was highly variable (−5.4-8.4 times faster, Figure 3).In both cohorts, the near-surface temperature trends are significantly greater than the corresponding regional air temperature trends (two-sample t-tests: p < 0.001 and p < 0.001, respectively).In boxplots, boxes span the first through third quartiles with the median trend represented by the thick center line with whiskers adding the ±1.5 interquartile range to the boxes.Open circles are individual points that fall outside that span.

Median Buoyancy Frequency
Lakes identified as dimictic by data providers had significantly higher MBF than those identified as polymictic (two-sample t-test, p < 0.001; Figure 4a), confirming the expectation that MBF is a good index of mixing regime.Polymictic lakes (lower MBF) had greater mean lake temperature increases than dimictic lakes (higher MBF) for both the 1975 (r = 0.37, p < 0.001; Figure 4b) and 1985 cohorts (r = 0.20, p = 0.003; Figure 4c).Box plots showing the trends in lake near-surface temperature for NENA lakes and in regional air temperature for both the 1975 cohort (left panel) and 1985 cohort (right panel).In both cohorts, the near-surface temperature trends are significantly greater than the corresponding regional air temperature trends (two-sample t-tests: p < 0.001 and p < 0.001, respectively).In boxplots, boxes span the first through third quartiles with the median trend represented by the thick center line with whiskers adding the ±1.5 interquartile range to the boxes.Open circles are individual points that fall outside that span.

Median Buoyancy Frequency
Lakes identified as dimictic by data providers had significantly higher MBF than those identified as polymictic (two-sample t-test, p < 0.001; Figure 4a), confirming the expectation that MBF is a good index of mixing regime.Polymictic lakes (lower MBF) had greater mean lake temperature increases than dimictic lakes (higher MBF) for both the 1975 (r = 0.37, p < 0.001; Figure 4b) and 1985 cohorts (r = 0.20, p = 0.003; Figure 4c).

Boosted Regression Tree Models
Observed trends in thermal metrics were positively correlated with modeled trends using BRT (Figure S4).Lake cohort (1975 or 1985) explained very little of the variation in thermal metric trends (<3%) across all models.BRT models performed best when explaining variability in density trends (r = 0.67), deepwater temperature trends (r = 0.66) and temperature difference trends (r = 0.61; Figure 5, Figure S4).Across all models, latitude, depth index and Secchi depth had the largest relative influence on long-term thermal metric trends; however, the most important explanatory variable in each model differed depending on the thermal metric trends being modeled (Figure 5).

Boosted Regression Tree Models
Observed trends in thermal metrics were positively correlated with modeled trends using BRT (Figure S4).Lake cohort (1975 or 1985) explained very little of the variation in thermal metric trends (<3%) across all models.BRT models performed best when explaining variability in density trends (r = 0.67), deepwater temperature trends (r = 0.66) and temperature difference trends (r = 0.61; Figure 5, Figure S4).Across all models, latitude, depth index and Secchi depth had the largest relative influence on long-term thermal metric trends; however, the most important explanatory variable in each model differed depending on the thermal metric trends being modeled (Figure 5).Relative influence of each of the explanatory variables on the trends in slopes of temperature and stratification for NENA lakes from boosted regression tree analysis.For the response slopes: Surf. is the near-surface temperature, Deep. is deepwater temperature, Mean is the depth-weighted mean temperature, Diff. is the temperature difference between near-surface and deepwater, BF is the max.buoyancy frequency and Density is the density gradient.The explanatory variables are indicated by the colors in the legend on the right and are latitude (Lat.), elevation (Elev.),distance to coast (Dist.),fetch (Fetch), depth index (Depth) and Secchi depth (Secc).r and RMSE represent the Pearson correlation coefficients and root mean square error, respectively, between the modeled and actual Sen's slopes.The bars may not add to 100%; the cohort category accounted for <3% of the relative influence with each response variable.
Secchi depth was the most important explanatory variable for near-surface temperature trends and temperature difference trends (Figures 5 and 6a,b).We identified positive non-linear relationships between water clarity and these two metrics, with a step change between 5-and 7-m Secchi depth (Figure 6a).Geographic variables were the most important explanatory variables for some metrics: distance to coast for deepwater temperature trends (Figure 6c) and latitude for both mean temperature (Figure 6d) and max.buoyancy frequency trends (Figure 6e).Lakes nearer to the coast (distances <60 km) had cooling deepwater temperatures, while lakes farthest from the coast (distances >100 km) had weak and variable trends (both warming and cooling) in deepwater temperature (Figure 6c).However, lakes ranging from ~60 to 100 km from the eastern seaboard generally had warming deepwater (Figure 6c).Similarly, elevation was also an important explanatory variable for deepwater temperature trends (Figure 5) with maximum deepwater warming occurring slightly above median elevations for the lakes in this study (~150 m asl).Higher latitude lakes had warming mean temperatures and increasing strength of stratification; lower latitude lakes had little change or cooling mean temperatures and more slowly increasing strength of stratification (Figure 6d,e).The depth index was the most important explanatory variable for density gradient trends (Figure 6f); lakes with depths relative to transparency between 2 and 4 m 0.5 had the highest increasing trends in density gradient (Figure 6f).Relative influence of each of the explanatory variables on the trends in slopes of temperature and stratification for NENA lakes from boosted regression tree analysis.For the response slopes: Surf. is the near-surface temperature, Deep. is deepwater temperature, Mean is the depth-weighted mean temperature, Diff. is the temperature difference between near-surface and deepwater, BF is the max.buoyancy frequency and Density is the density gradient.The explanatory variables are indicated by the colors in the legend on the right and are latitude (Lat.), elevation (Elev.),distance to coast (Dist.),fetch (Fetch), depth index (Depth) and Secchi depth (Secc).r and RMSE represent the Pearson correlation coefficients and root mean square error, respectively, between the modeled and actual Sen's slopes.The bars may not add to 100%; the cohort category accounted for <3% of the relative influence with each response variable.
Secchi depth was the most important explanatory variable for near-surface temperature trends and temperature difference trends (Figures 5 and 6a,b).We identified positive non-linear relationships between water clarity and these two metrics, with a step change between 5-and 7-m Secchi depth (Figure 6a).Geographic variables were the most important explanatory variables for some metrics: distance to coast for deepwater temperature trends (Figure 6c) and latitude for both mean temperature (Figure 6d) and max.buoyancy frequency trends (Figure 6e).Lakes nearer to the coast (distances <60 km) had cooling deepwater temperatures, while lakes farthest from the coast (distances >100 km) had weak and variable trends (both warming and cooling) in deepwater temperature (Figure 6c).However, lakes ranging from ~60 to 100 km from the eastern seaboard generally had warming deepwater (Figure 6c).Similarly, elevation was also an important explanatory variable for deepwater temperature trends (Figure 5) with maximum deepwater warming occurring slightly above median elevations for the lakes in this study (~150 m asl).Higher latitude lakes had warming mean temperatures and increasing strength of stratification; lower latitude lakes had little change or cooling mean temperatures and more slowly increasing strength of stratification (Figure 6d,e).The depth index was the most important explanatory variable for density gradient trends (Figure 6f); lakes with depths relative to transparency between 2 and 4 m 0.5 had the highest increasing trends in density gradient (Figure 6f).(e) max.buoyancy frequency; and (f) density gradient and the explanatory variable with the highest relative influence for each thermal metric for NENA lakes.The lines represent the individual fits for the explanatory variable while holding all other explanatory variables at their median value.Note, the y-axis scales are modified for ease of presentation.
Of the 19 studies that report global, regional or single-lake LSWT trends, none reported average LSWT cooling over a similar or longer time period than this study (Table 2).Of the 13 of these studies that also reported corresponding air temperature trends, 85% indicated that LSWT was warming more rapidly than air temperature over the same time period (Table 2).Torbick et al. [9] covered a similar geographic region using satellite data from 1984 to 2012 and found faster warming rates for both LSWT and regional air temperature than we report in this study, but similarly showed LSWT was warming more rapidly than air temperature in this region (Table 2).The lines represent the individual fits for the explanatory variable while holding all other explanatory variables at their median value.Note, the y-axis scales are modified for ease of presentation.

Meta-Analysis
In comparison to three prominent studies of globally-distributed LSWT, our average near-surface trends warmed much more rapidly than the average reported in Kraemer et al.
Of the 19 studies that report global, regional or single-lake LSWT trends, none reported average LSWT cooling over a similar or longer time period than this study (Table 2).Of the 13 of these studies that also reported corresponding air temperature trends, 85% indicated that LSWT was warming more rapidly than air temperature over the same time period (Table 2).Torbick et al. [9] covered a similar geographic region using satellite data from 1984 to 2012 and found faster warming rates for both LSWT and regional air temperature than we report in this study, but similarly showed LSWT was warming more rapidly than air temperature in this region (Table 2).
Table 2. Lake surface water temperature (LSWT) trends and corresponding air temperature (AT) in this study of northeastern North American lakes (top row) and 19 other studies of lake surface temperature trends ranging from a global suite of lakes to a single lake.Presented trends are mean • C decade −1 ± SE unless otherwise noted."Type" indicates how measurements were taken ("IS" = in situ, "SAT" = satellite).The number following in parentheses indicates the number of lakes included in the study."Location" indicates if the study was global, regional ("R") or single lake ("S"), with the 2-character country code and 3-character state or province name (regional studies) or the lake name (single-lake studies).An asterisk in some single-lake studies indicates multiple sites in the lake.See footnotes for additional details.

Discussion
NENA lakes had near-surface water trends that were consistently warming at rates that were higher than air temperature warming rates.This suggests that the region may be a warming hotspot for lakes, similar to other identified hotspot areas such as around the Laurentian Great Lakes and northern Europe [12].Surface trends were also accompanied by increased strength of lake stratification.Deepwater temperatures, however, did not warm in concert with their surface waters and showed a range of changes from warming to cooling.We found that clear lakes are warming most rapidly; however, lake depth and size showed more variable trends.Further, lakes 60-100 km inland had the greatest deepwater warming indicating the effect of the ocean and geography on moderating local climate and lake mixing.

Near-Surface Warming
The vast majority of our study lakes (>80% for both cohorts) exhibited surface water warming trends, consistent with single-lake, regional and global-scale studies reporting strong long-term surface temperature warming (Table 2; [46]).Warming surface waters during peak summer stratification in the NENA region are consistent with, and likely related to, increasing air temperatures in the region [3].Yet, in this and in past studies in other regions (Table 2), surface waters are warming at a faster rate than annual regional air temperatures [14].In the most extreme results we reviewed, six large lakes in California and Nevada showed surface warming at twice the rate of regional air temperature changes based on satellite imagery (1991-2008; [34]) relative to 50% faster in this study (Table 2).Air temperature trends alone underestimate LSWT in models of lake surface equilibrium temperature (e.g., [19,47] predict LSWT at 70-80% of air temperature warming rates), suggesting that such models may be missing key integrative factors.
Intensified surface water warming trends that exceed rising air temperature rates may be related to a variety of additional factors that affect light energy arriving at the lake surface, light energy that has entered the lake and energy radiated from the lake.Solar brightening [39,48] and decreases in cloud cover [49] lead to increased solar radiation reaching the lake surface.Schmid and Köster [39] found that 60% of LSWT in a central European lake was attributable to increasing air temperature, while the remaining 40% was attributable to increased solar radiation due to air quality improvements.In the future, LSWT warming rates may be more consistent with air temperature rates depending on future stabilization in air quality improvements [4,39].Similarly, milder winters, shortening ice cover and less snow cover on ice occurring in NENA and around the northern latitudes [50,51] will all contribute to lowered lake albedo, earlier stratification and longer duration of energy absorption, resulting in intensified summer warming [11,12,52].Changes in wind speed [53,54], water transparency [55,56] and hydrologic inflows [54] can all alter lake warming rates.
A small proportion of lakes had cooling near-surface temperatures (Figure 2, Tables S2 and S3).Increases in tree cover or shading at the edge of the lake or other land-use changes within the watershed might decrease solar radiation entering the lake [12].Water transparency can control near-surface temperature trends through controlling vertical light absorption at the surface [56].For example, less transparent lakes could have increased rates of surface absorption and outgoing radiation, especially at night [57], which would lead to decreasing surface warming trends, a rare occurrence in NENA lakes.Additionally, sunlight can only warm a fraction of the mixed layer in less transparent lakes [47], which limits light and heat penetration to sub-surface layers.In NENA in particular, additional regional drivers of increasing transparency (acidification and invasive mussels) or decreasing transparency (browning and land use change) could modify patterns in LSWT (see below).

Deepwater Trends
Across all NENA lakes, there was no significant overall change in deepwater temperatures; however, there was high variability among lakes with divergent patterns of warming and cooling.
While deepwater temperatures were related to shifts in near-surface waters (Figure S3), the explanatory power was low, suggesting additional controls on deepwater temperatures, with a general decoupling of air temperature trends from deepwater temperature trends.Similar to our study, a suite of lakes in Wisconsin had deepwater trends with high variability around a mean of zero; however, the range of deepwater trends in our study was larger [18].
In dimictic NENA lakes, deepwater temperature trends follow geographic patterns that likely reflect a number of variables, including underlying geology, mixing regime, conditions during spring mixing and precipitation.Distance to the Atlantic coast was important to deepwater trends with peak warming at 60-100 km inland (Figure 6c) in the Piedmont of the Appalachian Mountains (Figure S5).Similarly, elevation trends indicate a comparable spatial pattern with maximum deepwater trends along the eastern Piedmont.As we move inland from the coast, moving from the west to the eastern edge of the Appalachians, a transition from coastal plain to piedmont and highland physiography is apparent, and this transition is associated with changing ecoregions, increasing elevation, more topographic shading, due to the surrounding terrain, and more temperature variability [58].In this transition region, we expect maximum change in ice cover duration, mixing regime shifts and wind-driven turbulent mixing, especially for large lakes (e.g., [18]).In deep lakes that stably stratify, the deepwater trends are more reflective of late winter and early spring conditions during turnover and onset of stratification, including wind-driven mixing following ice-off, as opposed to annual or summer air temperatures [59,60].Underlying geology can affect deepwater temperature through modifications of groundwater fluxes [18] and sediment heat storage in the shallower zones of lakes [61].In addition to warming air temperatures, the NENA region has experienced the greatest increase in precipitation and extreme storm events in the USA [3].East of the Appalachians receives more precipitation; changes in the delivery of this surface or groundwater flow may be critical factors influencing deepwater temperature in this region, as has been observed in a deep, tropical lake [62].In small, sheltered, but stratified lakes where wind mixing is minimal, water transparency may be a key modifier of deepwater trends [22,57,60].Low water transparency limits heat penetration and restricts transfer of thermal energy to deeper waters [56]; however, in this study, transparency was not a clear driver of deepwater trends as a whole (Figure 5).Deepwater temperatures respond to different drivers or mediators of lake thermal structure than near-surface temperature and are generally disconnected from air temperature trends.

Thermal Stratification Trends
The strength of thermal stratification is increasing in lakes in this region and around the world, largely due to the greater warming of surface temperatures relative to the deepwater temperatures, resulting in steeper thermoclines and more stable stratification ([16]; Figure 2).In this study, all three stratification metrics (temperature difference, density difference and max.buoyancy frequency) indicate increasing strength of stratification over time in the majority of NENA lakes (Figure 2).These stratification trends reflect greater increases in surface water temperature relative to deeper waters.Additionally, warmer late winter and early spring air temperatures combined with earlier ice off and earlier stratification will increase thermocline stability over the course of the season, resulting in stronger peak stratification [19].The strength of thermal stratification is increasing more rapidly in lakes of intermediate depth (Figure 6f) with larger heat budgets and more resistance to increases in total thermal content throughout the water column than smaller lakes [19].While the smallest lakes are often polymictic, smaller dimictic lakes can be sheltered from wind and turbulent mixing, with minimal transference of the climate signal below the thermocline, resulting in increased stratification strength [18,59,63].Similarly, reduced mean annual wind speed in recent decades in this region [64] may lead to enhanced thermal stratification.In lakes with decreasing transparency, light attenuation, absorption and reflection and re-radiation of heat energy will all occur nearer the surface, preventing heat from reaching deeper waters and strengthening thermal stratification [57].
Latitude was an important explanatory variable for the strength of thermal stratification and mean lake temperature (Figure 6d,e and Figure S6).Lakes at higher latitudes (above ~44 • in this study) experienced the most rapid increases in the strength of thermal stratification and the most rapid mean lake warming rates.Although the span of latitude in this study is somewhat limited, higher latitudes tend to experience the greatest air temperature warming rates [3,65], which translate to warmer lake surface waters and greater density gradients between lake strata.Further, more rapidly warming air temperatures at higher latitudes may lead to changes in ice cover phenology, i.e., earlier ice break-up dates in particular, and through extended periods of stratification, rapid northern warming may increase the strength of stratification by the peak summer period.
Mixing regime explained mean temperature trends (Figure 4) because turbulent mixing can control thermal transfer from the air-water interface to deeper waters.Polymictic lakes have regular mixing events, which can distribute the heat throughout the water column and tightly couple air and water temperature throughout the water column, including deepwater temperature [66].In polymictic lakes, weak stratification resulted in more uniform thermal distribution and strong mean lake warming (Figure 4).In dimictic lakes, surface waters warm with resistance to heat transfer below the thermocline, which leads to increases in the strength of thermal stratification (Figure 2), yet results in weaker mean lake response (Figure 4).Dimictic lakes experiencing both near surface warming and deepwater cooling, with stable stratification, showed the strongest increases in the strength of thermal stratification.

Changes in Transparency in the NENA Region
Preceding and during the time period of our study , several lakes in the NENA region experienced a number of anthropogenic stressors that may either decrease or increase water transparency.In many NENA lakes, acidification via acid precipitation was common [67,68].Up to the early 1990s, this resulted in increasing water clarity due to the reduction of sensitive phytoplankton taxa [69] and lower concentrations of dissolved organic carbon [70].Following the Clean Air Act Amendment in 1990, acid deposition to lakes and watersheds decreased and led to increases in dissolved organic carbon (lake browning) and resulted in decreased water clarity, trends that are prevalent throughout the NENA region [71][72][73].Long-term browning trends are likely influenced by both recovery from acidification [67,68,74] and increased precipitation and storm events occurring in the NENA region [75,76].For example, Lakes Giles and Lacawac in northeastern Pennsylvania had steadily increasing dissolved organic carbon concentrations over 27 years.Giles had strong reductions in water transparency, which resulted in strong increases in surface temperature and the strength of thermal stratification, but this lake had significantly cooling deepwater temperatures (this study; [73]).In the NENA region, land-cover from 1973 to 2000 shifted with increased abandonment of agricultural areas, increased areas of developed land (e.g., urban, suburban, industrial, golf courses) and reduced forested areas [77].This may have led to high lake nutrient concentrations and reduced water clarity as found in many NENA lakes [78].On the other hand, invasive species, such as non-native zebra and quagga mussels (Dreissena spp.), invaded a few of the NENA lakes by the mid-1980s and can dramatically increase water clarity due to their efficient filtering of phytoplankton (review in [79]).For example, within two years of the discovery of the Dreissena spp.invasion, Oneida Lake in central New York experienced a dramatic and sustained 40% increase in water clarity [80].Not all lakes in the NENA region were subjected to all of these external stressors, and corresponding trends in transparency may be non-linear (e.g., Dreissena spp.invasion); however, these stressors may explain some of the remaining variability that we could not explicitly address in our study.We plan to examine these other drivers of transparency change in the subset of lakes where temporally explicit data on transparency trends exist, but are not currently compiled.

Ecological Implications
Extensive modifications to lake physical properties and heat distribution are likely causing a variety of ecosystem-level consequences with concomitant impacts on ecosystem services.Lake warming is often coupled with changes in ice phenology, where many lakes are experiencing a shortened ice cover duration in temperate regions [51,81,82].This may ultimately provide a positive feedback for lake warming [11] and result in a longer stratified period with shorter ice cover periods.A longer stratification period each year will exacerbate the extent and duration of hypolimnetic hypoxic or anoxic conditions [21,36,83].In many lakes, changes in the amount of vertical mixing can also contribute to changes in nutrient dynamics and ultimately lake productivity [84,85] and changing fisheries [86].Warming and increasingly anoxic deepwater can both result in increased cyanobacterial recruitment in species that both overwinter in and draw phosphorus from the sediments (e.g., [87]).Increasingly warm surface temperatures could favor harmful cyanobacterial algal blooms and transition lakes from less productive towards eutrophic [88,89].Oxygen, lake temperature and water transparency also regulate predator-prey interactions through the vertical distribution and migration of phytoplankton and zooplankton and thereby modify the availability of refugia from predation (e.g., [90][91][92]).Further, a compressed or non-existent hypolimnion, due to a deeper thermocline and increasingly hypoxic conditions, will have costly energetic consequences for cooland cold-water fishes [93].However, the opposite pattern could emerge in lakes with increasing hypolimnetic volume, assuming ample deepwater oxygen concentrations.In NENA lakes, fisheries are frequently supported by a mix of warm-, cool-and cold-water species.Long-term warming trends will favor warm water fish species (e.g., black bass) as has already been observed in NENA lakes, resulting in the dominance of warm-water species [94,95].We expect that the thermal changes documented in the NENA region are having and will continue to have profound ecological consequences for these systems.

Conclusions
Lakes are critically important sentinels of climate change, but lake characteristics within and across regions (e.g., transparency, mixing regime, geomorphology, geography and trophic status) moderate the magnitude of the signal from climate change.We identified not only a consistent climate signal from warming air temperatures in this region, with near-surface water warming and increasing strength of thermal stratification, but we also found variability in the signal for subsurface thermal stratification trends.Studies focusing on individual large lakes (e.g., Baikal, Superior) are critically important and will document changes occurring in the majority of freshwater by volume on Earth.However, most lakes are small relative to these large, well-studied lakes [1,96] and may respond differently to climate change (e.g., this study; [16,18,36]).Given the likely trajectory of climate warming and changing precipitation patterns [3,65], there is a continual need to track lake thermal structure as it interacts with other external forcings on lake ecosystems and watersheds and to consider the implications for lake biology, ecology, chemistry and ecosystem services.

Figure. 1 .
Figure. 1. Location of study lakes for the northeastern North America (NENA) study.The 1975 cohort of lakes (centroids) is depicted as dark blue circles (n = 85), and the 1985 cohort of lakes is depicted as light blue circles (n = 226).Note that the 1975 cohort of lakes includes the 1985 cohort of lakes (except for five lakes).

Figure 1 .
Figure 1.Location of study lakes for the northeastern North America (NENA) study.The 1975 cohort of lakes (centroids) is depicted as dark blue circles (n = 85), and the 1985 cohort of lakes is depicted as light blue circles (n = 226).Note that the 1975 cohort of lakes includes the 1985 cohort of lakes (except for five lakes).

Figure 2 .
Figure 2. Bean plots showing the magnitude of trends for each metric of lake thermal structure for NENA lakes as follows: Near-surface temperature (Near-surf.temp), Deepwater temperature (Deepwater temp.), mean temperature (Mean temp.),Temperature difference (Temp.diff.), maximum buoyancy frequency (Buoy.freq.), and density gradient (Dens.grad.).For each metric, the left half (light blue) of the bean plot represents the distribution of the 1975 cohort (1975-2014), and the right half (dark blue) represents the distribution of the 1985 cohort (1985-2014).The four metrics of temperature follow the left y-axis (°C y −1 ); the buoy.freq.follows the first y-axis on the right (cph y −1 : cycles per hour y −1 ); and the dens.grad.follows the second y-axis on the right (kg m −3 y −1 ).Asterisks indicate a significant difference from zero.

Figure 2 .
Figure 2. Bean plots showing the magnitude of trends for each metric of lake thermal structure for NENA lakes as follows: Near-surface temperature (Near-surf.temp), Deepwater temperature (Deepwater temp.), mean temperature (Mean temp.),Temperature difference (Temp.diff.), maximum buoyancy frequency (Buoy.freq.), and density gradient (Dens.grad.).For each metric, the left half (light blue) of the bean plot represents the distribution of the 1975 cohort (1975-2014), and the right half (dark blue) represents the distribution of the 1985 cohort (1985-2014).The four metrics of temperature follow the left y-axis ( • C y −1 ); the buoy.freq.follows the first y-axis on the right (cph y −1 : cycles per hour y −1 ); and the dens.grad.follows the second y-axis on the right (kg m −3 y −1 ).Asterisks indicate a significant difference from zero.

Figure 3 .
Figure 3. Box plots showing the trends in lake near-surface temperature for NENA lakes and in regional air temperature for both the 1975 cohort (left panel) and 1985 cohort (right panel).In both cohorts, the near-surface temperature trends are significantly greater than the corresponding regional air temperature trends (two-sample t-tests: p < 0.001 and p < 0.001, respectively).In boxplots, boxes span the first through third quartiles with the median trend represented by the thick center line with whiskers adding the ±1.5 interquartile range to the boxes.Open circles are individual points that fall outside that span.

Figure 3 .
Figure 3. Box plots showing the trends in lake near-surface temperature for NENA lakes and in regional air temperature for both the 1975 cohort (left panel) and 1985 cohort (right panel).In both cohorts, the near-surface temperature trends are significantly greater than the corresponding regional air temperature trends (two-sample t-tests: p < 0.001 and p < 0.001, respectively).In boxplots, boxes span the first through third quartiles with the median trend represented by the thick center line with whiskers adding the ±1.5 interquartile range to the boxes.Open circles are individual points that fall outside that span.

Figure 4 .
Figure 4. (a) The median buoyancy frequency (MBF) from the entire multi-year record of NENA lakes for each lake with the clear identification of polymictic or dimictic mixing regimes.Correlation of MBF with the mean temperature slope for (b) 1975 and (c) 1985 cohorts.The LOESS line has been added to visualize the relationship.In the boxplots, boxes span the first through third quartiles with the median trend represented by the thick center line with whiskers adding the ±1.5 interquartile range to the boxes.Open circles are individual points that fall outside that span.

Figure 4 .
Figure 4. (a) The median buoyancy frequency (MBF) from the entire multi-year record of NENA lakes for each lake with the clear identification of polymictic or dimictic mixing regimes.Correlation of MBF with the mean temperature slope for (b) 1975 and (c) 1985 cohorts.The LOESS line has been added to visualize the relationship.In the boxplots, boxes span the first through third quartiles with the median trend represented by the thick center line with whiskers adding the ±1.5 interquartile range to the boxes.Open circles are individual points that fall outside that span.

Figure 5 .
Figure 5.Relative influence of each of the explanatory variables on the trends in slopes of temperature and stratification for NENA lakes from boosted regression tree analysis.For the response slopes: Surf. is the near-surface temperature, Deep. is deepwater temperature, Mean is the depth-weighted mean temperature, Diff. is the temperature difference between near-surface and deepwater, BF is the max.buoyancy frequency and Density is the density gradient.The explanatory variables are indicated by the colors in the legend on the right and are latitude (Lat.), elevation (Elev.),distance to coast (Dist.),fetch (Fetch), depth index (Depth) and Secchi depth (Secc).r and RMSE represent the Pearson correlation coefficients and root mean square error, respectively, between the modeled and actual Sen's slopes.The bars may not add to 100%; the cohort category accounted for <3% of the relative influence with each response variable.

Figure 5 .
Figure 5.Relative influence of each of the explanatory variables on the trends in slopes of temperature and stratification for NENA lakes from boosted regression tree analysis.For the response slopes: Surf. is the near-surface temperature, Deep. is deepwater temperature, Mean is the depth-weighted mean temperature, Diff. is the temperature difference between near-surface and deepwater, BF is the max.buoyancy frequency and Density is the density gradient.The explanatory variables are indicated by the colors in the legend on the right and are latitude (Lat.), elevation (Elev.),distance to coast (Dist.),fetch (Fetch), depth index (Depth) and Secchi depth (Secc).r and RMSE represent the Pearson correlation coefficients and root mean square error, respectively, between the modeled and actual Sen's slopes.The bars may not add to 100%; the cohort category accounted for <3% of the relative influence with each response variable.

Figure 6 .
Figure 6.Relationship between boosted regression tree model fits (yellow points) for trends of: (a) the near-surface temperature; (b) temperature difference between near-surface and deepwater; (c) deepwater); (d) mean temperature;(e) max.buoyancy frequency; and (f) density gradient and the explanatory variable with the highest relative influence for each thermal metric for NENA lakes.The lines represent the individual fits for the explanatory variable while holding all other explanatory variables at their median value.Note, the y-axis scales are modified for ease of presentation.

Figure 6 .
Figure 6.Relationship between boosted regression tree model fits (yellow points) for trends of: (a) the near-surface temperature; (b) temperature difference between near-surface and deepwater; (c) deepwater); (d) mean temperature; (e) max.buoyancy frequency; and (f) density gradient and the explanatory variable with the highest relative influence for each thermal metric for NENA lakes.The lines represent the individual fits for the explanatory variable while holding all other explanatory variables at their median value.Note, the y-axis scales are modified for ease of presentation.

Table 1 .
Summary of lake variables from this northeastern North American lake study, 1985 cohort.The ranges for the 1975 cohort are consistent with those listed here.See the Materials and Methods for further details."BRT" indicates if the explanatory variables were included (Y) or not (N) in the boosted regression tree analysis.
Note: * Included as a transparency-independent depth metric; see the Materials and Methods.
since the 19th century; l April, May and June trend; m IS; whole-lake average lake surface water warming trends (LSWT); n IS; whole-lake average LSWT; o SAT; whole-lake average LSWT.