Lake Surface Water Temperature Derived from 35 Years of AVHRR Sensor Data for European Lakes

: Lake surface water temperature (LSWT) is an important parameter with which to assess aquatic ecosystems and to study the lake’s response to climate change. The AVHRR archive of the University of Bern offers great potential to derive consistent LSWT data suited for the study of climate change and lake dynamics. To derive such a dataset, challenges such as orbit drift correction, non-water pixel detection, and homogenization had to be solved. The result is a dataset covering over 3.5 decades of spatial LSWT data for 26 European lakes. The validation against in-situ temperature data at 19 locations showed an uncertainty between ± 0.8 K and ± 2.0 K (standard deviation), depending on locations of the lakes. The long-term robustness of the dataset was conﬁrmed by comparing in-situ and satellite derived temperature trends, which showed no signiﬁcant difference. The ﬁnal trend analysis showed signiﬁcant LSWT warming trends at all locations (0.2 K/decade to 0.8 K/decade). A gradient of increasing trends from south-west to north-east of Europe was revealed. The strong intra-annual variability of trends indicates that single seasonal trends do not well represent the response of a lake to climate change, e.g., autumn trends are dominant in the north of Europe, whereas winter trends are dominant in the south. Intra-lake variability of trends indicates that trends at single in-situ stations do not necessarily represent the lake’s response. The LSWT dataset generated for this study gives some new and interesting insights into the response of European lakes to climate change during the last 36 years (1981–2016).


Introduction
Water temperature is an important parameter within aquatic ecosystems.The metabolic rate of organisms and the activity of bacteria and algae is governed by water temperature [1][2][3][4].Changing water temperature can alter the balance and the species composition of an aquatic ecosystem and even alter its trophic state [5][6][7].Lake dynamic processes such as thermal stratification and mixing dynamics are firmly related to water temperatures [8,9].The energy balance of large lakes is dominated by heat exchange processes at the surface [8].Lake surface water temperature (LSWT) governs energy and mass exchange between the water body and the troposphere.It is also an important parameter in mesoscale meteorology and local climate models [10][11][12].LSWT can be seen as an indicator of regional climate change [13], and several recent publications analyzed climate signals within lakes by investigating LSWT trends [14][15][16][17].Recognizing its importance, LSWT was also added to the list of essential climate variables (ECV) by the Global Climate Observing System (GCOS), which needs to be monitored [18].With these criterions, a preselection adapted from Politi et al. [19] was used first, followed by a final manual selection.This led to a selection of 26 lakes (Table 1), with latitudes reaching from 31.53°N (Djorf Torba reservoir in Algeria) to 69.04°N (Inarijärvi in Finnland), and longitudes from 21.13°W (Thingvallavatn in Iceland) to 61.90°E (Lake Onega in the Russian Federation).The exact locations, where the LSWT data was resampled for the final time series, are shown in Table The sampling points were preferentially placed in the middle of large undisturbed open water surfaces, as the priority was to get the best possible conditions for an accurate measurement.Whereas, for the validation with in-situ data, sometimes different locations within the lakes were chosen (Table 2) to be as close as possible to the in-situ measurement location.These validation points were therefore often close to the shores, which had the effect of decreasing the availability and reliability of the satellite data while increasing the comparability with the in-situ measurements.With these criterions, a preselection adapted from Politi et al. [19] was used first, followed by a final manual selection.This led to a selection of 26 lakes (Table 1), with latitudes reaching from 31.53 • N (Djorf Torba reservoir in Algeria) to 69.04 • N (Inarijärvi in Finnland), and longitudes from 21.13 • W (Thingvallavatn in Iceland) to 61.90 • E (Lake Onega in the Russian Federation).The exact locations, where the LSWT data was resampled for the final time series, are shown in Table The sampling points were preferentially placed in the middle of large undisturbed open water surfaces, as the priority was to get the best possible conditions for an accurate measurement.Whereas, for the validation with in-situ data, sometimes different locations within the lakes were chosen (Table 2) to be as close as possible to the in-situ measurement location.These validation points were therefore often close to the shores, which had the effect of decreasing the availability and reliability of the satellite data while increasing the comparability with the in-situ measurements.Table 2. List of in-situ measuring locations from which data was available for this study.Most of the measurements are either hourly or daily averaged temperatures from automatic stations.In the case of Müritz, Lake Trichonida, and the older measurements from Lake Balaton, the measurements are sporadically and manually taken at specific locations on the lake.The average measurement density varies from 12 measurements per year at Lake Trichonida (monthly measurements) to 8760 measurements per year at Bregenz, Lake Constance (hourly measurements).

The Data
The AVHRR data archive from the University of Bern [20] was used to retrieve LSWT time series.The archive contains around 130,000 scenes from the European subset (Figure 1) and includes data from all NOAA and MetOp satellites with an AVHRR-2 or AVHRR-3 sensor on board (Figure 2).The data in the archive was preprocessed (see [20] for details) and included a probability cloud mask after Musial et al. [21] for each scene.The data of NOAA-15 was omitted, because the data from this platform showed significantly increased amounts of erroneous pixels and scanlines.The data from MetOp-1 was also omitted due to inconsistencies during the pre-processing.

The Data
The AVHRR data archive from the University of Bern [20] was used to retrieve LSWT time series.The archive contains around 130,000 scenes from the European subset (Figure 1) and includes data from all NOAA and MetOp satellites with an AVHRR-2 or AVHRR-3 sensor on board (Figure 2).The data in the archive was preprocessed (see [20] for details) and included a probability cloud mask after Musial et al. [21] for each scene.The data of NOAA-15 was omitted, because the data from this platform showed significantly increased amounts of erroneous pixels and scanlines.The data from MetOp-1 was also omitted due to inconsistencies during the pre-processing.In blue, the data that was used for the present study, and, in red, the data that was not processed due to the elevated amount of inconsistent data.
The in-situ data available for this study is shown in Table 3 are, mainly, two types of measurements: homogeneous averaged measurements on a daily or hourly basis, and sporadic instant measurements at a defined point in time.For the measurements at lake Trichonida and Lake Balaton (Keszthely old), only the date (day) without exact timestamp (hh:mm) was available; therefore 12:00 local time was assumed.All the measurements are representing bulk temperatures close to the surface (taken at 0.1-1.0m).Further, most of the northern lakes are ice covered during a certain amount of time per year.For these lakes, the amount of data per year (N/Yr. in Table 2) is considerably below of the expected 365 measurements (daily averaged data) per year.
The in-situ data were filtered for outliers and physically inconsistent values, as well as checked for irregularities such as systematic temperature changes due to modification of the equipment or the location.In the case of the INRA (Institute Nationale de la Recherche Agronomique) measurement at Lake Geneva, a systematic temperature decrease of roughly 0.4 °C was detected after autumn 2007, probably due to a dislocation of the station.Therefore, only the data from 1991 to 2007 from the INRA station was used.
It is assumed that the best validation performance will result from comparison with the hourly averaged and instant measurements, whereas for the daily measurements, the error between an individual satellite measurement and the daily average temperature can theoretically be as high as the amplitude of the diurnal temperature cycle.In blue, the data that was used for the present study, and, in red, the data that was not processed due to the elevated amount of inconsistent data.

