A simple lumped model to convert air temperature into surface water temperature in lakes

Introduction Conclusions References


Introduction
Water temperature is crucial for lakes physical, chemical and biological dynamics: temperature is the primary driver of the vertical stratification, and thus directly affects vertical exchanges of mass, energy and momentum within the water column.Water temperature plays a key role in influencing the aquatic ecosystem of lakes, which usually adapts to a specific range of physical and environmental conditions.As a matter of fact, water temperature can affect both the chemical (e.g., dissolved oxygen concentration) and biological (e.g., fish growth) processes occurring in the water body (e.g., Wetzel, 2001).Recent studies demonstrate that lakes are highly sensitive to climate, and their physical, chemical and biological properties respond rapidly to climate-related changes (Adrian et al., 2009).In the light of these considerations, it is evident that any significant modification to current environmental conditions may influence the limnic system, with direct impacts on the composition and richness of its ecosystem (MacKay et al., 2009).There are indeed several reasons to look for a reliable tool to have information about the dependence of water temperature on the various factors influencing the heat balance of the lake compartments.
Water temperature in lakes follows complex dynamics and is the result of a combination of different fluxes, whose sum is often small compared to the single terms (e.g., Imboden and Wüest, 1995).Therefore, relatively small errors in the estimate of the single contributions may result Published by Copernicus Publications on behalf of the European Geosciences Union.

S. Piccolroaz et al.: Air to water temperature in lakes
in a significantly large error in the evaluation of the net heat flux.This is particularly true for the well-mixed surface layer, usually termed as epilimnion during stratified conditions, which experiences strong oscillations at a variety of temporal scales: from short (hourly and daily) to long (annual and interannual) up to climatic (decades to centuries).Closing the heat balance correctly at the different scales and predicting the future trend of surface water temperature is therefore challenging, and not always possible (e.g., if meteorological data are not sufficient).As a consequence, some hydrodynamic lake models prescribe surface water temperature as surface boundary condition instead of computing the net heat flux at the water-atmosphere interface (e.g., Goudsmit et al., 2002;Piccolroaz and Toffolon, 2013;Piccolroaz, 2013).In general, large uncertainties are associated with the estimates of the various heat exchange components; however, the variables involved in the different processes are either not all independent from each other or do not present strong interannual variations, suggesting that some simplifications can be possibly adopted.For instance, shortwave solar radiation substantially depends on the latitude of the lake and on cloudiness, with the former presenting a rather regular annual trend and the latter being important mainly at short timescales (from hourly to weekly).Deep-water temperature typically changes on timescales much longer than surface water, and thus the heat exchanged with the hypolimnion can be reasonably assumed as constant in many situations.Differently, air temperature is a significant index of the overall meteorological conditions, and can be reasonably assumed as the main variable influencing the heat balance of the surface layer of the lake (Livingstone and Padisák, 2007).
Thankfully, long-term, high-resolution air temperature observational datasets are in general available, both for historical periods adopted to calibrate general circulation models (GCMs) and regional climate models (RCMs), and for future periods where air temperature is a variable commonly derived from GCM or RCM projections.Conversely, water temperature measurements are far less available, and future projections could be only obtained through the adoption of predictive models fully coupled with atmospheric and land surface models, which at the present stage is not a common practice (MacKay et al., 2009).In order to overcome these limitations (i.e., scarce availability and difficult estimation), several simple models have been formulated that use air temperature (widely accessible both for past and future periods) to derive surface water temperature of lakes.
Regression models (see Sharma et al., 2008, for a review) are typically adopted for this scope, but their use may be questionable, especially when it is necessary to extrapolate temperature values beyond the maximum (or minimum) limit of the measured time series.This is often the case in climate change studies, where the regression relationships built upon current climate condition are applied to estimate surface water temperature for different climate change scenarios, with the possibility that the projected air temperature may cover a wider interval of values.Regression-type models, either linear or non-linear, have been successfully applied to estimate the temperature of rivers and streams, giving rise to a rich literature (e.g., Kothandaraman and Evans, 1972;Crisp and Howson, 1982;Webb et al., 2003;Benyahya et al., 2007;Morrill et al., 2005).Nevertheless, these models cannot straightforwardly be extended to the case of lakes, especially for those water basins that have a significant seasonal hysteresis.As a matter of fact, the variety of processes of heat exchange across the lake surface and the thermal inertia of the water mass cause an annual phase lag between air and water temperatures, which is hard to consider in regressions.In many cases, simplistic linear regressions are adopted for the conversion, following the assumption of a direct monotonic relationship between air temperature and surface water temperature (e.g., Shuter et al., 1983;Livingstone and Lotter, 1998;Livingstone et al., 1999), which do not allow for capturing the hysteresis cycle.In other cases, seasonal hysteresis is solved by estimating different seasonal regression relationships, one for each branch of the hysteresis loop (e.g., one for the ascending and another for the descending branch) (e.g., Webb, 1974), or by using linear regressions to estimate the monthly means of surface water temperature from the monthly means of measured air temperature data (McCombie, 1959).
Besides regression analysis, water temperature of lakes can be estimated by means of process-based numerical models (e.g., Arhonditsis et al., 2004;Fang and Stefan, 1999;Peeters et al., 2002;Martynov et al., 2010), possibly coupled with an atmospheric model (e.g., Goyette and Perroud, 2012;Martynov et al., 2012) aimed at including the mutual interaction between water and atmosphere.These kinds of models can provide exhaustive information about the thermal structure of lakes, and accurately characterize the different energy fluxes involved in the lake temperature dynamics.The major drawback of the process-based models is the requirement of detailed time series of meteorological data in input (e.g., wind speed, humidity, cloudiness, etc., in addition to air temperature), which are often not available or not accurate enough.
In order to overcome the limitations of traditional approaches (both regression-and process-based models), semiempirical models based on physical principles may represent a valid alternative, having the key advantage of requiring fewer data in input than deterministic models, whilst preserving a clear physical basis.Recently, Kettle et al. (2004) proposed a simple empirical model to estimate mean daily water temperature, using only air temperature and the theoretical clear-sky solar radiation as input information.The model is based on the sensible heat exchange model of Rodhe (1952) (see also Bilello, 1964), and implicitly accounts for the main heat exchange processes through four parameters.The model performs well, but its application is limited to specific periods of the year with nearly uniform stratification conditions (it has been applied to small lakes in Greenland from late June to early September, long after ice melts, when water temperature is always above 4 • C).However, it does not account for the seasonal evolution of the thermal structure of the lake, and hence of the mixing depth (i.e., the depth of the epilimnion), which determines the volume of water responding to meteorological forcing, and has a significant influence on the seasonal patterns of atmosphere-lake heat exchange.
In the attempt to reliably estimate the cycle of surface water temperature of lakes from air temperature measurements/projections only, both under past, current and projected climate conditions, a simplified model has been developed.Such a model is primarily based on the energy balance between atmosphere and lake surface water (Fig. 1), but avoids the need to take into account all heat budget terms explicitly.A simple parameterization of the seasonal evolution of the mixing depth is included in the model equations, which only depend upon air temperature.The key objective of the present work is thus the definition of a modeling framework that could allow for a consistent description of the physical principles governing lake surface temperature, and ensures general applicability of the model (e.g., over the entire year).
The paper is structured as follows.In Sect. 2 the governing equations are presented, and the model is formulated.The heat budget is presented in detail in Appendix A. Section 3 provides a general description of the data used in this study.The results concerning model calibration and validation are presented in Sect. 4 with reference to the different versions of the model and the various dataset used in this work.Results are discussed in Sect.5, and finally we draw the main conclusions in Sect.6.

