Diel Surface Temperature Range Scales with Lake Size

Ecological and biogeochemical processes in lakes are strongly dependent upon water temperature. Long-term surface warming of many lakes is unequivocal, but little is known about the comparative magnitude of temperature variation at diel timescales, due to a lack of appropriately resolved data. Here we quantify the pattern and magnitude of diel temperature variability of surface waters using high-frequency data from 100 lakes. We show that the near-surface diel temperature range can be substantial in summer relative to long-term change and, for lakes smaller than 3 km2, increases sharply and predictably with decreasing lake area. Most small lakes included in this study experience average summer diel ranges in their near-surface temperatures of between 4 and 7°C. Large diel temperature fluctuations in the majority of lakes undoubtedly influence their structure, function and role in biogeochemical cycles, but the full implications remain largely unexplored.


Introduction
Temperature is one of the most fundamental drivers of ecosystem structure and function. It affects rates and equilibria positions of chemical reactions [1] and rates of metabolic processes, especially amongst poikilothermic aquatic organisms [2][3][4]. In lakes, temperature has a pervasive effect on a large range of physical, chemical and biological attributes and processes and influences the physical structure, rates of photosynthesis and respiration [5], biological growth rates [6], organic carbon mineralization [7], greenhouse gas emissions [8][9][10], organism size [11], the timing of phenological events [12], the likelihood of toxic cyanobacterial blooms [13], and the available habitat for fish species [14]. Quantification of surface temperature variation, and the factors that control it, are therefore of paramount importance in understanding lake behaviour and function.
Lake temperatures have been measured around the world for many years at weekly to monthly frequencies, enabling an understanding of causes and wide-ranging ecological consequences of seasonal, annual and decadal temperature changes (e.g. [15]). Knowledge of diel temperature cycles is also important for, among other things, calculating biogeochemical reaction rates and gas fluxes accurately, elucidating systematic differences between different lakes and for determining whether the predominance of day-time measurements has biased limnological understanding. Until recently, the continuous high-frequency measurements required to resolve diel cycles have not been sufficiently wide-spread to allow a systematic examination of diel temperature changes to be undertaken. However, the recent establishment of scientific networks dedicated to the collaborative analysis of high-frequency data, such as GLEON (http://www.gleon.org/) and NETLAKE (https://www.dkit.ie/netlake/), have provided an opportunity for a large-scale analysis of diel temperature variability. Thus, we have collated data from 100 lakes, that are deep enough to stratify, across four continents (S1 Table, S1 Fig), for which high-frequency temperature measurements were available, to determine the extent, and causes, of variation in diel surface temperature cycles.
Upper water temperature is controlled by heat exchange across the air-water interface, which in turn is determined predominantly by incident solar radiation, cloud cover, air temperature, relative humidity, and wind speed [16]. It is also influenced by the depth of the upper mixed layer and light attenuation in the water column (e.g. [17]). To understand the overarching controls on diel temperature cycles we have, therefore, examined the influence of four integrating variables that might be expected to have a strong effect: i) latitude, which determines variation in solar radiation flux via day-length and insolation; ii) altitude, which affects air temperature via the adiabatic lapse rate; iii) lake area, which is the primary influence on thermocline depth within a lake [18]; and iv) in-lake attenuation of photosynthetically active radiation (PAR, spectral band 400 to 700 nm), which influences the vertical distribution of incoming solar radiation within a lake and is also known to influence the vertical thermal structure [19].