LSWT Processing Chain
The in-situ data available for this study is shown in Table 3 are, mainly, two types of measurements: homogeneous averaged measurements on a daily or hourly basis, and sporadic instant measurements at a defined point in time.For the measurements at lake Trichonida and Lake Balaton (Keszthely old), only the date (day) without exact timestamp (hh:mm) was available; therefore 12:00 local time was assumed.All the measurements are representing bulk temperatures close to the surface (taken at 0.1-1.0m).Further, most of the northern lakes are ice covered during a certain amount of time per year.For these lakes, the amount of data per year (N/Yr. in Table 2) is considerably below of the expected 365 measurements (daily averaged data) per year.
The in-situ data were filtered for outliers and physically inconsistent values, as well as checked for irregularities such as systematic temperature changes due to modification of the equipment or the location.In the case of the INRA (Institute Nationale de la Recherche Agronomique) measurement at Lake Geneva, a systematic temperature decrease of roughly 0.4 • C was detected after autumn 2007, probably due to a dislocation of the station.Therefore, only the data from 1991 to 2007 from the INRA station was used.
It is assumed that the best validation performance will result from comparison with the hourly averaged and instant measurements, whereas for the daily measurements, the error between an individual satellite measurement and the daily average temperature can theoretically be as high as the amplitude of the diurnal temperature cycle.

LSWT Processing Chain
Figure 3 shows a scheme of the main LSWT retrieval and processing steps performed for this study.For accurate LSWT retrievals, also the pre-and post-processing of the satellite data is very important, e.g., only reliable clear sky observations with accurate geocoding should be retained to ensure that the resulting temperatures correspond to the features we actually investigate.Conservative cloud and land masking are responsible for the elimination of non-water pixels.Parameters like view angles and spatial consistency of LSWT are additional indicators for the quality and reliability of the measurements.All this information is retained in a quality mask associated with the final LSWT product.
important, e.g., only reliable clear sky observations with accurate geocoding should be retained to ensure that the resulting temperatures correspond to the features we actually investigate.Conservative cloud and land masking are responsible for the elimination of non-water pixels.Parameters like view angles and spatial consistency of LSWT are additional indicators for the quality and reliability of the measurements.All this information is retained in a quality mask associated with the final LSWT product.The pre-processing of the AVHRR data used for this study is described in Hüsler et al. [20] and provides georeferenced AVHRR scenes, as well as the associated cloud mask after Musial et al. [21].The coefficients used for this study are time-independent and lake-specific coefficients, empirically matched to radiative transfer model results as described in Lieberherr et al. [22].The details of the final LSWT retrieval is described within the following subchapters.

Georeferencing Test
Even on well geocoded AVHRR scenes, the accuracy of the geolocation can vary within the same scene.Thus, an additional georeferencing test was added to the pre-processing described in Hüsler et al. [20].During the extractions of the sub-sets from the scene, an algorithm verifies and adjusts the local geocoding for each target location.The test is basically an edge detection and matching algorithm.First, an edge detection with a Laplace filter is applied to the subset of the scene, and to the subset of a reference image.Then, the correlation for each possible positioning of the AVHRR subset within the reference subset is computed, resulting in a correlation map.In the case of a correct georeferenced image, a correlation peak is found in the center of the correlation map.Otherwise, depending on the position of the correlation peak, a latitude and/or longitude shift is performed before the subset of the lake is extracted.To determine if the correlation maximum corresponds to the correct location, it is tested if the correlation peak is unique and of significant amplitude.The geolocation check is performed with AVHRR Band 2 data for daytime scenes, and with AVHRR Band 4 data for night-time scenes.The pre-processing of the AVHRR data used for this study is described in Hüsler et al. [20] and provides georeferenced AVHRR scenes, as well as the associated cloud mask after Musial et al. [21].The coefficients used for this study are time-independent and lake-specific coefficients, empirically matched to radiative transfer model results as described in Lieberherr et al. [22].The details of the final LSWT retrieval is described within the following subchapters.

Georeferencing Test
Even on well geocoded AVHRR scenes, the accuracy of the geolocation can vary within the same scene.Thus, an additional georeferencing test was added to the pre-processing described in Hüsler et al. [20].During the extractions of the sub-sets from the scene, an algorithm verifies and adjusts the local geocoding for each target location.The test is basically an edge detection and matching algorithm.First, an edge detection with a Laplace filter is applied to the subset of the scene, and to the subset of a reference image.Then, the correlation for each possible positioning of the AVHRR subset within the reference subset is computed, resulting in a correlation map.In the case of a correct georeferenced image, a correlation peak is found in the center of the correlation map.Otherwise, depending on the position of the correlation peak, a latitude and/or longitude shift is performed before the subset of the lake is extracted.To determine if the correlation maximum corresponds to the correct location, it is tested if the correlation peak is unique and of significant amplitude.The geolocation check is performed with AVHRR Band 2 data for daytime scenes, and with AVHRR Band 4 data for night-time scenes.

The LSWT Retrieval
A state-of-the-art split-window method was used to retrieve surface water temperature for this study.The idea is to make use of the different absorption characteristics of water vapor within differential spectral windows of the thermal infrared [23].Two algorithms are currently used in various applications: the Multichannel Sea Surface Temperature (MCSST) algorithm [24] assumes a linear relationship between the atmospheric transmittance of the signal and the difference in absorption of the two channels.In contrast, the Non-Linear Sea Surface Temperature (NLSST) algorithm introduced in 1990 for NOAA SST products takes the nonlinear nature of atmospheric absorption into account, which increases the accuracy especially for daytime measurements [25,26].Recent publications show no significant improvement when using the NLSST instead of the MCSST for LSWT retrieval [27][28][29].Therefore, here the MCSST approach described by Equation (1) with empirically determined split window coefficients is used to derive LSWT from AVHRR data.
a 0,1,2,3 Split-window coefficients BT 4,5 Brightness temperature from AVHRR band 4 and 5 θ Satellite view zenith angle The coefficients where derived for each location (lake specific) and satellite, based on the findings of Lieberherr et al. [22].It was assumed that the use of matchup data from a 2-year period is sufficient to obtain the infinite temporal validity of the coefficients.The used matchup data was generated with a radiative transfer model, RTTOV-11 [30], using atmospheric profiles of water vapor and temperature from the ERA-INTERIM [31] dataset.

Non-Water Pixels Elimination
Cloud and land masking is crucial for satellite-based temperature retrievals of surface lake water.Clouds are opaque to thermal infrared radiation, and even a slightly cloud-contaminated signal (sub-pixel sized clouds) can result in a considerably colder measurement of the temperature.Also, for land contaminated pixels, the signal can be biased, as the land surface temperature sometimes differs by several degrees from the water temperature.Therefore, strict cloud and land masking are necessary to eliminate affected pixels.On the other hand, this reduces the amount of available data to be used for further processing.An ideal cloud and land masking would eliminate all affected pixels without any false positives.As, of course, this is not feasible, a good trade-off tailored to the application must be found.For the cloud detection within this study, the probability cloud mask developed by [21] is used.Additionally, in order to improve the performance over inland-waterbodies and during night-time, some additional tests are performed.These additional tests are described in [32].The masking of land pixels is done by applying a land-water mask (LWM).
The LWM is derived from freely available Open Street Map (OSM) data, which was found to be more accurate in geolocation than several other land classification datasets.The accuracy of OSM data was the subject of various studies, e.g., [33] reported an accuracy of ±3 m (SD) in an urban area in Germany, and [34] found an RMSE of 15 m in a remote rural area in Portugal.Further, the reliability of OSM data is considered good, since it is based on in-situ GPS measurements and is verified and reviewed by the OSM community.For easily accessible and frequently visited objects such as Lake Constance, the community is very active, and the data is highly reliable.When rasterizing the OSM-shapefile to the AVHRR pixel grid, the LWM is extended by 500 m (~1/2 pixel) into the lake to exclude lake pixels where the signal is potentially land contaminated.The low number of lakes (26) allowed us to check the thematic and geographic accuracy manually during the generation of the LWM raster files from the OSM water shapes.