Formulation of the model
The net heat flux H net in the surface layer of a lake results from the combination of the different fluxes entering and exiting the upper water volume (see Fig. 1).The main heat exchanges occur at the interface between the epilimnion and atmosphere, and between the epilimnion and deep water (i.e., hypolimnion).According to the simplifications discussed in Appendix A, a simplified version of the net heat flux can be expressed as where t is time, t yr is the duration of the year expressed in suitable time units, T a and T w are air and water temperature (expressed in [ • C] for simplicity), respectively, and c i (i from 1 to 5) are coefficients with a physical correspondence, whose definition is detailed in Appendix A. Note that Eq. ( 1) accounts for a sinusoidal annual forcing term with amplitude c 1 and phase c 2 (which result from a combination of the seasonal fluctuations of solar radiation and of sensible and latent heat fluxes), an air-water exchange term c 4 (T a − T w ) (where  The inset shows the location of lake Superior in North America.c 4 is a coefficient that is primarily ascribable to a transfer function of sensible heat flux), a constant term c 3 and a residual correction c 5 T w dependent on the water temperature (the last two terms combined together basically account for the contribution of the latent heat flux).The only meteorological variable explicitly included in the model is T a , while the remaining meteorological forcing (e.g., wind speed, solar radiation, humidity and cloudiness, which besides air temperature are the major factors controlling the heat budget of lakes) are inherently accounted for in the model's parameters.In particular, the formulation of the model implicitly accounts for the seasonal patterns of these external forcing terms through the data-driven calibration of the parameters.Consistent with the main aim of the model (i.e., to reproduce the evolution of T w at long timescales, e.g., monthly, annual, interannual), higher frequency fluctuations are not considered.
Considering the upper layer of the lake, the volumeintegrated heat equation can be expressed as follows where ρ is the water density, c p is the specific heat at constant pressure and V and A are the volume and the surface area of the layer, respectively.Both V and A can be left undetermined in the analysis if we define the depth of the wellmixed surface layer (i.e., the epilimnion thickness) as follows The depth D typically depends on the stratification of the water column, and is characterized by a clear seasonal behavior.In order to include this essential feature, a suitable parameterization of D in time is required.The most appropriate choice is to estimate D as a function of the thermal stratification, and thus of the vertical temperature gradient.As a first approximation, the strength of the stratification can be evaluated as proportional to the difference between S. Piccolroaz et al.: Air to water temperature in lakes the surface water temperature T w and a reference value T r .T r is representative of deep-water temperature, and thus can be suitably chosen depending on the thermal regime of the lake (for a classification of lakes refer to, e.g., Hutchinson and Löffler, 1956;Lewis, 1983).In the case of cold monomictic lake (i.e., never over the temperature of maximum density T ρ,max ≈ 4 • C, stably stratified in winter and circulating in summer; in these lakes the thermal stratification is referred to as inverse since water temperature at the surface, T w , is colder than in the hypolimnion), T r can be assumed as the maximum surface temperature registered during the year.In the case of warm monomictic lakes (i.e., always above 4 • C, circulating in winter and stably stratified in summer; in these lakes the thermal stratification is referred to as direct since water temperature at the surface, T w , is warmer than in the hypolimnion), T r can be assumed as the minimum surface temperature registered during the year.Finally T r can be assumed equal to the temperature of maximum density T ρ,max ≈ 4 • C for the case of dimictic lakes (i.e., inversely stratified in winter, stably stratified in summer and circulating twice a year at the transition between the two states, at about 4 • C).In all cases, when the water column is nearly isothermal (i.e., T w → T r ), close to the onset of the seasonal turnover, the thermal stratification weakens and, as a consequence, the surface mixed layer D reaches its maximum thickness D r .In general, we assume that the stronger the stratification (i.e., |T w − T r | 0), the thinner the surface mixed layer D.
For the period of direct stratification (T w > T r ), the evolution of D is described using the simple exponential decay law where τ warm [ • C] is the inverse of the decay rate (the subscript "warm" refers to the case T w > T r : direct stratification) and D r indicates the maximum thickness of the mixed layer.
With the aim to consider the variation of D when the lake is inversely stratified (i.e., T w < T r , the subscript "cold"), a modified version of Eq. ( 4) has been derived, where τ cold [ • C] and τ ice [ • C] are the inverse of decay rates.
In principle, τ cold is not necessarily equal to τ warm , since the evolution of D below and above T r is possibly different.In addition, the second term in the exponential function has been introduced to account for the potential formation of the ice cover at the surface.As T w tends to 0 with where the model parameters p i (i = 1, 3, 4, 5) are the coefficients c i present in Eq. ( 1) divided by the product ρc p D r , p 2 = c 2 /t yr , p 6 = τ warm , p 7 = τ cold , p 8 = τ ice , and δ = D/D r is the normalized depth, whose seasonal evolution is schematically represented in Fig. 2 for the case of dimictic lakes (monomictic regimes are particular cases of the dimictic regime, which represents the most general case).
In conclusion, we propose a semi-empirical lumped model, which solves the temporal evolution of surface water temperature of lakes, using only air temperature as input forcing.The model requires the calibration of eight physically based parameters, whose range of variation can be reasonably estimated according to their definitions (see Appendix A).

