Deeper waters are changing less consistently than surface waters in a global analysis of 102 lakes

Globally, lake surface water temperatures have warmed rapidly relative to air temperatures, but changes in deepwater temperatures and vertical thermal structure are still largely unknown. We have compiled the most comprehensive data set to date of long-term (1970–2009) summertime vertical temperature profiles in lakes across the world to examine trends and drivers of whole-lake vertical thermal structure. We found significant increases in surface water temperatures across lakes at an average rate of + 0.37 °C decade−1, comparable to changes reported previously for other lakes, and similarly consistent trends of increasing water column stability (+ 0.08 kg m−3 decade−1). In contrast, however, deepwater temperature trends showed little change on average (+ 0.06 °C decade−1), but had high variability across lakes, with trends in individual lakes ranging from − 0.68 °C decade−1 to + 0.65 °C decade−1. The variability in deepwater temperature trends was not explained by trends in either surface water temperatures or thermal stability within lakes, and only 8.4% was explained by lake thermal region or local lake characteristics in a random forest analysis. These findings suggest that external drivers beyond our tested lake characteristics are important in explaining long-term trends in thermal structure, such as local to regional climate patterns or additional external anthropogenic influences.


A) Lake Information and Metadata
Table S1: Detailed information and metadata for lakes included in the analysis. Lake data include lake name, location, time period(s) included in analysis, mean depth increment of sampling (m), mean number of profiles per year, month of approximate peak thermal stability, thermal region classification, latitude (°), longitude (°), elevation (m above sea level), surface area (km 2 ), maximum depth (m), average Secchi depth (m), average chlorophyll-a (μg L -1 ), and average dissolved organic carbon (DOC; mg L -1 ). Chemical and limnological measurements represent averages from recent years, and do not capture any long-term changes in the variable.
Data can be viewed in the provided Supplementary Dataset online.

B) Single Profile Method Comparison
Due to time and resource demands, many lakes are not sampled frequently throughout the year, and often not at all during the ice-free season. Hence, using one summer profile to represent stable summertime stratification can maximize utility of data from infrequently samples lakes in a way that is useful for long-term trend analysis. This methodology enabled us to utilize a much larger number of lakes, which are generally sampled twice a month or less, and allows for analyses that account for broader range of geographic, morphological, and chemical variability in lakes. From this dataset, nearly half of the total of 102 lakes would have been excluded from the analysis if more frequent summer sampling was required.
A single summer profile is an important sentinel for capturing lake thermal structure and, in general, provides a similar signal of change compared to metrics averaged across one or several months during summertime. The concept is that including these other periods of lower thermal stability will likely result in a weaker signal due to a smaller signal to noise ratio. We used a sub-sample of the 102 lakes with frequent summer sampling to analyse if a single summer profile from within the period of maximum thermal stratification adequately captured general patterns in summer thermal structure.
To assess this methodology, a sub-sample of lakes with a minimum of 15 years of data including at least one profile from each of three summer months were used in this analysis. A sub-sample of n = 51 lakes all located in the Northern Hemisphere met these criteria. We compared the long-term trends from the single profile selection method with those from three individual summer months (June, July, and August), and with the aggregated summer month trend (June, July, and August = "JJA" average). We used the same five thermal response metrics included in the main analysis (surface water temperature, deepwater temperature, mean water column temperature, density difference, and thermocline depth).
Thermal metrics were calculated from all available summer profiles and were then aggregated by individual month or JJA summer average. For each metric for each lake, we calculated the trends based on individual month data, the JJA summer average data, and the single profile data over time from 1990-2009 using Sen's slope. We used a 1-way ANOVA blocked by lake for each thermal metric to assess the difference in Sen's slope between the time frames, with a significance level of α = 0.05. If the ANOVA was statistically significant, we followed with a Tukey's honest significant difference test to determine which time frames showed significantly different slopes across lakes.
Our results show that there is no detectable difference in the overall trends when using different time frames for four of the five thermal metrics (Fig. S1). Only density difference showed a significant difference in trends from these different time frames (Fig. S1d, p = 0.01).
For this metric, trends based on July data (A) were significantly greater compared to those based on August data (B). Trends from June (AB), summertime average (AB), and the single profile method (AB) were statistically indistinguishable from each other. This may represent the integrated signal of stability that is often greatest in July in many Northern Hemisphere lakes, where changes in phenology of this metric may also be important to consider. For the metrics considered here, the single profile method is highly comparable to other summertime aggregation methods for assessing overall long-term trends in thermal structure. Figure S1. Comparison of Sen's slope estimates when using data from the single profile selection method, three individual summer months, and the JJA summer average for 51 Northern Hemisphere lakes. Only density difference trends (d) had significant differences across the time periods when blocked by lake (p = 0.01). Letters at the bottom of panel (d) indicate Tukey's honest significant differences where only periods without the same letter differ significantly from each other. There were no significant differences across time periods for the other four metrics: (a) surface water temperature trends (p = 0.17); b) deepwater temperature trends (p = 0.14); (c) mean water column temperature trends (p = 0.43); (e) thermocline depth trends (p = 0.34). All statistical analyses are blocked by lake, so differences across lakes (which is generally expected) are not factored into the statistical output presented here.