Quality Levels
The reliability of each measurement (pixel) is expressed by a quality score reaching from Q0 to Q6; Table 3 recapitulates the six quality levels.The tests for the quality levels are cumulative, meaning that for each level, the tests of the previous quality levels have been passed too.Pixels with physically inconsistent values (LSWT > 35 • C or LSWT < −5 • C), as well as isolated pixels, are assigned to Q0; all other clear-sky pixels and water pixels are classified at least QPixels with a low spatial standard deviation (similarity to surrounding pixels) get QPixels outside the potential sun glint zone are assigned QQ4 and Q5 are achieved if the view zenith angle (VZA) is lower than 55 • or 45 • , respectively.Q6 is assigned if the georeferentiation is considered highly reliable.This means that during the georeferencing test, the feature (lake) must be detected with certainty.The results showed that at Q5 and Q6 (VZA < 45 and the enhanced georeferentiation test), the available pixels decrease strongly without considerably reducing the uncertainty.Finally, Q4 was used for the trend analysis during this study, as it was considered to be the best tradeoff between accuracy and data availability.

Time Series Generation
Once all the AVHRR scenes are processed, the 1D temperature timeseries for each specific target point are extracted from the 2D LSWT product.For each target point, pixels within a radius of 3 km are extracted.All invalid pixels and all pixels with a quality below the aimed quality level are removed.If at least 3 pixels remain within the 3 km radius around the target, the average of the valid pixels is retained.The thus-sampled temperature value is then integrated into the raw time series with the timestamp of the satellite scene accorded.
Before using this raw time series for validation and for trend analysis, an outlier detection adapted for such time-series is applied.The main problem with basic outlier detection on temperature time series is that the yearly and the diurnal temperature cycle already induces a high amount of variance in the data.For example, a measurement of 5 • C at Lake Constance in the middle of August should be detected as an outlier, as it is 10 • C to 20 • C below the expected temperature for August.However, since temperatures of 5 • C are usual in January, the measurement is not conspicuous when performing straightforward outlier detection.For robust outlier detection on a temperature time series, the known temperature cycles such as the yearly temperature cycle (YTC) must be removed first.In our case, this is done by computing the mean annual curve from the raw time series and then subtracting it from the respective values of the time-series.The YTC-corrected time series can then be analyzed.For the outlier detection, we have used a robust outlier detection method described in Leys et al. [35] using the mean absolute deviation and the median.
For the analysis of trends, the raw time series is not suitable, as neither the data density nor the overpass-time where the temperature is measured is stable through time.This can induce additional trends.In Figure 4, one can nicely see the evolving data density through the years: In the 1980s, only one satellite at the time was available, whereas in 2013 there were five.Also, the heterogeneity of overpass-times in general, as well as the shift in overpass time due to orbit drift of the individual satellites, can be also be seen in Figure 4.
To overcome the data density issue, monthly and seasonal averages of the YTC anomaly are computed.The YTC anomaly is computed with reference to lake-specific YTC issued form the whole period of available LSWT data .The averaging over months or seasons leads to a homogeneous distribution of the temperature data suited for trend analysis.By using the YTC anomaly instead of the LSWT, we aim at removing an eventual bias induced by uneven distribution of the measurements within the averaging period (month or season).
Finally, to solve the overpass-time problem, a diurnal temperature cycle (DTC) model was used to adjust each measurement to the daily mean temperature, as described in the following paragraphs.
anomaly instead of the LSWT, we aim at removing an eventual bias induced by uneven distribution of the measurements within the averaging period (month or season).
Finally, to solve the overpass-time problem, a diurnal temperature cycle (DTC) model was used to adjust each measurement to the daily mean temperature, as described in the following paragraphs.

Diurnal Temperature Cycle (DTC) Correction
The NOAA satellites are subject to relatively strong orbit drifts [36,37].This means that through the years of operation, the overpass time at the same location is shifted.Figure 4 shows all scenes and recording times used during this study.One can see that the overpass time of each satellite changes over the years of operation, leading to shifts of over five hours.As lake surface water temperature is changing throughout a day, an overpass time at 1 p.m. would lead to higher average temperatures than an overpass time at 8 am.Therefore, drifting orbit leads to an artificial trend in long-time records when looking at only one platform.When it comes to time series over several platforms, as is the case here, the impact is hard to predict.Some authors claimed that the impact is negligible (e.g., [32]); others correct for the DTC (e.g., [38]).
For this study, we used a DTC correction based on the DTC model described in Equations ( 2)-(5) after Duan et al. [39].The model was first developed by Inamdar et al. [40] for land surface temperature (LST) and tested against in-situ data with various land cover types, including water (Lake Tahoe).Duan et al. [39] improved the model by deducing the width ω (Equation ( 4)) over the half-period of the cosine term from the thermal diffusion equation.This model was successfully implemented and used for daytime LSWT retrieval by Pareeth et al. (e.g., [38]).

Diurnal Temperature Cycle (DTC) Correction
The NOAA satellites are subject to relatively strong orbit drifts [36,37].This means that through the years of operation, the overpass time at the same location is shifted.Figure 4 shows all scenes and recording times used during this study.One can see that the overpass time of each satellite changes over the years of operation, leading to shifts of over five hours.As lake surface water temperature is changing throughout a day, an overpass time at 1 p.m. would lead to higher average temperatures than an overpass time at 8 am.Therefore, drifting orbit leads to an artificial trend in long-time records when looking at only one platform.When it comes to time series over several platforms, as is the case here, the impact is hard to predict.Some authors claimed that the impact is negligible (e.g., [32]); others correct for the DTC (e.g., [38]).
For this study, we used a DTC correction based on the DTC model described in Equations ( 2)-( 5) after Duan et al. [39].The model was first developed by Inamdar et al. [40] for land surface temperature (LST) and tested against in-situ data with various land cover types, including water (Lake Tahoe).Duan et al. [39] improved the model by deducing the width ω (Equation ( 4)) over the half-period of the cosine term from the thermal diffusion equation.This model was successfully implemented and used for daytime LSWT retrieval by Pareeth et al. (e.g., [38]). with T d s and T n s are the day and night time surface temperatures, t is the time, T 0 the residual temperature at sunrise, T a is the temperature amplitude, δT is the difference in temperature between T 0 and T (t→∞) , ω is the width over the half-period of the cosine term, t sr is the sunrise time, t m is the time where daily maximum temperature is reached, t s is the time where free attenuation starts, and k is the attenuation constant.
During the here-presented study, instead of adjusting the satellite derived LSWT to a temperature equivalent of a standard acquisition time (12:00 UTC), we adjusted the temperature to the daily average temperature, computed from the modeled DTC curves.For each lake and for each month, typical DTC anomaly curves were computed by fitting the model to the available satellite LSWT data.The final LSWT was then adjusted to the daily average temperature by subtracting the temperature anomaly.The use of this approach improved the validation results within the here-presented study (see Section 4.2).