Methods
Near-surface (~1 m), in situ, water temperature measurements of high-frequency (data resolution ranged from minutely to hourly, and accuracy of between 0.001°C and 0.2°C; specific details (S2 Table) of the thermistors used in each lake are available upon request) were collated from 100 temperate and boreal lakes (totalling more than 200 thousand measurements) on four continents. No specific permissions were required for any of the lakes studied in this investigation and the data used were collected previously for other studies. The field studies did not involve endangered or protected species. The lakes varied in surface area between about 2.5 x 10 3 m 2 and 1.6 x 10 8 m 2 , light attenuation of PAR (K d ) between 0.08 m -1 and 5.7 m -1 , in altitude between -211 m a.s.l. and 2464 m a.s.l., in (absolute) latitude between 32.817°and 59.846°, and in maximum depth between 2 m and 256 m. Of the 100 lakes, 74 had direct or indirect measurements of light attenuation, and 24 had both meteorological measurements (wind speed, solar radiation, air temperature and relative humidity) and temperature profiles (S1 Table). Unfortunately, the majority of the lakes included in this investigation did not have information on the temporal variations in K d thus we only had single values for each lake. This may be problematic, as K d varies through time with, among other things, algal production which can have a large influence on the thermal dynamics of lakes (e.g. [20][21][22]), but similar to other large-scale studies (e.g. [23]) we were restricted to single values in this study.
The diel temperature range (DTR) of the near surface water was calculated as the difference between the maximum and minimum daily surface temperature for each lake. For the 24 lakes with meteorological data and depth-resolved temperature measurements, an approximate theoretical DTR (ΔT 0 ) was also calculated: where Δt is the heating period for a given day (taken to be 12 h, or 43200 s), A 0 is lake surface area (m 2 ), Q zmix is the net heat flux to the upper mixed layer (J m -2 s -1 ), ρ 0 is the density of the surface water (kg m -3 ), C Pw is the specific heat of water at constant pressure (4186 J kg -1°C-1 ), and V zmix is the volume of water (m 3 ) within the upper mixed layer (z mix , m). The depth of the upper mixed layer was defined as the first depth where the temperature difference was estimated to be greater than 0.2°C relative to the temperature located at 1 m (e.g. [17]). To determine the volume of water within the upper mixed layer, bathymetric maps were required. Hypsographic curves for each lake were extracted from these maps, where available, or from GPS/depth-sounder data. For locations where neither of these datasets where available, the lakes were assumed to have a conical shape constrained by surface area and maximum depth (e.g. [23]). The amount of surface heating which influences the upper mixed layer, Q zmix , was estimated following the methods detailed in Woolway et al. [17], by using the 'Lake Heat Flux Analyzer' program [24].

Statistical methods
To investigate the controls of the DTR, we used a Generalised Additive Model (GAM) with a gamma error distribution and the logarithm link function. Specifically, the GAM was used to examine the relationship between the average summer DTR and the explanatory variables. Four predictor variables (available for 74 lakes) were used to explain the variation in DTR: lake surface area (A 0 ), K d , latitude (φ), and altitude (h). Absolute latitude was used, calculated as distance from the equator i.e. irrespective of hemisphere, and a logarithmic transformation was applied to lake surface area. The light attenuation coefficient was converted to a percent transmission per metre, I z , with the following formula: We used a GAM of the following form: where μ i is the expectation of the response Y i (μ i E(Y i )) for the ith observation, g is the logarithm link function, the inverse of which maps values from the linear predictor, η i , on to the scale of the response. η i consists of a constant term, β 0 , plus four smooth functions, f j , j = {1, 2, 3, 4}, one function per covariate considered.
To implement model selection we used a double penalty approach [25], which adds a penalty on the null space of the smoother to the usual wiggliness penalty used to select the smoothness of the f j . By penalizing the null space of the smoother in addition to the wiggliness penalty, individual smooths can be completely removed from the model, thus providing a principled means of model selection. Smoothness parameters and estimates of the model coefficients (i.e. β 0 and coefficients for the f j ) were calculated using restricted maximum likelihood [26]. Models were fitted using the R programming language [27] with the mgcv package (version 1.8-3; [28]) using the select = TRUE option for the double penalty. The method used to formulate the test statistics and p-values follows Wood [28] and is a test against the null hypothesis of zero effect for the jth smooth.
In addition to the method described above, we also followed a multi-model inference approach to determine the relative importance of each of the predictor variables in determining the DTR. The R package MuMIn [29] was used to select the best possible combination of predictor variables (= 2 N , where N is the number of predictors) contained within the GAM. Models were compared using an adjusted Akaike Information Criterion (AICc). AICc is a measure of model performance, which compares the maximum-likelihood estimate of models, adjusted for increasing complexity. The model with the lowest AICc is considered to exhibit the best performance of the set tested. All models with AICc values to within four of the model with the lowest AICc value were then selected as a 'confidence set', thus including all possible models possessing a considerable level of empirical support. The confidence set of models was used to derive relative importance values for each explanatory variable. Relative importance, which represents the probability of a variable being present in the best-performing model for a particular predictor, was calculated in MuMIn using the relative Akaike weights of models within the confidence set.
Ecological and biogeochemical consequences of the diel surface temperature range We made simple calculations to illustrate the potential biases of different diel temperature cycles on estimates of temperature-dependent ecological processes. For these calculations, we assumed a mean temperature of 20°C (i.e. the mean summer temperature for the 100 lakes included in this investigation) and computed synthetic temperature cycles with amplitudes of 0.5, 3.5, and 7.5°C (equalling a diel range of 1, 7, and 15°C, respectively). These synthetic temperature cycles were computed using simple sine waves with the appropriate amplitudes (e.g. Asin(2π(1/24)t)+T mean ; where t = 1:24, A is amplitude, and T mean is the mean temperature). It is anticipated that the diel temperature cycles of lakes will not always follow a perfect sinusoid, but these temperature cycles are used for illustration here. Using these synthetic cycles we calculated the temperature-dependent CO 2 solubility from the equation in Rebsdorf et al. [30] and an atmospheric CO 2 partial pressure of 400 ppm; and oxygen solubility was calculated from equations in Mortimer [31]. Q 10 values (i.e. the factor by which a biological rate is increased by 10°C rise in temperature) were calculated for a value of 2, typical of photosynthesis and respiration and a value of 4, typical of methane emissions from lakes. Therefore, we estimated a diel cycle in each of these temperature-dependent processes and then calculated the differences between the estimated value at 20°C (i.e. the mean temperature) and that estimated at the time of maximum temperature during the diel cycle (e.g. 23.5°C for a diel cycle with an amplitude of 3.5°C). Our objective here is to illustrate the potential bias resulting from measurements taken at different times during the day. Accounting for diel temperature variability can potentially lead to large differences in estimates of temperature-dependent ecological processes.