C) Strength of Thermal Stratification Metric Comparisons
It is important to analyse complex data such as these thoroughly, especially as different approaches or metrics can result in different temporal trends [1] . There are many possible ways to measure thermal structure and stability given the high volume of vertical temperature profile data used for this analysis. For example, Kraemer et al. [1] showed that metrics of stability and stratification do not always indicate the same overall patterns, even when calculated from the same set of lakes. To assess potential differences in metrics of thermal stratification, we explicitly compared six metrics of strength of thermal stratification that can potentially capture different aspects of change over time. We compared these different metrics of strength of thermal stratification in this set of lakes on both interannual and decadal scales, and addressed the differences in the measurements in order to select the most useful metric of strength of thermal stratification for our main analysis.
We calculated six metrics of strength of thermal stratification for each single summer profile for each lake. The temperature readings were interpolated or binned to 0.5 m increments, as described in the methods of the main text. The span of data was each lake's complete individual record since these intercomparisons are all within-lake, but required a minimum of 15 years of data, leading to a total of n = 101 lakes included in this comparison. The six metrics of thermal stratification under consideration were: 1) Temperature difference ("Temp. Diff."; °C): the surface water temperature reading at 2 m minus the deepwater temperature taken at the deepest consistently-recorded depth (varies by lake).
3) Relative thermal resistance to mixing ("RTR", unitless): the ratio of density difference between the deepwater minus the surface water temperature divided by the density difference between 4°C and 5°C [2,3] . See equation (1) in main text. 4) Schmidt stability* ("SS*", units = J m -2 ): the energy required to fully mix the lake vertically. Note that the lack of the inclusion of bathymetric data means that SS* represents the amount of work necessary to bring the water column to uniform density; calculated using the R package "rLakeAnalyzer" [4] . 5) Summed buoyancy frequency ("Sum BF"; cycles hour -1 ): the summed Brunt-Väisälä frequency from each 0.5 m depth interval in the full water column; calculated using the R package "rLakeAnalyzer" [4] . 6) Maximum buoyancy frequency ("Max BF"; cycles hour -1 ): the maximum Brunt-Väisälä frequency at the seasonal thermocline; calculated using the R package "rLakeAnalyzer" [4] .
The time series for each metric was individually standardized within-lake using z-scores.
We used non-parametric Kendall correlations to assess the similarity of all the different metrics for all years within each lake over the full available data record (n = 101 correlations for each metric pairing). To assess the correlation between long-term trends in these different metrics, we calculated Sen's slope from the above metrics spanning the full available data record for all lakes. We computed non-parametric Kendall correlations with the resulting trend estimates across all combinations of the six metrics, where each lake served as a replicate.
The correlations indicated that the magnitude and direction of differences between years within lakes and the long-term trends across lakes were broadly similar between metrics of thermal stratification. Overall, the interannual variability within lakes in the various metrics of the strength of thermal stratification were highly correlated (Fig. S2). Five of the metrics were highly correlated interannually, with median Kendall τ ranging from 0.75 (SS* vs. Temp. Diff.) to 1.00 (Dens. Diff. vs. RTR; Fig. S2). However, maximum buoyancy frequency was consistently less correlated with the other five metrics, with median Kendall τ ranging from 0.37 (vs. SS*) to 0.47 (vs. Temp. Diff.; Fig. S2). There was also very high variability across lakes for maximum buoyancy frequency with the other five metrics, occasionally even indicating a small, negative correlation (Fig. S2).
Similarly, the correlation of the long-term trends in the metrics across lakes was generally high (Table S2). But again, maximum buoyancy frequency had the lowest set of correlation coefficients with the other five metrics, ranging from as little as τ = 0.27 (vs. SS*) to only τ = 0.35 (vs. Temp. Diff. and vs. Dens. Diff.). All other comparisons had a correlation coefficient of τ = 0.65 or greater, indicating similar long-term signals across this set of global lakes (Table S2).  Table S2. Correlations between long-term trends of six different metrics of lake thermal stratification. Kendall correlation coefficients (τ) of the long-term trends in n = 101 lakes between each pair of metrics of thermal stratification. Note that the diagonals are identity (τ = 1).