Trend Analysis
The linear trends were computed using the homogenized time series at the 26 locations (Table 1) described in the previous chapter.Non-parametric trend analysis was performed using the Thiel-Sen slope to quantify the trend and the Mann-Kendal test to determine its significance.
Several studies do compute only summer trends (e.g., [17]).This has the advantage of not having to deal with ice onset and breakup cycle.However, Winslow et al. [41] showed that seasonal trends such as a summer trend do not well represent the effects of climate change on the lakes.During this study, trends for all four seasons, as well as annual trends, were computed.When it comes to the definition of the seasons, there are mainly two ways described in literature: (1) summer defined as June, July, August (JJA), e.g., [42][43][44]; (2) summer defined as July, August, September (JAS), e.g., [17,32,41,45].However, in the case of the lakes studied here, the influence of spring warming on the JJA period and the influence of fall cooling on the JAS period are disproportionately high (see Figure 5).We concluded that both JAS and JJA were not suitable for "summer only trend" in a restricted sense.For this study, we overcame this issue by taking the middle way; thus, the periods start and end in the middle of the months: Summer is defined as the period from 15 June to 14 September (JJAS), autumn as from 15 September to 14 December (SOND), winter from 15 December to 14 March (DJFM), and, finally, spring from 15 March to 14 June (MAMJ).Choosing these periods, the characteristics of the seasons are more pronounced (spring = warming, fall = cooling, summer = hot, and winter = cold) looking at the annual temperature cycle (Figure 5).

Validation of Retrieved Lake Surface Water Temperature
The main performance indicator for the retrieval validation during this study is the standard deviation (SD), whereas the root mean square error (RMSE) and the bias are considered as secondary performance indicators.As we are primarily interested in temperature patterns such as temporal trends or spatial gradients, the uncertainty between two distinct measurements in time or space (SD) needs to be as small as possible, whereas the systematic part of the uncertainty (bias) is of lower

Validation of Retrieved Lake Surface Water Temperature
The main performance indicator for the retrieval validation during this study is the standard deviation (SD), whereas the root mean square error (RMSE) and the bias are considered as secondary performance indicators.As we are primarily interested in temperature patterns such as temporal trends or spatial gradients, the uncertainty between two distinct measurements in time or space (SD) needs to be as small as possible, whereas the systematic part of the uncertainty (bias) is of lower importance.In fact, the bias affects all measurements in the same way, and therefore does not mask patterns such as spatial gradients and temporal trends.For the LSWT retrieval, the variable part of the uncertainty (SD) is mainly induced by variability in atmospheric composition, weather conditions, small-scale lake dynamic processes, and sensor noise, whereas the systematic part of the uncertainty (bias) is induced by factors such as skin-bulk temperature difference, in-situ sensor location, or sensor-related issues.
Further, it must be mentioned that during validation, we are comparing two measurements of a physical quantity against each other.This means that none of the two measurements (in-situ and satellite) do represent the truth to the full extent.This got even worse when considering the fact that they technically do not measure exactly the same quantity: The in-situ temperature measurement deduces the water temperature at a specific point in space (location and depth), whereas the satellite measurement estimates the temperature of the uppermost skin layer, averaged over a certain spatial extent, which is defined by its resolution.The exact amount of uncertainty issued from the described difference in measurement technique and measured quantity could not be determined during this study.Also, the correction methods for the skin to bulk difference as proposed by Donlon [46] or Minnet [47] for SST retrieval did not improve the validation performance: The method, applied to LSWT, only affects the systematic bias and has no effect on the variance of the error.
In Figure 6, one can see the difference when validating the same satellite data with data from two different in-situ locations close to each other.This relativizes the significance of such validation results, even though they still give a general indication about how well the retrieval works.However, it reminds us to be careful with the in-situ data, and to not take it as the finite truth.In Figure 6, one can see the difference when validating the same satellite data with data from two different in-situ locations close to each other.This relativizes the significance of such validation results, even though they still give a general indication about how well the retrieval works.However, it reminds us to be careful with the in-situ data, and to not take it as the finite truth.Both show a good performance and a high number of matchups despite the in-situ locations close to the shore.However, the correlation with the satellite data is significantly higher in Bregenz (a) than in Lindau (b).For both validations, the same satellite data is used.
Table 4 shows the validation results at all in-situ stations.There seem to be two main factors influencing the validation performance.First, there is the latitude, as it seems that with higher cloud cover rate and longer periods of ice cover during the cold season, also the quality of the retrieval decreases.On the other hand, the data type of the in-situ data seems to play an important role too: There are three types of in-situ data: daily averaged, hourly averaged, and manually taken instant measurements.The validation against instant and hourly averaged measurements yields the best a b Figure 6.The validation with hourly averaged measurements at Lake Constance for the station in Bregenz (a) and for the station in Lindau (b).Both show a good performance and a high number of matchups despite the in-situ locations close to the shore.However, the correlation with the satellite data is significantly higher in Bregenz (a) than in Lindau (b).For both validations, the same satellite data is used.
Table 4 shows the validation results at all in-situ stations.There seem to be two main factors influencing the validation performance.First, there is the latitude, as it seems that with higher cloud cover rate and longer periods of ice cover during the cold season, also the quality of the retrieval decreases.On the other hand, the data type of the in-situ data seems to play an important role too: There are three types of in-situ data: daily averaged, hourly averaged, and manually taken instant measurements.The validation against instant and hourly averaged measurements yields the best results, as it allows for close temporal matchups between satellite overpass and in-situ measurements (matchup window ±1 h).

Table 4.
Validation results for each in-situ station.The data is ordered by latitude from south to north.The measurement type describes the averaging of the in-situ measurement: daily averaged, hourly averaged, and instant measurements.Instant measurements are single measurements that were manually taken during measurement campaigns.As the exact timestamps for the instant measurements where not known at Lake Balaton (Kesztehly old) and Lake Trichonida, 12:00 local time is assumed.Amongst the instant measurements, the validation at Müritz yields the best results (SD ± 0.82 K, Figure 7a), but due to its sporadic sampling frequency only 66 matchups where available.The performance of the instant measurements at Lake Balaton (SD ± 1.44 K, Figure 7b) and Lake Trichonidas (SD ± 1.04 K) is lower, but at these locations the exact in-situ measurement timestamp was unknown.The validation with hourly averaged measurements has considerably higher number of matchups, and the performances reach from ±0.95 K (Lake Constance, Bregenz new, Figure 6a) to ±1.40 K (Lake Constance, Lindau, Figure 6b).The validation with daily averaged in-situ data shows, generally, a lower performance than with instant and hourly measurements.The SD reaches from very low ±1.07K (Geneva, INRA Figure 8a) to rather increased values of ±2.05 K at Lake Inarijärvi (Figure 8b).