Study site
In order to apply the model described in Sect.2, only two series of data are required: air temperature as input forcing, and surface water temperature for calibration purpose.A sufficiently long dataset (i.e., more than one year) is an essential prerequisite to perform a robust model calibration and validation procedure.Moreover, a long-term dataset provides a clear picture of the possible interannual temperature variability, thus allowing for the identification of a set of parameters that is appropriate to investigate long-term climate dynamics.
The model has been tested on Lake Superior (surface area: 82 103 km 2 ; volume: 12 000 km 3 ; maximum depth: 406 m), the largest of the five Great Lakes of North America (see Fig. 3), and the largest freshwater basin on Earth by surface area.Lake Superior is a dimictic lake: the temperature of the epilimnion is warmer than 4 • C in summer and cools below 4 • C in winter.While the surface water temperature varies seasonally, the temperature of the hypolimnion is almost constant over the year at about 4 • C. Twice a year, in December and in June, surface water reaches this temperature, and thus the thermal stratification weakens.Under these conditions, and in the presence of a sufficiently strong wind blowing at the surface, the entire lake can mix (i.e., lake's turnover).The inset shows the location of lake Superior in North America.

Fig. 2.
Seasonal evolution of the dimensionless thickness δ of the surface well-mixed layer for the general case of a dimictic lake.
Long-term temperature data (both for air and surface water) have been obtained from the National Data Buoy Center (NDBC) and from the Great Lakes Environmental Research Laboratory (GLERL), which are part of the National Oceanic and Atmospheric Administration (NOAA).In particular, the NDBC provides historical meteorological and oceanographic data for a network of offshore buoys and coastal stations (i.e., Coastal Marine Automated Network -C-MAN) that is installed all over the world, while the GLERL, through the CoastWatch program, releases daily digital maps of the Great Lakes' surface water temperature and ice cover (i.e., the Great Lakes Surface Environmental Analysis -GLSEA).
Concerning the NDBC dataset, two different stations have been used in this work: (a) 45004 -Marquette, an offshore mooring buoy that provides water temperature measured at 1 m below the water surface; and (b) STDM4 -Stannard Rock, a C-MAN station installed on a lighthouse that provides air temperature series measured at about 35 m above the lake surface.These two stations have been chosen from the many that are available for Lake Superior (both offshore buoys and C-MAN) because of their central location (see Fig. 3) and long-term data availability, but time series registered in other stations present similar behavior (not presented here).The observational dataset cover a 27 yr long period, from 1985 to 2011, and consists of measurements with 1 h temporal resolution.Since NDBC buoys in the Great Lakes are removed during winter to prevent damage from ice, no measurements are available at the 45004 -Marquette during winter months (except for 1991), while the STDM4 -Stannard Rock measurements do not show systematic gaps.
Concerning GLERL dataset, daily temperatures have been used for the period 1994 to 2011.Data refer to the daily lake average surface water temperature obtained from NOAA polar-orbiting satellite imagery.The series does not present systematic gaps (missing data, see Table 1, are concentrated in the initial warm-up year, and hence do not contribute to the evaluation of the model efficiency; see Sect.4.1), thus providing surface water temperature also in winter, which, on the contrary, is almost completely uncovered by the NDBC dataset.A mismatch between NDBC and GLERL datasets is The inset shows the location of lake Superior in North America.visible in the rising limb of the annual cycle of temperature (i.e., between April and July; see Fig. 8a), which is likely to be a consequence of the different spatial scales of the two series of data: while the NDBC dataset represents surface water temperature measured nearly at the center of the basin, the GLERL dataset provides values averaged over the whole lake.In the latter case, the spatial variability of surface water temperature (e.g., in spring, lake water heats from the shores towards the offshore deeper zones) is intrinsically included in the estimates, thus determining smoother annual cycles of temperature.Despite this discrepancy, Schwab et al. (1999) compared GLERL data with measurements at some of the NDBC buoys, finding overall good agreement.In particular, for the case of the 45004 -Marquette buoy used in this work, the mean difference between the two datasets for the period 1992-1997 is less than 0.28 • C, the root-mean-square error is 1.10 • C and the correlation coefficient is 0.96.As customary, the available datasets have been divided into two parts: the first part, containing around two-thirds of the available data, is used for model calibration and sensitivity analysis, while the second part, containing the remaining one-third, is used for model validation.Missing data in the water temperature series have not been replaced (they do not contribute to the evaluation of the efficiency of the model); on the other hand, gaps in the air temperature series have been reconstructed with estimates obtained as an average of the available data in the same day over the corresponding period (i.e., calibration or validation).The datasets used in this work are listed in Table 1, together with their main statistics.
The differential Eq. ( 6) has been solved numerically by using the Euler explicit numerical scheme (see, e.g., Butcher, 2003), with a daily time step (concerning NDBC data, mean daily temperatures have been preliminary calculated from the original data).

Sensitivity analysis and model calibration
Inverse modeling of complex systems, as such as those encountered in hydrological applications, is an inherently illposed problem, as the information provided by observational data is insufficient to identify the parameters without uncertainty.In a typical situation many different combinations of the parameters may provide similar fitting to the observational data.For example, even a simple model with only four or five parameters to be estimated may require at least 10 hydrographs for a robust calibration (e.g., Hornberger et al., 1985).This identification problem can be alleviated by reducing the number of parameters used in the model, for example through sensitivity analysis, which is the typical methodology used for this purpose (e.g., Majone et al., 2010).
In this work we perform sensitivity analysis by using generalized likelihood uncertainty estimation (GLUE), a methodology proposed by Beven and Binley (1992) that requires the identification of a validity range for each parameter, a strategy for sampling the parameter space and finally a likelihood measure to be used in order to rank the different parameter sets.We carried out 100 000 000 Monte Carlo model realizations using uniform random sampling across specified parameter ranges selected according to physical limitations of the model's parameters.Indeed, the physical meaning of the parameters allowed for a reasonable definition of the possible range of variability of each of them.Finally, we used as a likelihood measure the Nash-Sutcliffe model efficiency coefficient, E, which is a widely used metric adopted in hydrological applications (e.g., Nash and Sutcliffe, 1970;Majone et al., 2012, and many others): where n is the number of data, σ 2 e and σ 2 o are the variance of the residuals and of the observations, respectively, T w,i and T w are the observed and simulated surface water temperature at time t i , and T w is the average of T w,i .Note that the residual is defined as the difference between the observational data and the model's prediction, and a parameter set identifies a point in the space of parameters.The Nash-Sutcliffe efficiency index ranges from −∞ to 1.An efficiency equal to 1 (E = 1) corresponds to a perfect match between measured and simulated values, whilst E = 0 indicates that the model prediction is as accurate as the mean of observations.Efficiency values lower than 0 (E < 0) occur when σ 2 e is larger than σ 2 o , and thus when the mean of observations is a better estimator than the model itself.
Since its introduction in 1992, GLUE has found wide applications and it is recognized as a useful methodology for uncertainty assessment in many fields of study especially in non-ideal situations (e.g., Beven, 2006).Nevertheless, the goal of this work is not to adopt a complete, informal Bayesian approach to estimate uncertainty of model predictions but rather to assess the impact of changes in uncertain parameter values on model output.Hence, the purpose is to set up a general and effective strategy to select which are the most sensitive parameters of the model, and to define which should be suitably adjusted in a calibration process.Indeed, the GLUE methodology is a powerful and effective tool that can be also used for model calibration besides uncertainty estimation procedures and sensitivity analysis.
One of the most acknowledged limitations of the GLUE methodology is the dependence on the number of Monte Carlo simulation, especially in the presence of complex models with high computational demand.However, in our case we were able to fully explore parameter response surfaces by adopting a significantly high number of realizations.Therefore, the use of the GLUE methodology, not only for parameter identifiability purposes but also as a calibration tool, appears appropriate.Just as a side note, 100 000 000 model runs over a period of 18 yr with a daily time step and adopting Intel(R) Xeon(R) CPU X5680 at 3.33 GHz took around 2 h; the code is written in Fortran 90.
Furthermore, in the attempt to test reliability and predictive capability of the model, a validation procedure was undertaken by running the model on validation datasets (see Table 1) by using the sets of parameters that maximize the efficiency E during the calibration periods (i.e., the set of parameters with the highest likelihood obtained through the GLUE methodology).
In the ensuing sections, the GLUE methodology is presented for different model configurations (from eight to four parameters), with reference to the calibration periods of the NDBC and GLERL datasets, respectively (see Table 1).
As a final comment, we point out that the first year of each time series is used as warm-up period (i.e., excluding the Table 2.Estimated model parameters for NDBC, GLERL and GLERL myr (mean year) simulations, and their estimated physical range of variation.period from the calculation of the Nash-Sutcliffe efficiency index E) in order to remove any transient effect due to the initial condition.

Eight-parameter model
As discussed in the previous section, the performance of the model has been tested in the framework of the GLUE methodology, using the time series of air and surface water temperatures provided by the NDBC center (see Table 1).A first-step sensitivity analysis has been carried out by solving Eq. ( 6) over the calibration period (18 yr, from 1985 to 2002) using randomly sampled sets of parameters.We recall here that each of the eight parameters has been allowed to vary over physically reasonable ranges of values, which have been previously determined through Eqs.(A12)-(A16) in Appendix A by considering the lake location properties (e.g., latitude, climate, typical temperatures) and by using coefficients available in the literature (e.g., Henderson-Sellers, 1986;Imboden and Wüest, 1995;Martin and Mc-Cutcheon, 1998).As far as the exponential decay laws are concerned, parameters p 6 and p 7 have been allowed to range within a wide interval comprised between 0 and 15 • C, and p 8 between 0 and 0.5 • C. As already mentioned in Sect.2, the range of p 8 has been set narrower (more stringent upper bound), aimed at confining the correction due to ice only when the lake is inversely stratified (e −4/0.5 = O(10 −4 ) 1).In the light of the results obtained from the first-step analysis, the ranges of variability of each parameter were narrowed, thus allowing for a detailed investigation of parameter space regions associated with high values of E. Subsequently, a second-step sensitivity analysis has been undertaken by sampling further 100 000 000 sets of parameters from the narrowed parameter ranges.
Figure 4 shows the dotty plots of the efficiency index E for each of the eight parameters of the model, corresponding to the narrowed ranges (second-step performance analysis).For the sake of clarity in the presentation of results, Fig. 4 (as well as Fig. 6) shows only the set of parameters with E larger than 0.8.Model simulation during the calibration period 1985-2002 using the best set of parameters is illustrated in Fig. 5, which shows noticeable agreement between simulated and observed values, with an efficiency index E > 0.9 (see Tables 2 and 3 for a summary of the results).A visual inspection of dotty plots of model efficiency can provide useful information on the identifiability of each parameter.According to Fig. 4, all the parameters are characterized by a good identifiability with the exception of parameters p 7 and p 8 .This can be explained if we consider that the NDBC dataset presents only a few values of water temperature during winter periods.In particular, no data (except for one single year out of 18) are available when surface water temperature approaches 0 • C, which corresponds to the period of the year when the parameter p 8 is relevant.We remark here that p 8 is the parameter associated to ice formation at the surface of the lake.Although to a minor extent, also p 7 does not show a clear identifiability, probably as a consequence of the importance that this parameter assumes only during winter time (it plays a role only when T w < 4 • C); that is, when water temperature measurements are not available.However, by analyzing the dotty plot an important information can be inferred: efficiency increases for higher values of p 7 , and approaches a nearly asymptotic high-efficiency trend when p 7 5 • C. Looking at the physical meaning of the parameter (see Eq. 7), this means that model performance improves as the mixed depth D approaches its maximum value D r when the lake is inversely stratified (T w < 4 • C).
In the light of these evidences, in the ensuing section the full eight-parameter version of the model has been simplified by neglecting parameters p 7 and p 8 .

From eight to four parameters
On the basis of the results discussed in the previous section, not all the model parameters seem to be significant and clearly identifiable.In particular, the parameter p 8 has been found to be insensitive to the model, and parameter p 7 provides an overall high performance over most of its variability domain (p 7 5 • C).The peculiar behavior of these parameters, together with the fact that both appear in the definition of the mixing depth D when the lake is inversely stratified, suggests that a simplification of the model may be possible by considering a different (simpler) expression for D. The parameter p 8 is not significant for the model, and can be easily neglected, thus eliminating the effect of ice formation.On the other hand, according to Eq. ( 7), high values of p 7 mean small decay rates of D, and thus thick mixing depths when the lake is inversely stratified.In the light of these considerations, we have derived a first simplified version of the model, where the mixing depth is assumed to be constant and at its maximum thickness D = D r when the lake is inversely stratified (T w < T r ).Thanks to this simplification p 7 and p 8 are removed and the number of parameters diminishes from eight to six.
As for the case of the full eight-parameter version of the model, the same sensitivity analysis described in Sect.4.2 has been carried out also for the simplified six-parameter version.The set of parameters presenting the highest efficiency index during the calibration period is summarized in Table 2, whilst dotty plots deriving from the application of GLUE methodology and the comparison between simulated and measured surface water temperature during the same period are not presented here for the sake of brevity.Indeed, results are essentially equivalent to those obtained using the full eight-parameter model, which is confirmed by the close similarity between the best set of parameters and the efficiency indexes obtained in the two cases (see Tables 2 and  3).The similarity of results supports the idea that the sixparameter model is a reasonable simplification, at least in the case considered herein where the winter data are not abundant (see also Sect.4.5 for further discussion).
The number of parameters can be further diminished from six to four by eliminating the parameter p 1 and, as a direct efficiency is presented with an orange dot.consequence, p 2 , in addition to p 7 and p 8 .This simplification is justified since three periodic terms appear in the model as characterized by an annual periodicity: the forcing term p 1 cos(2π(t/t yr − p 2 )), the exchange term p 4 (T a − T w ) and the residual correction p 5 T w .The simultaneous co-presence of all these terms may be considered redundant in those cases in which the annual cycles of T w and/or of T a − T w can be suitably approximated as sinusoids.Indeed, the sum of sinusoidal functions with the same frequency, but different amplitude and phase, yields another sinusoid with different amplitude and phase but the same frequency.Therefore, two sinusoids are sufficient, and the forcing term p 1 cos(2π(t/t yr − p 2 )) can be removed, relieving the overall annual variations on the periodic terms controlled by the model variable T w and the external forcing T a .Following this logic, the term p 5 T w could be neglected alternatively (on the contrary, p 4 (T a −T w ) cannot be neglected since it is the only term that includes information about the external forcing), but this assumption would remove only one parameter (p 5 ) instead of two (p 1 and p 2 ), thus making it less attractive.p 1 cos(2π(t/t yr −p 2 )) and p 5 T w cannot be neglected contemporaneously, because in this case the phase of the overall periodic term would be forced to that of the temperature difference.It is worth noting that in the four-parameter version of the model (retaining p 3 , p 4 , p 5 and p 6 ), the physical meaning of the parameters is distorted, as the processes that were accounted in p 1 are now included in p 4 and p 5 .
Figure 6 shows the dotty plots of the efficiency indexes E for the four-parameter version of the model, where only the parameters p 3 , p 4 , p 5 and p 6 are retained.In this case, since the number of random samplings has been kept unchanged (i.e., 100 000 000), the predictions of the Monte Carlo realizations appear much less sparse (i.e., denser dotty plots) if compared to the case of the eight-parameter version.Notice that all the parameters are characterized by high identifiability and the model does not present signs of overparameterization. Figure 5 shows the comparison between observed and simulated surface water temperatures during the calibration period 1985-2002 for the four-parameter and the full eightparameter models, respectively.The difference is small, and mainly localized during the winter period, when no surface water temperature data are available for comparison.During the rest of the year, when measurements are available and E index can be effectively calculated, the two solutions are comparable and the efficiency indexes are similar (just slightly lower for the simplified four-parameter version of the model; see Table 3

Model validation
The best set of parameters obtained during the calibration period for the different versions of the model (with eight, six and four parameters; see Table 2) have been used to run the model during the validation period 2003-2011 (see Table 1).In all the cases, simulations have been characterized by high efficiency indexes, comparable to those obtained with the best simulations during the calibration period (E 0.9; see Table 3).In Fig. 7 simulated water temperatures (for the versions with the eight and four parameters) are compared with observations showing overall very good agreement.Results confirm the reliability of the model as a valuable tool for surface water estimation over long-term periods with different model configurations.
Besides the evaluation of the surface water temperature, the model provides additional qualitative information regarding the annual evolution of the epilimnion thickness.In fact, as discussed in Sect.2, the model explicitly includes a simplified parameterization of the seasonal behavior of the mixing depth through Eq. ( 7).In particular, the normalized thickness of the epilimnion δ = D/D r is automatically determined once the parameters p 6 , p 7 and p 8 are defined.Furthermore, it is evident that if an estimate of the reference mixing depth D r were known, the actual thickness of the well-mixed layer D could be evaluated as well.With reference to the NDBC dataset, Fig. 8b shows the evolution of δ over the validation period 2003-2011 for the eight-and four-parameter versions of the model (continuous lines).In the first case (eight parameters) the fictitious increase of depth due to the presence of ice is evident (peaks at values greater than 1).Such an increase is related to the presence of a larger water volume involved in the heat balance (see Fig. 2).Thanks to this assumption the model accounts for the insulation effect due to the presence of ice, which may be even more significant when the ice surface is snow covered and the penetration of solar radiation is strongly attenuated.

Satellite data
We have seen that results are remarkable both using the full eight-parameter version and the simplified six-and fourparameter versions of the model.In particular, little difference has been found regarding the best set of parameters and the efficiency indexes obtained using the eight-and the six-parameter versions (see Tables 2 and 3).Therefore, one may infer that no significant advantages can be expected by using a more accurate expression for D during the winter period (eight-parameter version), instead of a constant value (six-parameter version).However, it is not possible to state this conclusion by simply analyzing results presented in Sects.4.2, 4.3 and 4.4.Indeed, no surface water temperature measurements are available during the winter period for the NDBC dataset except for the year 1991, and thus the model efficiency has not been tested during the period of inverse stratification.
Aimed at overcoming this limitation, the GLERL dataset has been used, which provides daily lake-averaged surface water temperature based on satellite imagery, and covers the whole year (see Table 1).The same GLUE procedure discussed in Sect.4.1 has been performed by adopting the GLERL dataset as reference surface water temperature data and by repeating the implementation details described in Sect.4.2.The simulations (hereafter referred to as GLERL) have been run over the calibration period 1994-2005, using as input forcing the air temperature data retrieved from the NDBC dataset (C-MAN station; see Table 1).The eight-, sixand four-parameter versions of the model have been tested, obtaining remarkable efficiency indexes (E > 0.95), which are higher with respect to the previous applications (i.e., using the NDBC dataset).A validation procedure has been conducted with reference to the period 2006-2011, confirming high performances of the model (E > 0.97).The parameter sets providing the highest efficiencies during the calibration period and the associated E values are given in Tables 2 and  0    3, while comparison between simulated and observed surface water temperature data for the eight-and four-parameter versions of the model are shown in Figs. 9 and 10 for the calibration and validation periods, respectively.Finally, the seasonal evolution of δ is shown in Fig. 8b for both versions of the model (dashed lines), and it is compared with those obtained from the application of the model to the NDBC dataset.Results are consistent, with the slight difference in the onset of summer stratification being due to the earlier increase of water temperature in the GLERL dataset compared to the NDBC dataset (see Fig. 8a).

Discussion
The physically based, semi-empirical model presented here has been shown to provide an accurate description of surface water temperature of lakes, with high values of Nash-Sutcliffe efficiency index E 0.9, and a root-mean-square error between observations and simulations of the order of 1 • C (see Table 3).This error in prediction capability is comparable to those obtainable using process-based numerical models (e.g., Fang and Stefan, 1996;Stefan et al., 1998), which, however, have the strong limitation of requiring highresolution weather data and the calibration of numerous internal parameters.
The close agreement between measurements and model estimates is further confirmed in Figs.11a and 11b, which illustrate the parity diagrams for monthly averaged surface water temperature during the calibration and validation periods of GLERL simulation, respectively.No systematic deviation (bias) is observed, and the dispersion along the diagonal does not exhibit significant trends.Both these characteristics are confirmed by the small values of mean error (ME) and rootmean-square error (RMSE) listed in Table 3. Figures 11a and  11b also illustrate that the model is able to adequately describe interannual fluctuations, as is indicated by the range of variability of monthly averaged temperatures associated to the coldest (March, blue dots) and warmest (August, red dots) months.This evidence is also confirmed by Figs. 9 and 10, where the model coherently reproduces the occurrence of relatively colder (e.g., 2004) and warmer (e.g., 1998) periods.
So far, the model has been tested with long-term series of data (NDBC: 27 yr; GLERL: 18 yr); however, long-term records are often not available, or are characterized by significant gaps due to missing data.Instead, it is relatively easier to have access to mean annual cycles of temperature (both of surface water and air), whose determination also represents a valuable strategy to overcome the possible lack of data.Therefore, a conversion model that could be calibrated on mean annual cycles, and successively applied over longterm periods without compromising the correct estimation of the interannual fluctuations, would represent a valuable tool.For this purpose, the mean annual cycle of surface water temperature has been derived from GLERL data during the calibration period 1994-2005, and the corresponding cycle of air temperature from the NDBC dataset (C-MAN station).A Monte Carlo sensitivity analysis (hereafter referred to as GLERL myr , the subscript "myr" standing for mean year) has been carried out following the same procedure adopted in the previous sections, but using mean annual cycles of air      and water temperature as forcing and reference data, respectively.In order to eliminate the influence of initial conditions, the temperature cycles have been replicated for two years with the first one used as a "warm-up".Results obtained by adopting the parameters providing the highest efficiency (see Table 2) are presented in Fig. 12, which shows the hysteresis cycles between air and surface water temperatures derived from measurements and model estimates (eight-and fourparameter versions).A very high efficiency index (E 1.0) is achieved, and both versions of the model are able to satisfactorily capture the seasonal pattern of thermal hysteresis.The parameter sets providing the highest efficiency during the calibration process with the mean annual temperature data and the associated E value are summarized in Tables 2  and 3.
Subsequently, a validation procedure has been conducted for GLERL myr during the period 2006-2011 (the same as GLERL simulations; thus results can be compared).Results are characterized by remarkable efficiency indexes (E 0.97), only slightly lower than the values obtained with the simulations presented in the previous sections.Indeed, the model calibrated on the mean year is able to capture the interannual variabilities well, producing remarkable results not dissimilar from those shown in Fig. 10 (for this reason they are not presented here).Furthermore, parameter values are significantly similar to those obtained calibrating the model with the whole 12 yr series of data (GLERL simulation; see Table 2).
In the light of the results presented in this section, we assert that the model can be calibrated and adopted using data of different origins (measurements at buoys and coastal stations, satellite estimates) and nature (long-term series of data, mean annual cycle of temperature).This conclusion is corroborated by the excellent results (not shown here for the sake of brevity) of the performance analysis (entirely comparable to those presented in the present paper) obtained using different datasets:  at the 45004 -Marquette offshore buoy station at only 4 m from the lake surface; and (c) water temperature measured at a different offshore buoy (the 45001 -Hancock).Furthermore, in all cases, even if the calibration is performed considering mean annual cycles of temperature, the model suitably captures the interannual variations that are likely to occur.On the basis of this evidence, we can assert that this simple model may be used with different air temperature datasets as input and, unlike process-based models, it can be calibrated using any water temperature dataset, independent of its physical representativeness (e.g., point measurements vs. spatial averages).Therefore, in principle, air temperature series provided by GCMs and RCMs can be used as well.In this regard, the model is particularly attractive for climate change impact studies, since predictions of air temperature are usually more reliable and available than meteorological variables (e.g., Gleckler et al., 2008).Based upon these considerations, Piccolroaz (2013) exploited the same approach to reproduce the current status and to predict future modifications of surface water temperature of Lake Baikal (Siberia).

Conclusions
In this work a simple, physically based model has been developed to estimate surface water temperature from air temperature.In particular, we show that our modeling framework is able to reproduce the observed water temperature data with limited information on external meteorological forcing over long timescales, ranging from monthly to interannual.
Starting from the zero-dimensional heat budget, we derived a simplified first-order differential equation for water temperature forced by a few terms representing the combined effects of the seasonally varying external terms and the exchange terms explicitly dependent on the difference between air and water temperatures.Assuming annual sinusoidal  cycles of the main heat flux components, eight parameters have been identified, which can be calibrated if temporal series of air and surface water temperature are available.Such a calibration is supported by the physical interpretation of the parameters, which provides reasonable initial conditions for the parameters ranges.
The relative importance of the model's parameters have been evaluated by using the GLUE methodology.Thanks to this analysis we were able to identify and neglect parameters that, under different conditions, appears less significant in the model formulation, leading to two simplified versions retaining six and four parameters, respectively.
The model has been applied to the case of Lake Superior (USA-Canada) with reference to different types of datasets, and all the versions of the model have shown to perform well in reproducing the measured water temperature data.This model has proved to be robust and able to suitably simulate the lake's response to meteorological forcing, including interannual variability, representation of the variability of the epilimnion thickness and the inverse stratification process that typically occurs in dimictic lakes.In our view, Air2Water represents a valuable alternative tool to correlation models, which require the same data in input as our model but are not able to address some fundamental processes (e.g., the hysteresis cycle between air and water temperature).Furthermore, it can be used in place of full process-based models when meteorological data are not sufficient for their proper application.In principle, the simple model presented here is likely to be effectively applied to lakes with different characteristics, although some inconsistencies could arise in those cases where the assumptions on which the model formulation has been based (see Appendix A) are no longer valid (e.g., S. Piccolroaz et al.: Air to water temperature in lakes tropical lakes characterized by intense evaporation, basins in which the through-flow is consistent, lakes located in regions where the variability of meteorological forcing is significant at sub-annual frequency).
In the light of these results, the model can represent a valuable tool in climate change impact studies, allowing for predictions of future trends of lake surface water temperature, given future projections of air temperature only.Finally, it is worth noting that if the model is calibrated using air temperature series from climate models (global or regional scale) and measured records of water temperature (lake scale), a downscaling operation is implicitly implemented in the conversion procedure (Piccolroaz, 2013).

Simplified heat fluxes
Indicating with H the generic heat flux per unit surface [W m −2 ], defined as positive when it is directed towards the considered layer, the net flux is assessed accounting for the following main terms: where H s is the net shortwave radiative heat flux due to solar radiation (considering only the incoming radiation that is actually absorbed), H a is the net longwave radiation emitted from the atmosphere toward the lake, H w is the longwave radiation emitted from the water, H e is the latent heat flux (due to evaporation/condensation processes), H c is the sensible heat flux (due to convection), H p is the heat flux due to precipitation onto the water surface, H i is the effect of the throughflow of water by inlets and outlets, and H d is the heat flux exchanged with deep water.Figure 1 shows a schematic representation of the heat exchanges at the epilimnion-atmosphere and epilimnion-hypolimnion interfaces.All the components of Eq. (A1) are analyzed in detail below to point out the main variables and physical parameters involved in the heat exchange process.
The incident shortwave solar radiation approximately follows a sinusoidal annual cycle.Considering the shortwave reflectivity r s (albedo), which is a function of the solar zenith angle and of the lake surface conditions (e.g., water wave height), the net solar radiation H s reads where t is time; t yr is the duration of a year in the units of time considered in the analysis; and s 1 , s 2 , s 3 are coefficients that primarily depends on the latitude and the shadowing effects of the local topography.The effects of cloud cover, which could be accounted for by means of empirical relationships, are not explicitly considered in the present analysis.
Longwave radiation terms are calculated according to the Stefan-Boltzmann law, yielding the following formulations H a = (1 − r a ) a σ (T K + T a ) 4 , (A3) where r a is the longwave reflectivity, generally assumed to have a constant values (Henderson-Sellers, 1986); a and w are the emissivities of atmosphere and lake surface, respectively; σ is the Stefan-Boltzmann constant (5.67 10 −8 W m −2 K −4 ), T K = 273.15K; and T a and T w are the temperatures of air and water expressed in Celsius [ The emissivity w is essentially constant and close to unity, as water is nearly a black body, while a is more variable and depends on several factors among which the most important are air temperature, humidity and cloud cover (Imboden and Wüest, 1995).Air and water temperatures can be decomposed into a reference value representative of the specific case study (T a and T w ) and a fluctuation (T a and T w ).Hence, considering that T a /(T K + T a ) and T w /(T K + T w ) are small parameters, the longwave fluxes Eqs. ( A3) and (A4) can be linearized using a Taylor expansion as H w − w σ T K + T w 4 1 + 4 T w where ˜ a = (1−r a ) a .By choosing T a = T w = T (hence T a = T a −T and T w = T w −T ), the terms H a and H w can be easily combined to yield the following equation: The sensible (H c ) and latent (H e ) heat fluxes are calculated through bulk semi-empirical relations that can be derived from turbulence theory (Henderson-Sellers, 1986): H e = α e (e a − e w ) , (A9) where α c [W m −2 K −1 ] and α e [W m −2 hPa −1 ] are transfer functions primarily depending on wind speed and other meteorological parameters, e a is the vapor pressure in the atmosphere and e w is the water vapor saturation pressure at the water temperature (both in [hPa]).The ratio α c /α e is known as the Bowen coefficient, and is often taken to be constant (≈ 0.61 hPa K −1 ) (Imboden and Wüest, 1995).The saturated water pressure e w is a function of temperature, and can be

Fig. 1 .Fig. 2 .Fig. 3 .
Fig. 1.Main heat exchange affecting the surface layer.For the description of the single terms refer to Appendix A.

Fig. 1 .
Fig. 1.Main heat exchange affecting the surface layer.For the description of the single terms, refer to Appendix A.

Fig. 1 .Fig. 2 .Fig. 3 .
Fig. 1.Main heat exchange affecting the surface layer.For the description of the single terms refer to Appendix A.

Fig. 2 .Fig. 3 .
Fig. 2. Seasonal evolution of the dimensionless thickness δ of the surface well-mixed layer for the general case of a dimictic lake.

Fig. 3 .
Fig. 3. Lake Superior with the location of the NDBC stations (45004 -Marquette and STDM4 -Stannard Rock) used in this work.The inset shows the location of lake Superior in North America.

Fig. 4 .
Fig. 4. Dotty plots of efficiency indexes (E) for the 8-parameters model during the calibration period 1985-2002 (NDBC simulation).Highest efficiency is presented with an orange dot.

Fig. 5 .
Fig. 5. Comparison between simulated and observed surface water temperature during the calibration period 1985-2002 (NDBC simulation).Simulated curves refer to the full 8-parameter and the simplified 4-parameter models, respectively.Observed air temperature data are also presented with cyan line.

Fig. 4 .
Fig. 4. Dotty plots of efficiency indexes (E) for the eight-parameter model during the calibration period 1985-2002 (NDBC simulation).Highest efficiency is presented with an orange dot.

Fig. 5 .
Fig. 5. Comparison between simulated and observed surface water temperature during the calibration period 1985-2002 (NDBC simulation).Simulated curves refer to the full 8-parameter and the simplified 4-parameter models, respectively.Observed air temperature data are also presented with cyan line.

Fig. 5 .Fig. 6 .
Fig. 5. Comparison between simulated and observed surface water temperature during the calibration period 1985-2002 (NDBC simulation).Simulated curves refer to the full eight-parameter and the simplified four-parameter models.Observed air temperature data are also presented with cyan line.

Fig. 7 .
Fig. 7. Comparison between simulated and observed surface water temperature during the validation period 2003-2011 (NDBC simulation).Simulated curves refer to the full 8-parameter and the simplified 4-parameter models, respectively.Observed air temperature data are also presented with cyan line.

Fig. 6 .
Fig. 6.Dotty plots of efficiency indexes (E) for the four-parameter model during the calibration period 1985-2002 (NDBC simulation).Highest efficiency is presented with an orange dot.

Fig. 7 .
Fig. 7. Comparison between simulated and observed surface water temperature during the validation period 2003-2011 (NDBC simulation).Simulated curves refer to the full 8-parameter and the simplified 4-parameter models, respectively.Observed air temperature data are also presented with cyan line.

Fig. 7 .
Fig. 7. Comparison between simulated and observed surface water temperature during the validation period 2003-2011 (NDBC simulation).Simulated curves refer to the full eight-parameter and the simplified four-parameter models.Observed air temperature data are also presented with cyan line.

Fig. 8 .
Fig. 8. Evolution of the dimensionless depth δ over the period 2003-2011: comparison between results for the full 8-and simplified 4parameter versions of the model obtained from NDBC and GLERL simulations.a) Comparison of observed water temperature time series during the period 2003-2011 for the NDBC and GLERL dataset, respectively.b) Evolution of the dimensionless depth δ over the period 2003-2011 for the full 8-and simplified 4-parameter versions of the model (NDBC and GLERL simulations).

Fig. 9 .
Fig. 9. Comparison between simulated and observed surface water temperature during the calibration period 1994-2005 (GLERL simulation).Simulated curves refer to the full 8-parameter and the simplified 4-parameter models, respectively.Observed air temperature data are also presented with cyan line.

Fig. 8 .
Fig. 8. Evolution of the dimensionless depth δ over the period 2003-2011: comparison between results for the full eight-and simplified fourparameter versions of the model obtained from NDBC and GLERL simulations.(a) Comparison of observed water temperature time series during the period 2003-2011 for the NDBC and GLERL dataset.(b) Evolution of the dimensionless depth δ over the period 2003-2011 for the full eight-and simplified four-parameter versions of the model (NDBC and GLERL simulations).

Fig. 8 .
Fig. 8. Evolution of the dimensionless depth δ over the period 2003-2011: comparison between results for the full 8-and simplified 4parameter versions of the model obtained from NDBC and GLERL simulations.a) Comparison of observed water temperature time series during the period 2003-2011 for the NDBC and GLERL dataset, respectively.b) Evolution of the dimensionless depth δ over the period 2003-2011 for the full 8-and simplified 4-parameter versions of the model (NDBC and GLERL simulations).