Results and Discussion
The average DTR varied seasonally across the 100 lakes investigated in this study, being largest in summer and lowest in winter (Fig 1A). We subsequently restricted our analysis to summer as this is the season during which insolation and temperature normally peak along with most key biogeochemical processes. Mean summer DTR's differed vastly (i.e. nearly 30-fold) between lakes (Fig 1B). Differences appeared unrelated to geographical proximity, as illustrated by the much greater diel cycle of Jekl Bog, Wisconsin, relative to that of neighbouring Sparkling Lake (Fig 1C). In fact, the temperature range in a single day in Jekl Bog was close to the entire annual temperature range of Sparkling Lake. The non-significant influence of geographical proximity has also been shown by Woolway et al. [17] for five lakes in the English Lake District and was also found for the 100 lakes studied here by investigating the relationship between the distances among lakes to differences in the diel changes in water temperature.
Summary output from the statistical model described in Eq 3 illustrates that lake area and I z are identified as the most important covariates related to the mean summer DTR ( Table 1). The altitude and latitude of these lakes appear to have little relationship with the DTR, especially the latter, whose contribution in the model had been shrunk considerably and is effectively zero. The fitted smooth functions illustrate how the effect on the response varies over the observed range of the covariate (Fig 2). The model suggests that the DTR is affected by I z only in the least transparent waters, at values of I z < 30%, but is affected by lake surface area throughout. Thus, surface area was by far the most significant explanatory variable, although I z was also significant, while latitude and altitude were not significant when lake area and I z were taken into account. This result was in agreement with that of the multi-model inference method (Table 2). However, calculated AICc weights indicated that, for each of the predictor variables under consideration, no single statistical model received overwhelming support for explaining the differences in the DTR among the lakes. Rather, sets of top-ranking models received similar levels of support. More importantly though, the top models consistently included a lake surface area effect (Table 2), and for the candidate set of models, lake surface area had the highest importance, followed by I z , latitude and altitude. This is not to say, however, that latitude and altitude would not influence the DTR if a larger distribution of lakes were examined. However, for the lakes studied here, their influence was not statistically significant.
Our analysis demonstrated clearly that for the 100 lakes analysed in this investigation, lake area is the principal determinant of the DTR. Therefore, we used a separate model to describe the relationship between the DTR and lake surface area alone. Similar to Eq 3, a GAM with a gamma error distribution was chosen to model the DTR as a function of lake surface area. For this model, 100 lakes were used. Lake area alone explained over 80% of the variation in the mean summer DTR (R 2 = 0.83, P < 0.001, n = 100), which was only slightly lower than that explained (81.5%) by the model that also included the other predictor variables (i.e. I z , altitude, and latitude). The DTR was found to increase sharply with decreasing lake surface area (Fig 3). This was also shown by our theoretical approximation of the DTR (Eq 1; R 2 = 0.88 between theoretical and observed DTR), which was calculated by utilising the further subset of 24 lakes for which continuous high-resolution meteorological and vertical temperature profile data were available, and assuming that heat absorbed within the upper mixed layer is equally Table 1. Summary output from the fitted statistical model. Summary of the model used to describe the influence of surface area (A 0 ), the percent transmission per metre (I z ), altitude above sea level (h), and latitude (φ), as shown in Eq 3, on the diel surface temperature range. EDF is the effective degrees of freedom for the spline representing each covariate. Ref. DF is the reference degrees of freedom used in the statistical test of "no effect" for each smooth. F is the test statistic and p the approximate p-value of the test. I z is the percent transmission per meter. distributed. The relationship between lake area and the theoretical DTR followed a pattern that was statistically indistinguishable from that calculated with observed data (RMSE = 0.6°C; Fig  3). This suggests that systematic variation of upper mixed depth with lake size is the dominant factor determining the DTR in these 100 lakes, far outweighing geographical influences on heat fluxes. This was also demonstrated by Woolway et al. [17] and is illustrated in Eq 1, which indicates that the DTR would be proportional to the reciprocal of the mean mixed depth. Mixed depth is known to increase with lake area [32], but diminishingly so as area increases. The results therefore indicate that, for small lakes, the change of upper mixed depth with lake area has more influence on the DTR than does the geographical variation in daily heat flux, suggesting fundamental differences in how surface heating is distributed in small lakes.
To determine the lake surface areas at which a significant increase in the DTR occurs, a finite difference method was used to estimate the first derivatives of the fitted model (i.e. the model with only lake area as a predictor variable). This estimates the rate of change and, specifically, demonstrates when the rate of change is distinguishable from zero given the uncertainty in the model. The first derivative of the model identifies one clear region of statistically significant change in the DTR, which occurs at 3.2 x 10 6 m 2 (Fig 4). This period of significant change in the DTR was then superimposed on the observed data (e.g. Fig 3).
As well as exploring the influence of lake surface area, which was calculated to have the most significant effect on the DTR, we also tested other morphological characteristics of the lakes, such as volume and maximum depth, which are frequently covariant with lake surface Table 2. Summary output from the multi-model inference approach. The relative contributions of surface area (A 0 ), the percent transmission per metre (I z ), altitude above sea level (h), and latitude (φ) are shown. Confidence set of models ranked according to their adjusted Akaike Information Criterion (AICc) statistic.  area. However, these were all statistically less significant in describing variation in the DTR. For example, by using maximum depth instead of lake area, our statistical model explained only 73% of the variations in mean summer DTR; substantially less than the original model form. Due to the high correlation between lake volume and lake surface area (R 2 = 0.96), substituting volume into Eq 3 did not alter the model output by much, although the deviance explained by the model was slightly lower. However, re-calculating the theoretical DTR (Eq 1) using lake volume instead of V zmix demonstrated that the RMSE of the theoretical calculation increased from 0.6°C (i.e., using V zmix ) to 1.53°C (i.e., using the entire lake volume), indicating that the DTR was more strongly related to the volume of the upper mixed layer than the lake volume.