Location
was unknown.The validation with hourly averaged measurements has considerably higher number of matchups, and the performances reach from ±0.95 K (Lake Constance, Bregenz new, Figure 6a) to ±1.40 K (Lake Constance, Lindau, Figure 6b).The validation with daily averaged in-situ data shows, generally, a lower performance than with instant and hourly measurements.The SD reaches from very low ±1.07K (Geneva, INRA Figure 8a) to rather increased values of ±2.05 K at Lake Inarijärvi (Figure 8b).The validation at Lake Geneva INRA shows good performance and a high number of matchups through the whole year.At Inarijärvi, however, the weakest validation performance of all is found.Inarijärvi is located north of the polar cycle and has a long period of ice cover per year.Further, the in-situ measurement station is located within a very complex shore geometry with many small islands.

Validation of the DTC Correction and Trends
For the validation with daily in-situ data, the available satellite LSWT is compiled to daily average temperatures.To reduce the impact of the daily temperature cycle (DTC), a DTC correction is applied to the satellite retrievals.The validation with DTC-corrected data results in a decrease of SD by 0% to 15.6% depending on the lakes, and in an average decrease of 6% or 0.11 K, respectively (Table 5).
The trends from the in-situ stations with a period covered of at least 25 years were compared to the trends computed from LSWT retrievals.This was done for both the DTC-corrected satellite data and for the non-adjusted satellite data.Table 6 shows results from this comparison.At three of the six locations, the DTC adjusted trends are considerably closer to the in-situ trends (Inarijärvi, Lappajärvi, and Oulujärvi).At two locations (Pyhäselka and Peipsi), the DTC adjustment has almost The validation at Lake Geneva INRA shows good performance and a high number of matchups through the whole year.At Inarijärvi, however, the weakest validation performance of all is found.Inarijärvi is located north of the polar cycle and has a long period of ice cover per year.Further, the in-situ measurement station is located within a very complex shore geometry with many small islands.

Validation of the DTC Correction and Trends
For the validation with daily in-situ data, the available satellite LSWT is compiled to daily average temperatures.To reduce the impact of the daily temperature cycle (DTC), a DTC correction is applied to the satellite retrievals.The validation with DTC-corrected data results in a decrease of SD by 0% to 15.6% depending on the lakes, and in an average decrease of 6% or 0.11 K, respectively (Table 5).
The trends from the in-situ stations with a period covered of at least 25 years were compared to the trends computed from LSWT retrievals.This was done for both the DTC-corrected satellite data Table 7. Yearly and seasonal trends at all studied locations ordered from South to North.The covered period is slightly below the maximum of 36 years, depending on availability of cloud-free scenes at each location.The grayed-out values are trends (K/dec) in which Kendall p-value > 0.They are not retained to compute the average trends showed on the last line.The strongest seasonal trend of each lake is marked with an asterisk.Blank values signify that there were not enough values to compute trends, mostly due to ice cover during the cold period.All annual trends where significant and positive, thus indicating a warming of the lakes water temperature (Table 7).The strongest annual warming trend is observed at Lake Onega with 0.76 K/dec, whereas the weakest trend is to be found at Djorf Torba Reservoir with 0.17 K/dec.The trend at Lake Onega is remarkable, as there is no similar trend within the studied lakes; the next highest trend is 0.54 K/dec for lake Pyhäjärvi.The trend at Djorf Torba has to be looked at with care, as there has been a strong siltation going on since the drought period in the 1990s: The volume loss was estimated at 35% in 2014, with the siltation still going on at a rate of 0.9% per year [48].
Figure 9 shows details of the temperature anomaly in annual resolution of lakes representing geographical extremes (N-S, W-E).The significance of a long time series including years from 1981 onwards is clearly shown with increased frequency of blue bars during the earlier years (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995) in comparison to the latest years with higher temperatures.Especially the warm period of the last five years is obvious, independently of the geographical location of the lakes.However, these high temperatures are more pronounced in northern and eastern areas, whereas the lake in Iceland (Thingvallavatn) shows reduced temperature increase, probably due to the oceanic influence.Also, one of the southern lakes (Trichonida, Figure 9) shows a reduced increase of temperature, which is probably also influenced by the proximity to the Mediterranean Sea.In all studied lakes, we can identify spatial gradients of the magnitude of temperature trends: The trends are increasing from south (Lake Trichonida 0.21 K/decade, Figure 9a) to north (e.g., Inarijärvi 0.50 K/decade, Figure 9b), and also from west (e.g., Thingvallavatn 0.38 K/decade, Figure 9c) to east (e.g., Lake Onega 0.76 K/decade, Figure 9d).The gradient is strongest on the diagonal from the south-west to north-east, whereas the magnitude of the trends shows a significant correlation with the latitude (R 2 = 58%, Figure 10a) and a rather weak correlation with longitude (R 2 = 24%, Figure 10b).When looking at the strongest seasonal trends in Table 7, one can see that at lower latitudes the winter trends and some spring trends are dominant.This is typically the case for the subalpine lakes such as Lake Geneva, Lake Garda, and Lake Constance, which all show dominant winter trends.When moving north, first the spring (e.g., Balaton, Müritz, and Vättern) and then the summer trends (e.g., Pyhäjärvi, Päijänne, and Onega) start to dominate.Finally, for the highest latitudes, the autumn trends are clearly dominating (all lakes north to Lake Onega).
Looking at the north-south gradient of the seasonal trends (Figure 11), we see that winter trends seem almost unaffected by the latitude in contrast to the other seasonal trends.Especially, autumn trends are very sensible to latitude: At low latitudes, temperature trends in autumn are less pronounced than winter or spring trends, whereas at northern latitudes they are the strongest seasonal trends.The West-East gradient of the seasonal trends is less significant, and probably influenced by the transition from oceanic to continental climate.However, this could not be investigated in more detail with the available data.When looking at the strongest seasonal trends in Table 7, one can see that at lower latitudes the winter trends and some spring trends are dominant.This is typically the case for the subalpine lakes such as Lake Geneva, Lake Garda, and Lake Constance, which all show dominant winter trends.When moving north, first the spring (e.g., Balaton, Müritz, and Vättern) and then the summer trends (e.g., Pyhäjärvi, Päijänne, and Onega) start to dominate.Finally, for the highest latitudes, the autumn trends are clearly dominating (all lakes north to Lake Onega).
Looking at the north-south gradient of the seasonal trends (Figure 11), we see that winter trends seem almost unaffected by the latitude in contrast to the other seasonal trends.Especially, autumn trends are very sensible to latitude: At low latitudes, temperature trends in autumn are less pronounced than winter or spring trends, whereas at northern latitudes they are the strongest seasonal trends.When looking at the strongest seasonal trends in Table 7, one can see that at lower latitudes the winter trends and some spring trends are dominant.This is typically the case for the subalpine lakes such as Lake Geneva, Lake Garda, and Lake Constance, which all show dominant winter trends.When moving north, first the spring (e.g., Balaton, Müritz, and Vättern) and then the summer trends (e.g., Pyhäjärvi, Päijänne, and Onega) start to dominate.Finally, for the highest latitudes, the autumn trends are clearly dominating (all lakes north to Lake Onega).
Looking at the north-south gradient of the seasonal trends (Figure 11), we see that winter trends seem almost unaffected by the latitude in contrast to the other seasonal trends.Especially, autumn trends are very sensible to latitude: At low latitudes, temperature trends in autumn are less pronounced than winter or spring trends, whereas at northern latitudes they are the strongest seasonal trends.The West-East gradient of the seasonal trends is less significant, and probably influenced by the transition from oceanic to continental climate.However, this could not be investigated in more detail with the available data.The West-East gradient of the seasonal trends is less significant, and probably influenced by the transition from oceanic to continental climate.However, this could not be investigated in more detail with the available data.