Fig. 9 .
Fig. 9. Comparison between simulated and observed surface water temperature during the calibration period 1994-2005 (GLERL simulation).Simulated curves refer to the full 8-parameter and the simplified 4-parameter models, respectively.Observed air temperature data are also presented with cyan line.

Fig. 9 .
Fig. 9. Comparison between simulated and observed surface water temperature during the calibration period 1994-2005 (GLERL simulation).Simulated curves refer to the full eight-parameter and the simplified four-parameter models.Observed air temperature data are also presented with cyan line.18 S. Piccolroaz et al.: Air to water temperature in lakes

Fig. 10 .Fig. 11 .
Fig. 10.Comparison between simulated and observed surface water temperature during the validation period 2006-2011 (GLERL simulation).Simulated curves refer to the full 8-parameter and the simplified 4-parameter models, respectively.Observed air temperature data are also presented with cyan line.

Fig. 10 .
Fig. 10.Comparison between simulated and observed surface water temperature during the validation period 2006-2011 (GLERL simulation).Simulated curves refer to the full eight-parameter and the simplified four-parameter models.Observed air temperature data are also presented with cyan line.

Fig. 10 .Fig. 11 .
Fig.10.Comparison between simulated and observed surface water temperature during the validation period2006 -2011 (GLERL simulation) (GLERL simula- tion).Simulated curves refer to the full 8-parameter and the simplified 4-parameter models, respectively.Observed air temperature data are also presented with cyan line.

Fig. 11 .
Fig. 11.Parity diagram for monthly averaged surface water temperature (eight-parameter version of the model): (a) calibration and (b) validation period of the GLERL simulation.Blue dots refer to March, red dots to August and gray dots to the remaining months of the year.

Fig. 12 .
Fig. 12.Comparison of the hysteresis cycles between daily air and surface water temperatures, as derived by the data and by the 8-and 4-parameters versions of the model.Hysteresis cycles refer to the mean year, calculated over the period 1994-2005, using GLERL and NDBC data for Tw and Ta, respectively (GLERLmy simulation).

Fig. 12 .
Fig. 12.Comparison of the hysteresis cycles between daily air and surface water temperatures, as derived by the data and by the eightand four-parameter versions of the model.Hysteresis cycles refer to the mean year, calculated over the period 1994-2005, using GLERL and NDBC data for T w and T a , respectively (GLERL myr simulation).

Table 1 .
Summary of the datasets adopted in this study, and their main statistics.