Model
Our dataset demonstrated that lake area exerted a significant influence on the DTR in lakes smaller than 3.2 x 10 6 m 2 (Figs 3 and 4). Small lakes vastly outnumber large lakes globally [33] and it has been estimated that 77% of lakes are of the order of 1 x 10 4 m 2 or smaller. Our results therefore emphasise that, in contrast to many of the world's most studied lakes that tend to be large and thus strongly buffered thermally, the majority of lakes worldwide undergo marked diel variations. The DTR on individual days within the summer season can be significantly greater still. Fig 1C, for example, shows that the small Jekl Bog (surface area 2.5 x 10 3 m 2 ) had a DTR that can extend beyond 15°C on some days.
Given the importance of temperature in influencing biogeochemical processes and the distribution of many organisms [34], our findings highlight the likely importance of a lake's dimensions in determining its ecological structure and function. Several studies have emphasised the role of lake size in shaping lake behaviour [35], including carbon cycling [8,36] and species diversity [37], but the possibility that diel temperature variation, mediated by lake size, may be important in these differences and more generally has not previously been recognised. For example, large surface diel temperature cycles may have important implications for the assessment of temperature-dependent biogeochemical cycles influenced by lakes. In small lakes these processes need to be studied in a manner that accurately resolves the influence of temperature at diel timescales.
Many ecological and biochemical processes are non-linearly dependent on temperature, as illustrated by the ubiquitous use of Q 10 values for representing rate changes. Literature values derived from analysis of global data sets suggest Q 10 values of about 1.6 and 2.5 for photosynthesis and respiration respectively [5] and about 4 for methane emissions [10]. Where diel temperature cycles are large, but not resolved in data collection, errors in assessing the magnitude of processes may, therefore, also be large. If, for example, a DTR is only 1°C, as is typical of large lakes, estimation of daily means of solubility of O 2 and CO 2 and rates of processes with a Q 10 of 2 or 4 from single point measurements could lead to errors of approximately 1% for O 2 and CO 2 , and up to 4 and 7% for processes with a Q 10 of 2 or 4, respectively (Fig 5). On the other hand, for small lakes, with diel cycles of, say, 7°C, errors could be up to 7 and 9% in the near-surface solubility of O 2 and CO 2 , respectively, and up to 28 and 62% for the rate of processes with a Q 10 of 2 and 4, respectively (Fig 5). On extreme days, such as those illustrated for Jekl Bog (Fig 1C) where diel temperature cycles can vary by up to 15°C, estimated rates for processes with a Q 10 of 4 from a single measurement could be in error by over 180%. Not including diel variations in temperature within gas flux calculations could even lead to erroneous estimates of the direction of gas flux as well as the magnitude. In addition, single point measurements taken at variable times of day may not provide sufficiently representative values to make cross-lake comparisons or robust assessments of long-term trends in surface temperature (e.g. [38]). This is a key area of research that requires attention when conducting large-scale comparisons in lake temperature trends. This work, therefore, underlines the importance of conducting measurements at an appropriate scale in order to capture the full range of response of dynamic systems, such as small lakes, and highlights the need to increase research effort on these common but understudied systems.
We have shown, for the first time, dramatic differences between the diel upper water temperature cycles from a large distribution of lakes. While a number of factors can potentially influence the DTR, our observation that a single variable, surface area, exerts such a dominant influence on surface temperature variation will be of value in scaling up observations of impacts on temperature-dependent processes, and in elucidating systematic differences between large and small lakes. Light attenuation was also shown to influence the DTR, albeit much less strongly than to lake area. Light attenuation has, however, been shown not only to increase with decreasing lake size [39,40], but also to have more effect on the depth of the epilimnion in smaller lakes than in larger ones [41], for which, except in the very clearest lakes, most of the solar radiation is likely to be absorbed in the mixed layer. Thus, even the influence of light attenuation on the DTR will be mediated by lake area. Increases in DOC have been experienced in many lakes over recent decades, probably as a result of recovery from acidification [42], and climate change is forecast to increase DOC concentration further [43]. Similarly, ongoing cultural eutrophication and remediation will alter phytoplankton concentrations in lakes and therefore also their light climate [20]. Accordingly, the DTR may change in the future as lakes become more or less transparent than today and mixing depths alter correspondingly. Similarly, as wind provides mixing energy for lakes and alters turbulent fluxes of heat then future changes that affect wind speed, such as development or deforestation (e.g. [44]), also have the potential to alter the DTR.

Conclusions
The diel range in lake surface temperature has a potentially major bearing on lake biogeochemical and ecological processes but, until the recent proliferation of high-frequency temperature measurements, large-scale analysis of diel temperature variations was not possible. We show that in summer the average diel temperature range can reach 7°C in small lakes (15°C on individual days) and that the magnitude of the diel range decreases strongly with increasing lake area. This has the potential to be a major source of uncertainty for current estimates of important temperature-dependent ecological and biogeochemical processes in lakes.  Table. General characteristics of the lakes studied in this investigation. Shown are the names of each lake, along with their calculated diel temperature range and mixed depth, their latitude (φ), surface area (A 0 ), altitude above sea level (h), maximum depth (z m ), and light attenuation coefficient (K d ). Details of the measurements available for each lake are also provided. (DOCX) S2 Table. Contact name for each lake included in this investigation. Contact information for any of the below listed data contributors can be obtained from R. Iestyn Woolway (riwool-way@gmail.com). (DOCX) publish, or preparation of the manuscript. This does not alter our adherence to PLOS ONE policies on sharing data and materials'. This work also benefited from participation in the Global