Discussion
The validation shows that with careful processing, the difficulties related to the AVHRR sensor such as orbit drift and small georegistration errors and others could be overcome.The importance of the in-situ measurements has been emphasized during the validation, and it was shown that the quality of the in-situ measurement is significantly influencing the performance of the validation.This indicates that the satellite retrieval has reached an accuracy level at which more accurate in-situ measurements are needed to further assess the quality of the LSWT retrieval.Also, the derived trends could be well validated with the in-situ measurements.
Concerning the significance test of the trends, only significant trends (Kendal p-value < 0.05 Table 7) were retained for further analysis.In contrast to the yearly trends, not all seasonal trends fulfilled this criterion.In the case in which the Kendall p-value is above 0.05, there are two plausible interpretations: (1) There is a strong trend, but due to elevated variance in the data and/or due to short periods and low amount of data, the trend cannot be determined with confidence.With this, erroneous estimates of trends can be detected (e.g., winter trend at Lake Päijiänne Table 7); (2) There is no trend or a weak trend in the data, and the p-value is increased, because we are close to the null-hypothesis.In this case, we eliminate a potentially correct value of weak trend close to 0 K/decades such as winter and autumn trends from lake Sanguinet (Figure 12).Thus, it is important to emphasize that with the used method, we can detect significant trends, but we cannot detect trends near 0 K/decade with certitude.This explains the discrepancy when comparing the annual trend with the valid seasonal trends from Table 7, in which our intuition is that the average of the seasonal trends is more or less equal to the yearly trend.

Discussion
The validation shows that with careful processing, the difficulties related to the AVHRR sensor such as orbit drift and small georegistration errors and others could be overcome.The importance of the in-situ measurements has been emphasized during the validation, and it was shown that the quality of the in-situ measurement is significantly influencing the performance of the validation.This indicates that the satellite retrieval has reached an accuracy level at which more accurate in-situ measurements are needed to further assess the quality of the LSWT retrieval.Also, the derived trends could be well validated with the in-situ measurements.
Concerning the significance test of the trends, only significant trends (Kendal p-value < 0.05 Table 7) were retained for further analysis.In contrast to the yearly trends, not all seasonal trends fulfilled this criterion.In the case in which the Kendall p-value is above 0.05, there are two plausible interpretations: (1) There is a strong trend, but due to elevated variance in the data and/or due to short periods and low amount of data, the trend cannot be determined with confidence.With this, erroneous estimates of trends can be detected (e.g., winter trend at Lake Päijiänne Table 7); (2) There is no trend or a weak trend in the data, and the p-value is increased, because we are close to the nullhypothesis.In this case, we eliminate a potentially correct value of weak trend close to 0 K/decades such as winter and autumn trends from lake Sanguinet (Figure 12).Thus, it is important to emphasize that with the used method, we can detect significant trends, but we cannot detect trends near 0 K/decade with certitude.This explains the discrepancy when comparing the annual trend with the valid seasonal trends from Table 7, in which our intuition is that the average of the seasonal trends is more or less equal to the yearly trend.Trends are sensitive to different factors such as the period in which a trend is determined, the homogenization of the data, and the type of trend (seasonal, yearly, etc.).Therefore, the comparison with trends from individual lakes in other studies is not meaningful.However, on the other hand, more general statements can be looked at in detail.For example, Schneider and Hook [17] mention elevated summer trends in Northern Europe of up to around 0.8 K/dec.Our results suggest the same tendency of stronger trends in north-east Europe, whereas the strongest trends are rather to be found in autumn than in summer.The approximatively 0.8 K/dec in summer is only reached in Lake Onega, whereas at Inarijärvi and at the Lokka Reservoir such elevated trends are occurring during autumn.
Another observation from Schneider and Hook [17], that lakes in Northern Europe are warming more rapidly than local air temperature, was confirmed (see e.g., Lake Vättern Figure 13a).According to our results, this phenomenon is not limited to northern Europe: it can also be observed at some lakes in central Europe such as Lake Constance (Figure 13b).However, the difference in trends between local air temperature and LSWT is reduced when moving south towards warmer lakes.Trends are sensitive to different factors such as the period in which a trend is determined, the homogenization of the data, and the type of trend (seasonal, yearly, etc.).Therefore, the comparison with trends from individual lakes in other studies is not meaningful.However, on the other hand, more general statements can be looked at in detail.For example, Schneider and Hook [17] mention elevated summer trends in Northern Europe of up to around 0.8 K/dec.Our results suggest the same tendency of stronger trends in north-east Europe, whereas the strongest trends are rather to be found in autumn than in summer.The approximatively 0.8 K/dec in summer is only reached in Lake Onega, whereas at Inarijärvi and at the Lokka Reservoir such elevated trends are occurring during autumn.
Another observation from Schneider and Hook [17], that lakes in Northern Europe are warming more rapidly than local air temperature, was confirmed (see e.g., Lake Vättern Figure 13a).According to our results, this phenomenon is not limited to northern Europe: it can also be observed at some lakes in central Europe such as Lake Constance (Figure 13b).However, the difference in trends between local air temperature and LSWT is reduced when moving south towards warmer lakes.Additionally, for Lake Bolsena (Figure 13c), the air temperature trend tends even to be more pronounced than the LSWT trends.However, the reasons for the strong LSWT trends compared to local air temperature were not further investigated here.
The seasonal trends we detected show a high intra-seasonal variability within the same lake (Table 7), and therefore we can state that none of the seasons are representative for the overall trend within a lake.This is perfectly aligned with the main statement of Winslow et al. [41], who claim that summer trends do not fully represent the effects of climate change on lake temperatures.When we compare the yearly temperature cycle deduced from the first and the last 15-year period of the dataset, we can observe that the temperature evolution is very heterogeneous over the year, and also very different from lake to lake.For example, at Lake Bolsena and Lake Constance (Figure 14a,b), we observe significant warming during winter and spring.However, then there are also some short periods in which the recent years were colder than the first years of the dataset (mid-summer at Lake Constance, end of summer and start of autumn at Lake Bolsena).In contrast to this, at Lake Vättern and Lake Onega (Figure 14c,d), a moderate to significant warming during the whole year is observed, whereas this is more pronounced in summer and autumn.Lake Onega shows a significantly earlier and more pronounced warming in spring and during summer, resulting in a remarkable summer temperature peak that is higher and earlier than in the first years of the dataset.Additionally, for Lake Bolsena (Figure 13c), the air temperature trend tends even to be more pronounced than the LSWT trends.However, the reasons for the strong LSWT trends compared to local air temperature were not further investigated here.
The seasonal trends we detected show a high intra-seasonal variability within the same lake (Table 7), and therefore we can state that none of the seasons are representative for the overall trend within a lake.This is perfectly aligned with the main statement of Winslow et al. [41], who claim that summer trends do not fully represent the effects of climate change on lake temperatures.When we compare the yearly temperature cycle deduced from the first and the last 15-year period of the dataset, we can observe that the temperature evolution is very heterogeneous over the year, and also very different from lake to lake.For example, at Lake Bolsena and Lake Constance (Figure 14a,b), we observe significant warming during winter and spring.However, then there are also some short periods in which the recent years were colder than the first years of the dataset (mid-summer at Lake Constance, end of summer and start of autumn at Lake Bolsena).In contrast to this, at Lake Vättern and Lake Onega (Figure 14c,d), a moderate to significant warming during the whole year is observed, whereas this is more pronounced in summer and autumn.Lake Onega shows a significantly earlier and more pronounced warming in spring and during summer, resulting in a remarkable summer temperature peak that is higher and earlier than in the first years of the dataset.Further, also within the same water body, trends can change significantly.This result is in agreement with the results of [16], in which intra-lake heterogeneity of trends was studied in large lakes.At lake Peipus (Figure 15), the yearly trend in the open waters in the north indicates a warming of 0.24 K/decade, whereas the warming rate almost doubles (0.44 K/decade) within the narrower part of the lake around Mehikoorma, in which also an in-situ station is located.This shows that the location where the temperature is sampled has an impact on the resulting trends.However, the example at lake Peipus shows us also that trends computed with the data of an in-situ measurement taken at a specific location are not necessarily representative of the temperature trend of the entire lake.Further, also within the same water body, trends can change significantly.This result is in agreement with the results of [16], in which intra-lake heterogeneity of trends was studied in large lakes.At lake Peipus (Figure 15), the yearly trend in the open waters in the north indicates a warming of 0.24 K/decade, whereas the warming rate almost doubles (0.44 K/decade) within the narrower part of the lake around Mehikoorma, in which also an in-situ station is located.This shows that the location where the temperature is sampled has an impact on the resulting trends.However, the example at lake Peipus shows us also that trends computed with the data of an in-situ measurement taken at a specific location are not necessarily representative of the temperature trend of the entire lake.2) with the trend in the open waters of the northern part (see coordinates in Table 1).

Conclusions
During this study, we could show that the AVHRR LAC data (1 km) archive at University of Bern can be used to yield good quality LSWT of an exceptional long-time span of over 35 years.This offers the possibility to study lake dynamics related to climate change.The quality of the retrieved LSWT could be confirmed by validation considering in-situ measurements, but it could also be shown that the quality of the validation data is essential.
A homogeneous LSWT dataset from 26 European lakes covering the period from 1981 to 2016 was created.With this, we were able to compute LSWT trends and to compare average yearly temperature cycles of different periods, of different geographic locations.
For all locations, significant LSWT trends indicating a warming were found.A gradient of increasing warming trends was revealed when moving from the south-west to the north-east of Europe.Whereas winter trends are not significantly affected, autumn trends are very sensitive to the location of the lake, especially to the latitude.Looking at the dominant seasonal trends per lake, we observe a shift from stronger winter trends in the south to dominant autumn trends in northern lakes, passing by a mixture of strong spring and summer trends in-between.
Generally, the strong interannual heterogeneity of the trends was detected in all lakes, showing that using only a defined period of the year to determine the temperature trend is not representative.It could be confirmed that the different seasons contribute unequal to the yearly trends, and that trend analysis based on summer seasons only does not show the complete reality.
Also, intra-lake heterogeneity of trends was detected.The location of measurements influences the trend directly, implying that a one-to-one comparison of trends between in-situ and satellite data has to be handled with care.We could show that temperature trends computed with data from insitu measurements from individual stations are not suited to represent trends of an entire lake.This emphasizes the importance of satellite-based lake monitoring, as it adds significantly to the understanding of the big picture.
The comparison with local air temperature trends showed that many lakes have a stronger warming rate than local air temperature.This effect was shown to be more significant in colder lakes and in higher latitudes, e.g., Lake Vättern shows a stronger increase of LSWT than air temperature, but for Lake Bolsena in Italy the LSWT trend is even slightly weaker than for local air temperature.2) with the trend in the open waters of the northern part (see coordinates in Table 1).

Conclusions
During this study, we could show that the AVHRR LAC data (1 km) archive at University of Bern can be used to yield good quality LSWT of an exceptional long-time span of over 35 years.This offers the possibility to study lake dynamics related to climate change.The quality of the retrieved LSWT could be confirmed by validation considering in-situ measurements, but it could also be shown that the quality of the validation data is essential.
A homogeneous LSWT dataset from 26 European lakes covering the period from 1981 to 2016 was created.With this, we were able to compute LSWT trends and to compare average yearly temperature cycles of different periods, of different geographic locations.
For all locations, significant LSWT trends indicating a warming were found.A gradient of increasing warming trends was revealed when moving from the south-west to the north-east of Europe.Whereas winter trends are not significantly affected, autumn trends are very sensitive to the location of the lake, especially to the latitude.Looking at the dominant seasonal trends per lake, we observe a shift from stronger winter trends in the south to dominant autumn trends in northern lakes, passing by a mixture of strong spring and summer trends in-between.
Generally, the strong interannual heterogeneity of the trends was detected in all lakes, showing that using only a defined period of the year to determine the temperature trend is not representative.It could be confirmed that the different seasons contribute unequal to the yearly trends, and that trend analysis based on summer seasons only does not show the complete reality.
Also, intra-lake heterogeneity of trends was detected.The location of measurements influences the trend directly, implying that a one-to-one comparison of trends between in-situ and satellite data has to be handled with care.We could show that temperature trends computed with data from in-situ measurements from individual stations are not suited to represent trends of an entire lake.This emphasizes the importance of satellite-based lake monitoring, as it adds significantly to the understanding of the big picture.
The comparison with local air temperature trends showed that many lakes have a stronger warming rate than local air temperature.This effect was shown to be more significant in colder lakes and in higher latitudes, e.g., Lake Vättern shows a stronger increase of LSWT than air temperature, but for Lake Bolsena in Italy the LSWT trend is even slightly weaker than for local air temperature.
The dataset generated within this study will be further investigated, and it has the potential to reveal further findings about the effects of climate change on European lakes.An extension of the dataset to more lakes in Europe would also be of great use for further investigations, especially regarding the east-west gradient of temperature trends, for which evidence was found.Also, a closer look at the influence on temperature trends of morphologic parameters such as the average depth of lakes could be the topic of further investigations.
For an extension of the dataset to other lakes in Europe, the size of the lakes is the only limiting factor, whereas for an extension on a global scale, the availability of AVHRR LAC data would be limiting.The application of the method to more recent sensors with higher resolution would be interesting for lake dynamic studies, but it is limited to a more recent time-period due to the availability of adequate satellite data.
Using open street map data for the of the LWM (land water mask) worked well.However, in the case of a larger number of lakes, and when taking in account the changes in time of global surface water occurrence, a global water mapping dataset such as presented in Pekel et al. 2016 [49] could be a more robust alternative for the generation of the LWM.

Figure 1 .
Figure 1.Advanced Very High Resolution Radiometer (AVHRR) scenes subset available within the University of Bern data archive with the lakes covered by the present study.Blue dots are all lakes for which Lake Surface Water Temperature (LSWT) data was retrieved and analysed; the green dots are the location where validation with in-situ data was performed.

Figure 1 .
Figure 1.Advanced Very High Resolution Radiometer (AVHRR) scenes subset available within the University of Bern data archive with the lakes covered by the present study.Blue dots are all lakes for which Lake Surface Water Temperature (LSWT) data was retrieved and analysed; the green dots are the location where validation with in-situ data was performed.

Figure 2 .
Figure2.Availability of AVHRR data by platform and by years.In blue, the data that was used for the present study, and, in red, the data that was not processed due to the elevated amount of inconsistent data.

Figure 2 .
Figure2.Availability of AVHRR data by platform and by years.In blue, the data that was used for the present study, and, in red, the data that was not processed due to the elevated amount of inconsistent data.

Figure 3 .
Figure 3. Schematic representation of the LSWT processing as used for this work.

Figure 3 .
Figure 3. Schematic representation of the LSWT processing as used for this work.

Figure 4 .
Figure 4. Overpass times of all AVHRR scenes in which Lake Geneva is visible.Each dot represents a scene; the different colors represents different platforms.

Figure 4 .
Figure 4. Overpass times of all AVHRR scenes in which Lake Geneva is visible.Each dot represents a scene; the different colors represents different platforms.

Figure 5 .
Figure 5.Comparison for different season periods (JJA, JJAS, JAS) in the yearly temperature cycle at lake Constance, Müritz, and Lake Inari.The middle row, representing the seasonal periods as used during this study is well adapted to the actual annual cycles, as the summer temperature peak is in the center of the period.

Figure 5 .
Figure 5.Comparison for different season periods (JJA, JJAS, JAS) in the yearly temperature cycle at lake Constance, Müritz, and Lake Inari.The middle row, representing the seasonal periods as used during this study is well adapted to the actual annual cycles, as the summer temperature peak is in the center of the period.

Figure 6 .
Figure6.The validation with hourly averaged measurements at Lake Constance for the station in Bregenz (a) and for the station in Lindau (b).Both show a good performance and a high number of matchups despite the in-situ locations close to the shore.However, the correlation with the satellite data is significantly higher in Bregenz (a) than in Lindau (b).For both validations, the same satellite data is used.

Figure 7 .Figure 7 .
Figure 7.The validation with sporadic instant measurement at Müritz (a) and at Lake Balaton (b).Both in-situ measurements where manually taken in open waters.The bi-monthly sampling frequency only allows for few matchups per year.Nevertheless, the measurements at Müritz yield the best performance thanks to close temporal matchups and the open waters.At Lake Balaton however, the exact timestamp was unknown (only day).Therefore, 12:00 local time was assumed, which led to less accurate temporal matchups and lower correlation.

Figure 8 .
Figure 8.The validation with daily in-situ measurement at Lake Geneva, INRA (a) and Inarijärvi (b).The validation at Lake Geneva INRA shows good performance and a high number of matchups through the whole year.At Inarijärvi, however, the weakest validation performance of all is found.Inarijärvi is located north of the polar cycle and has a long period of ice cover per year.Further, the in-situ measurement station is located within a very complex shore geometry with many small islands.

Figure 8 .
Figure 8.The validation with daily in-situ measurement at Lake Geneva, INRA (a) and Inarijärvi (b).The validation at Lake Geneva INRA shows good performance and a high number of matchups through the whole year.At Inarijärvi, however, the weakest validation performance of all is found.Inarijärvi is located north of the polar cycle and has a long period of ice cover per year.Further, the in-situ measurement station is located within a very complex shore geometry with many small islands.

Figure 9 .Figure 9 .
Figure 9. Temperature trends of Lake Trichonidas (a), Inarijärvi (b), Thingvallavatn (c), and Lake Onega (d).Blue bars indicate temperatures below the long-term mean, whereas the red bars show warmer years.The dashed line shows the linear trend, whereas the solid black line shows the locally weighted scatterplot smoothing (LOWESS) over a 2-years period.Lake Trichonidas (a) is the lake with the lowest latitude (except for Djorf Torba), whereas Inarivjärvi (b) is the lake at highest latitude.The same applies for the longitude where it goes from Thingvallavatn (c) in the west to Lake Onega (d) in the east.

Figure 10 .
Figure 10.South-north gradient (a) and the west-east gradient (b) of the yearly LSWT trend.The red dashed line is a simple linear model derived with least square regression, and R2 is the coefficient of determination.

Figure 11 .
Figure 11.South-North gradient of seasonal trends from European lakes.The dashed lines are simple linear models derived with least square regression of the data points.The reservoirs lakes are excluded from the graph.

Figure 10 .
Figure 10.South-north gradient (a) and the west-east gradient (b) of the yearly LSWT trend.The red dashed line is a simple linear model derived with least square regression, and R2 is the coefficient of determination.

25 Figure 10 .
Figure 10.South-north gradient (a) and the west-east gradient (b) of the yearly LSWT trend.The red dashed line is a simple linear model derived with least square regression, and R2 is the coefficient of determination.

Figure 11 .
Figure 11.South-North gradient of seasonal trends from European lakes.The dashed lines are simple linear models derived with least square regression of the data points.The reservoirs lakes are excluded from the graph.

Figure 11 .
Figure 11.South-North gradient of seasonal trends from European lakes.The dashed lines are simple linear models derived with least square regression of the data points.The reservoirs lakes are excluded from the graph.

Figure 14 .
Figure 14.Comparison of the Yearly Temperature Cycle from two different 15-years-periods at Lake Bolsena (a), Lake Constance (b), Lake Vättern (c), and Lake Onega (d).The red and the blue lines represent the yearly cycle of the two periods; the black line is the difference between the two periods.

Figure 14 .
Figure 14.Comparison of the Yearly Temperature Cycle from two different 15-years-periods at Lake Bolsena (a), Lake Constance (b), Lake Vättern (c), and Lake Onega (d).The red and the blue lines represent the yearly cycle of the two periods; the black line is the difference between the two periods.

Figure 15 .
Figure 15.Significant difference in trends within Lake Peipus when comparing the trend near the Mehikoorma in-situ station (see coordinates in Table2) with the trend in the open waters of the northern part (see coordinates in Table1).

Figure 15 .
Figure 15.Significant difference in trends within Lake Peipus when comparing the trend near the Mehikoorma in-situ station (see coordinates in Table2) with the trend in the open waters of the northern part (see coordinates in Table1).

Table 1 .
List of Lakes where LSWT was retrieved and trend analysis performed.The acronyms are the same as those used in Figure1.

Table 3 .
The six quality levels indicate which additional tests are passed.The levels are cumulative (e.g., Q4 means Q3 AND VZA < 55).The average data availability (compared to Q1) and the average standard deviation during the validation are shown for the different quality levels.
* Cloud and land mask applied, and physical inconsistency removed; ** passed if the detection of the feature is confirmed (additional georeferentiation test).