The Importance of Systematic Spatial Variability in the Surface Heat Flux of a Large Lake: A Multiannual Analysis for Lake Geneva

The spatiotemporal surface heat flux (SurHF) distribution over Lake Geneva, the largest lake in Western Europe, was estimated for a 7‐year period (2008–2014). Data sources included hourly maps of over‐the‐lake assimilated meteorological data from a validated numerical weather model and lake surface water temperature (LSWT) from satellite imagery. A set of bulk algorithms, previously optimized and calibrated at two locations in Lake Geneva, was used. Results indicate a systematic long‐term average spatial range of >40 Wm‐2 between different parts of the lake and little year‐to‐year variability. This variability is mainly due to topographically induced wind sheltering over parts of the lake, which in turn produces spatial variability in the sensible and latent heat fluxes. These results are supported by a systematic spatial heat content variability obtained from long‐term temperature profile measurements in the lake. During spring, a lower SurHF spatial range was evident. Unlike other seasons, the spring spatial variability of air‐water temperature differences and, to a lesser extent, the global radiation variability resulting from sheltering by the mountainous topography were the main drivers of the SurHF spatial variability. Analysis of the atmospheric thermal boundary layer showed stable conditions from March to early June and unstable conditions for the rest of the year. This regime change can explain the low SurHF spatial variability observed during spring. The results emphasize that spatial variability in meteorological and LSWT patterns, and consequently in the spatiotemporal SurHF data, should be considered when assessing the time evolution of the heat budget of large lakes.


Introduction
Wind forcing and surface heat flux (SurHF) are the main external parameters driving lake thermodynamics and hydrodynamics, which in turn influence the chemical and biological properties of aquatic systems (e.g., Bonvin et al., 2013;Churchill & Kerfoot, 2007;Finlay et al., 2015) and ecosystem dynamics (e.g., Bauersachs et al., 2015;Beaulieu et al., 2013;Binding et al., 2018). SurHF affects the total lake heat content and thus water temperatures, which are often taken as a climate change indicator (e.g., Adrian et al., 2009). An increase in lake surface water temperature (LSWT) over the past decades was observed (O'Reilly et al., 2015) and its spatial heterogeneity within lakes was reported as well (Woolway & Merchant, 2018). However, Ye et al. (2019) recently demonstrated that the lake heat content, which is mainly determined by the SurHF, is a more appropriate indicator for climate-induced warming than LSWT. Therefore, understanding the seasonal and monthly evolution of spatial SurHF variability over different parts of a large lake is important for characterizing the dynamics of many aspects of the functioning of a lake ecosystem.
The determination of SurHF is often based on one-point estimates. Relying on a one-point approach, however, can result in significant errors in the SurHF estimation of large lakes (Mahrer & Assouline, 1993) and consequently the lake heat content variability. A multiple point approach can capture some details of the SurHF spatial variability over lakes (e.g., Lee et al., 2014;Rahaghi et al., 2018;Verburg & Antenucci, 2010;Wang et al., 2014) but only provides a partial understanding of the overall energy exchange dynamics in large lakes. With increasing water body size, systematic variation of LSWT and the meteorological conditions over the lake surface can be expected and can potentially induce systematic larger spatial variability of SurHF between different parts of the lake. Even though identifying and quantifying this aspect is important for understanding the lake ecosystem dynamics, this topic has received little attention in the literature because, for most lakes, data from a sufficiently large number of points are scarce, and these points are often located close to shore.
Investigation of systematic spatial SurHF patterns over the whole surface of inland water bodies requires spatiotemporal LSWT data, over-the-lake meteorological data (both with sufficient detail), as well as multistation over-lake measurements or multiple temperature profiles for the determination of the overall energy balance variability. Spatially resolved meteorological data can be provided by numerical weather models (e.g., Kourzeneva et al., 2012;Mironov et al., 2012), while LSWT data are often available from satellite remote sensing (e.g., Politi et al., 2012;Riffler et al., 2015;Sima et al., 2013). Satellite data that are collected on a regular basis are suitable for long-term thermal analyses, although pixel sizes may be large (at the km scale) and clouds can reduce coverage. The combination of LSWT maps and meteorological data from stations located around the shore (e.g., Alcantara et al., 2010;Phillips et al., 2016) or numerical weather models (e.g., Moukomla & Blanken, 2017;Xue et al., 2015) are used to estimate the surface thermal patterns of inland water bodies. Spatial thermal patterns in lakes are also obtained by three-dimensional (3-D) hydrodynamic modeling (e.g., Bai et al., 2013;Xue et al., 2015). Unfortunately, these models rely on bulk temperatures for back radiation and nonradiative heat flux calculations, whereas LSWT is preferred (Minnett et al., 2011;Wilson et al., 2013). Lake Geneva, the largest lake in Western Europe, is a suitable study site to investigate the existence and importance of systematic spatial SurHF variability, since all the required data are available. Previous studies used one-point measurements (Vercauteren et al., 2009) or a two-point analysis (Rahaghi et al., 2018) to quantify the temporal surface heat exchange and heat content variation. However, the analysis of time series of satellite data indicated the existence of systematic warm/cold LSWT regions over Lake Geneva (Oesch et al., 2005(Oesch et al., , 2008) that could induce systematic spatial variability in SurHF over this large lake. Below, we quantify the impact of the spatiotemporal SurHF variability of Lake Geneva on its overall energy balance based on the optimized two-point calibrated SurHF model of Rahaghi et al. (2018), using the same data sets. Specifically, the following questions are addressed: (i) Does systematic spatial variability of the SurHF exist over Lake Geneva? If so, what is the typical monthly range of SurHF spatial variability between different parts of Lake Geneva using optimized and calibrated bulk formulas? Is it significant? and (ii) What are the major meteorological parameters controlling the dominant heat flux terms that determine spatial thermal patterns? In order to improve the statistical significance of the results and to investigate whether their year-to-year variability is important, this analysis was carried out for a 7-year period.
Note that certain figures and tables mentioned in the text (using the prefix S) appear in the Supporting Information section.

Study Site
Located between Switzerland and France, Lake Geneva (local name: Lac Léman) is a large, deep, crescentshaped perialpine lake with a mean surface altitude of 372 m. It is approximately 70 km long, with a maximum width of 14 km, a surface area of 582 km 2 , and a volume of 89 km 3 . The lake is composed of two basins: an eastern, large basin called the Grand Lac, with a maximum depth of 309 m, and a western, small, narrow basin called the Petit Lac, with a maximum depth of approximately 70 m ( Figure 1). It is a warm monomictic lake that is vertically stratified in density throughout the year mainly due to temperature. Occasionally, complete mixing takes place during very cold winters (CIPEL, last accessed 31 October 2019).
The lake is surrounded by the Jura Mountains in the northwest, and by the Alps in the south and, to a lesser extent, the northeast ( Figure 1). This topography creates a wide "corridor" through which two strong dominant winds pass over most of the lake surface, namely, the Bise, coming from the northeast, and the Vent, from the southwest (Figure 1). The eastern part of the Grand Lac basin, also referred to as the Haut Lac area, has surrounding mountains that can reach over 1,000 m above the lake surface within 2 km from shore ( Figure 1). Lemmin and D'Adamo (1996) showed that, on average, the western part of the Grand Lac and most of the Petit Lac experience higher wind speeds than the eastern part of the Grand Lac due to topographic sheltering by the high mountains surrounding the Haut Lac area.

Meteorological Data
Offshore meteorological data over Lake Geneva are not measured. However, the Swiss Federal Office of Meteorology and Climatology (MeteoSwiss, last accessed 31 October 2019) runs a numerical weather model, called COSMO-2 (COnsortium of Small-scale MOdeling, last accessed 31 October 2019), which provides hourly output on a 0.02°(~2.2 km) grid for a domain extending mainly over Switzerland that also includes our study site (Voudouri et al., 2017). The model has 60 vertical levels reaching 22 km into the atmosphere. It applies initial and boundary fields taken from the MeteoSwiss operational analyses archive at 0.0625°resolution. COSMO-2 data used in this study include spatiotemporal maps of wind speed (10 m above the lake), air temperature (2 m above the lake), relative humidity (2 m above the lake), cloudiness, global radiation, and air pressure. Model outputs are systematically verified by MeteoSwiss against in situ measurements from synoptic and automatic stations located in Switzerland and Europe. We used the assimilated outputs (called analysis data), based on observational data in all of Switzerland, rather than forecast data (Voudouri et al., 2017). Hourly COSMO-2 maps of the above-mentioned parameters were taken as the input for the lake SurHF estimation over the study period, that is, 1 March 2008 to 31 December 2014. . The Swiss coordinate system with km length-based units (CH1903, last accessed 31 October 2019) is used. The colorbar legend on the right indicates the water depth. SHL2 (red star, 46.45°N, 6.59°E) and GE3 (black triangle, 46.3°N, 6.22°E) are two monitoring stations in the lake that were used for the SurHF model calibration and validation. The orange dot shows the location of the Buchillon station (100 m offshore). The solid red line separates the two basins of Lake Geneva called the Petit Lac and the Grand Lac. The yellow dashed line approximately delimits the eastern part of the Grand Lac that is also referred to as the Haut Lac. Arrows indicate the direction of the strongest winds, Bise and Vent, that pass over most of the lake surface.
COSMO-2 distinguishes lakes from land by using a lake model for the momentum transfer calculation (Mironov, 2008), but COSMO-2 cannot spatially resolve LSWT. Although this might affect the meteorological patterns over the lake, short-term in situ measurements over Lake Geneva using a moving platform (Rahaghi et al., 2019), as well as long-term point measurements on this lake (Vercauteren et al., 2008) and other water bodies (e.g., Assouline et al., 2008;Solcerova et al., 2018), showed that no correlation between the patterns of air temperature and LSWT could be established. These studies indicated that the meteorological patterns over small to moderately large lakes are mainly driven by large-scale atmospheric patterns rather than by lake thermal structure. To further investigate the quality of the analysis COSMO-2 outputs for the Lake Geneva area, in situ measurements taken at four over-land (located within 2.4 km from the lake shore) and two over-lake meteorological stations were compared with the corresponding COSMO-2 data for the study period (Rahaghi et al., 2018). A high correlation coefficient (>0.96) between those measurements and the COSMO-2 data was found, with the exception of wind speed, which has higher local spatiotemporal variability. For wind speed, the cross-correlation between different stations was comparable for COSMO-2 data and in situ measurements, thus confirming the capability of the COSMO-2 model to also capture the large-scale wind patterns over Lake Geneva. In previous studies, the feasibility of the COSMO-2 data to realistically force hydrodynamic models of Lake Geneva was demonstrated by validation of model results against in situ measurements (Bonvin et al., 2013;Cimatoribus et al., 2018;Razmi et al., 2013Razmi et al., , 2014. In addition to the over-lake meteorological patterns, LSWT spatial variability can also have an effect on the SurHF estimates of large lakes (e.g., Mahrt & Hristov, 2017;Mahrt & Khelif, 2010;Rahaghi et al., 2019). Therefore, LSWT maps (section 2.3) are essential for a better quantification of lake SurHF. Riffler et al. (2015) determined LSWT for 25 lakes in and near the Alps from a long-term archive of Advanced Very High Resolution Radiometer (AVHRR) satellite imagery (~1-km pixel size). The satellite-based temperatures agreed well with the near-surface in situ measurements, corrected for the skin-bulk temperature difference (Minnett et al., 2011), with a bias and root mean square error within the range of −0.5 to 0.6°C and 1.0 to 1.6°C, respectively. This range of values corresponds favorably to another long-term LSWT calibration for Lake Geneva (Oesch et al., 2005). In the present study, we use the same data set as Riffler et al. (2015), that is, 4,384 images from 1 March 2008 to 31 December 2014, to retrieve LSWT. At times, these images can be incomplete due to cloud cover. Since pixel-wise spatially resolved hourly time series of LSWT were required for our calculations, a sequence of spatial and temporal interpolations was carried out on the satellite images. First, images with more than 70% lake coverage were selected, resulting in a total of 856 diurnal images and 308 nocturnal images. Missing pixels in these images were interpolated spatially using Barnes interpolation (Koch et al., 1983;Liston & Elder, 2006). In the next step, the hourly time series were produced using piecewise cubic Hermite polynomials. For a large surface area, the effect of diel temperature variation on the pixel-wise interpolated data is expected to be small (Woolway et al., 2016). A cross-validation analysis was performed on the interpolation algorithms by randomly and repeatedly removing some measured points ( Figure S1). A correlation coefficient of~0.99 and a root mean square error of 0.85°C were found when comparing actual measurements and the virtually interpolated points. Since LSWT and COSMO-2 maps have different spatial resolutions, that is,~1 and 2.2 km, respectively, cubic polynomial interpolation was applied to map the meteorological data into the LSWT satellite grid cells.

Spatiotemporal Variation of Meteorological Data and LSWT
Spatiotemporal maps of wind speed (U 10 ), air temperature (T a ), relative humidity (ϕ rel ), cloudiness (C), global radiation (Q sc ), air pressure, and LSWT (T w ) were used in this study. Previous studies demonstrated that SurHF is more sensitive to wind speed and LSWT than to air temperature variability and relative humidity (Austin & Allen, 2011;Rahaghi et al., 2019).
Monthly mean spatial anomalies of wind speed for March, June, September, and December, representative of the four seasons, spring, summer, fall, and winter, respectively, averaged over the study period (1 March 2008 to 31 December 2014) are plotted in Figure 2a (lake maps for all months are presented in Figure S2). In this study, seasons are defined as follows: spring (March to May), summer (June to August), fall (September to November), and winter (December to February). A distinct pattern is evident for all seasons, with the highest winds over the western part of the Grand Lac and often slightly lower winds over the Petit Lac.
Lowest wind speeds are systematically found over the Haut Lac area, the eastern part of the Grand Lac, because of the surrounding high mountains. On average, wind speeds over the western half of Lake Geneva are twice as high as those over its eastern half and are higher from September to March compared to April to August. The lake-wide mean wind speed also shows a seasonal pattern with high speeds measured during winter and low speeds during the summer.
The monthly mean LSWT distribution over the Lake Geneva surface for selected months, averaged over the period 1 March 2008 to 31 December 2014 (Figure 2b), shows that the temperatures in the Petit Lac are most of the time below the lake-wide mean LSWT whereas the temperature values are above the mean LSWT in the eastern part of the Grand Lac (see maps for all months in Figure S3). A typical pattern for the summer months (June to August) is seen in the June image of Figure 2b with a strong North-South temperature gradient across the Grand Lac, which is most pronounced in the western and central parts. During fall and winter, a band of below-average LSWT is observed along the southern shore of the Grand Lac due to shielding of solar radiation by mountains. Temperatures in the Haut Lac area of the lake remain above the lake-wide

Water Resources Research
average year-round whereas in the western part of the Grand Lac, they may be close to or slightly above the average. Lake-wide mean LSWTs follow a seasonal cycle that is characteristic of the midlatitude climate zone. The monthly mean spatial patterns of air-water temperature differences, which will be discussed below, are also plotted (Figures 2c and S4).
The temporal evolution of meteorological data and LSWT (labeled T w in the figures), smoothed with a 30day moving average, shows the mean annual cycle typical for LSWT, air temperature, and global radiation at this latitude, with maxima in the summer and minima in the winter ( Figure S5). The wind speed reaches its highest mean values during the fall to spring period. To further investigate the temporal variability of the spatial distribution of these parameters, an interquartile range (IQR) analysis based on the 25-75 percentile was carried out to determine the range of the actual values with respect to the spatial median value. To better visualize the results, the median time trend given in Figure S5 was subtracted, and to facilitate the interpretation of the IQR results, we included the 1-99 percentile, called P1-99 hereinafter, which is indicative for the total range.
The spatial distribution range of the above parameters over the whole lake, identified via IQR and P1-99 and smoothed with a 30-day running mean (Figure 3), reveals that spatial variability is important throughout the study period (2008 to 2014). Figure 3a gives the temporal evolution of LSWT spatial anomalies. The negatively skewed distribution, which is dominant in Figure 3a, is due to the colder mean LSWT values in the Petit Lac compared to the Grand Lac and the narrow cold band close to the southern shore of the lake (Figure 2b). The spatial variability in wind speed, U 10 ( Figure 3b), is considerable. The IQR band shows seasonal variations with the widest bands usually during the fall to spring period.
Comparing the median to the quartile values reveals that the wind speed distribution over the lake has a high proportion of low speeds, and therefore, the median is closer to the first quartile (Figures 3b and S5b). Significant peaks, mainly during winter, are seen in the total range of the P1-99 Figure 3. Time series of the spatial anomalies, IQR and P1-99, of lake surface water temperature (LSWT) and meteorological data for the study period (2008)(2009)(2010)(2011)(2012)(2013)(2014). The time series were smoothed with a 30-day running window, and the median (given in Figure S5) was removed. (a) LSWT (T w ), (b) wind speed (U 10 ), (c) air temperature (T a ), and (d) global radiation (Q sc ). These plots show the time evolution of the colored areas in Figure S5. Colors are identified in the legend in (c). The overbar~symbols indicate the lake-wide median values for the corresponding parameter.

10.1029/2019WR024954
Water Resources Research band ( Figure 3b). Air temperature has the narrowest IQR band with respect to the standard deviation, indicating that air temperatures over the lake are more evenly distributed ( Figure 3c). However, the total range (P1-99) suggests that often, in particular from fall to spring, significant temperature differences can be expected over Lake Geneva. The total range (P1-99) of spatial relative humidity variability is small ( Figure S6) with higher values mainly during winter. The IQR band of the global radiation is narrowest during winter and widest during summer ( Figure 3d). The total range is wide with more pronounced negative peaks mainly occurring during the summer season. This may be related to different types of cloud cover during the different seasons. Shielding by the high mountains along the southern shore in the central to eastern part of the Grand Lac may also contribute. The above results demonstrate the seasonal variability and confirm the significance of spatial variability of LSWT and those meteorological parameters that most strongly affect the SurHF dynamics over Lake Geneva.

Bulk Modeling of SurHFs
The net (total) heat flux at the air-water interface, SurHF, also called Q N [Wm -2 ], is given by where the right-hand side terms describe the flux due to absorbed solar shortwave radiation, Q sn , absorbed incoming longwave radiation from the sky, Q an , back longwave radiation, Q br , latent (evaporation), Q ev , and sensible (convection), Q co , heat fluxes. All the SurHF components are assumed positive downwards.
Recently, Rahaghi et al. (2018) evaluated various combinations of SurHF models for Lake Geneva based on the total heat content variation concept. They showed that, on the spatial and temporal scales of this study, advection of heat is consequential only within the nearshore regions of Lake Geneva and can be neglected in the rest of the lake. A two-point calibration and optimization (Rahaghi et al., 2018) was performed using long-term detailed water temperature profiles at two locations far from the shore (SHL2 and GE3 in Figure 1). The present analysis is based on the obtained optimal models, which are assumed to be valid over the entire lake. A summary of the calibrated formulas is provided in Text S1 and in Tables S1 and S2.
In order to evaluate the performance of the turbulent heat flux models (Table S2), the calculated values were compared with direct sensible and latent heat flux measurements obtained by Vercauteren et al. (2009) at the Buchillon station located 100 m off the northern shore of Lake Geneva (orange dot in Figure 1) using the eddy covariance technique. Since COSMO-2 outputs were not available for their measurement period (18 August to 26 October 2006) and the satellite pixel for their station close to shore was probably contaminated partly by land, we estimated the turbulent heat fluxes using their meteorological and LSWT measurements. The results ( Figure S7) indicated a high (>0.9) correlation coefficient between our estimates and the direct measurements. However, notable normalized (with mean) root mean square difference (NRMSD) values for our model results were found when compared with the direct measurements, that is, 48% for Q ev and 113% for Q co . The high NRMSD values are mainly due to considerable negative bias in the LSWT measurements (Vercauteren et al., 2009). Assuming a 10% bias in the LSWT data as a reasonable estimation (based on measured data by Vercauteren et al., 2009 using two different techniques) significantly improved our results, that is, NRMSDs of 25% for Q ev and 64% for Q co . These results are close to the reported NRMSD of 23% for Q ev obtained by applying the eddy covariance technique, which estimates Q ev by means of measured Q co and other meteorological data (Vercauteren et al., 2009).

Mean Annual Spatiotemporal Variability of SurHF and Its Terms
The lake-wide net SurHF and its terms, averaged over the study period (1 March 2008 to 31 December 2014) and smoothed with a 30-day running mean, are presented in Figure 4. The mean shortwave solar radiation, Q sn , follows the annual solar cycle (Figure 4a), with a maximum and minimum around the summer and winter solstices, respectively. Most years, the mean Q sn curve shows submonthly temporal fluctuations most often from the end of the spring to the summer, reflecting cloud cover variability. Atmospheric longwave radiation, Q an , is the largest positive contributor to SurHF with the second largest annual amplitude after 10.1029/2019WR024954 Water Resources Research solar radiation (Figure 4b). It is compensated for by the longwave back radiation from the lake, Q br ( Figure 4c). The Q an and Q br terms approximately follow a comparable periodic variation but are shifted in phase by half an annual cycle. The mean annual radiative heat fluxes were estimated to be 135.6 for Q sn , 292.7 for Q an , and −370.2 Wm -2 for Q br , respectively, for the period 1 January 2009 to 31 December 2014 (due to the missing data of January and February, 2008 results were excluded for the mean annual analyses).
The contributions of the nonradiative (turbulent) heat flux terms are smaller in magnitude compared to the radiative components, but they show larger variability in both time and space. On average, the evaporative heat loss (Q ev , Figure 4d) is lower during the spring and higher in the summer. The mean annual evaporative cooling for the period 2009 to 2014 is −45.0 Wm -2 . Convective heat loss (Q co , Figure 4e) is usually highest in winter when both wind speed and air-water temperature differences are high. In contrast to evaporation, which always cools the lake, convection may occasionally warm the lake during spring, mainly in April and May, when the lake-wide average LSWT can be lower than the air temperature. Overall, the mean annual convective cooling is −14.1 Wm -2 for the period 2009 to 2014.
The net SurHF (Figure 4f) follows an annual cycle with occasional exceptionally strong cooling (e.g., January 2012 and December 2013) and warming (e.g., May and June 2014) events (see the mean  Figure S8). Due to the phase shift in cooling and warming of the lake, the range of mean net SurHF from equation (1) Figure 4f) is higher than the range of each of the individual SurHF terms (<250 Wm -2 , Figures 4a to 4e). The spatial variability of IQR and P1-99 will be discussed below.

Spatial Variability of SurHF
The monthly spatially resolved variability of the net SurHF maps ( Figure 5) indicates that, on average, the SurHF in the western part of the Grand Lac is below the lake-wide mean value, whereas it is above in most of its eastern part. This nonuniformity in SurHF is accentuated during the cooling season, and the spatial pattern changes little in shape during that period. Thus, the eastern part of the Grand Lac cools less rapidly than its western part. A minimum SurHF spatial variation is found from March to May, with a spatial standard deviation, σ s , of~8 Wm -2 . In the Petit Lac, the SurHF is either above the lake-wide mean or close to it, and the least spatial variability is again seen during March to May.
Mean monthlyQ N and σ s values, as well as the mean SurHF spatial range (min/max difference), R s (Q N ), are listed in Table S3. A relatively high spatial standard deviation of >20 Wm -2 (equivalent to a temperature

10.1029/2019WR024954
Water Resources Research change of 0.88°Cy -1 in the whole water body of Lake Geneva) was found during February, October, and December (25%, 32%, and 17% of mean Q N , respectively). The results also indicate that during certain months, on average, some parts of the lake warm up, while other parts cool down. This behavior is pronounced during September and to a lesser extent in March. During these months, the mean lake-wide SurHF is low. Lofgren and Zhu (2000) estimated the monthly SurHF for Lake Huron using satellite LSWT and data from onshore meteorological stations. They found a maximum spatial standard deviation of~20 Wm -2 during August and December. A close-to-zero mean net SurHF was reported for March and September, in agreement with our results from Lake Geneva. Over Lake Constance, a total SurHF spatial range was estimated to be 60 Wm -2 during May 1989 (Schneider & Mauser, 1991). We obtain a value of 63 Wm -2 for that month over Lake Geneva. The effect of the spatial SurHF variation on the spatial lake heat content variation is addressed below.
The average spatial pattern of the net SurHF for 2009 to 2014 ( Figure 6) shows a range of >40 Wm -2 and has a mean spatial standard deviation of~13 Wm -2 . This variability emphasizes that a single-point analysis can lead to sizable errors in the estimation of the SurHF and hence the heat budget of a large water body (Rahaghi et al., 2018). A small negative mean spatiotemporal heat flux of −1 Wm -2 (indicated in Figure 6) is found for the entire period considered, that is, the lake slightly cooled during the period from 2009 to 2014, if we assume that SurHF dominates the lake's energy budget. This small value confirms the validity of the methods and formulas employed in this study, since a close-to-zero value is expected for the mean spatiotemporal SurHF of Lake Geneva. However, due to various sources of uncertainty, for example, errors in data retrieval and model assumptions, the obtained mean value (−1 Wm -2 ) must be interpreted with caution. The physical parameters controlling the obtained spatial SurHF anomalies and the effect of such variability on the estimation of lake heat content are discussed in the following sections.

Statistical Analysis of the Spatial Heat Flux Variability
In order to determine the contribution of the individual SurHF terms to the observed spatial SurHF anomalies ( Figures 5 and 6), we computed statistics relating the spatial variability of the net SurHF to that of the different SurHF terms. The correlation coefficient and root mean square difference (RMSD) between spatially resolved hourly maps of different SurHF terms and the total SurHF were calculated. The results were first smoothed with a 30-day running mean window and then averaged for each hour in the year over the corresponding time during the 6-year period from 1 January 2009 to 31 December 2014 (for leap years, the additional day in February was removed). Higher correlation coefficients and smaller RMSD values indicate a more significant contribution of the corresponding SurHF term to the spatial variability of the net SurHF. The radiative heat flux components, Q sn , Q an , and Q br , have relatively low correlation coefficients and high RMSDs (Figure 7) pointing to the low probability of spatial variability of these terms. Evaporation (latent heat flux), Q ev , patterns have the highest correlation coefficient (Figure 7a) and the lowest RMSD (Figure 7b) for all months, followed by convection, Q co . The narrow, color shaded areas (95% confidence intervals) of Q ev and Q co curves in Figure 7a indicate that the obtained high correlation coefficients are statistically significant. Therefore, spatial evaporation and convection variations are mainly responsible for the observed spatial anomalies of SurHF over Lake Geneva. This is in agreement with the findings in other large lakes, for example, Lake Superior (Xue et al., 2015) and Lake Huron (Lofgren & Zhu, 2000). However, there is more spatial variability in global radiation and air temperature over Lake Huron due to its larger surface area (~10 times greater than Lake Geneva). Therefore, the spatial variability of the net radiative flux on Lake Huron is more than two times higher than that of Lake Geneva. As is the case for Lake Geneva, significant spatial variability of evaporation was also reported for those lakes, as well as for other large lakes (e.g., Moukomla & Blanken, 2017;Verburg & Antenucci, 2010).
Figure 7 also indicates that the correlation and RMSD between the SurHF terms and net SurHF vary throughout the year. Evaporation has a high correlation (>0.8) with net SurHF from July to February,

Water Resources Research
which is reduced to about 0.7 in March and April. In March, and to a lesser extent in April, the RMSD curves converge due to the overall small spatial variability of the net SurHF ( Figure 5) and the various SurHF terms (not shown). The convective (sensible heat flux), Q co , component demonstrates a behavior similar to Q ev from November to February, while it deviates from Q ev from March to June. Generally, a change of the SurHF forcing regime during the March to June period is observed for the net SurHF spatial patterns compared to the rest of the year ( Figure 5) and for the analysis of SurHF components (Figure 7, the Q br correlation curve is a clear example). This will be further discussed below.
Time series of the IQR of the spatial variability of the net SurHF and its components with their median values removed ( Figure S9) show that the spatial variability of the radiative components ( Figure S9a) is small throughout the study period and hardly changes with time, thus confirming the correlation analysis results above. The spatial variability of SurHF in IQR is dominated by the spatial variability of the two turbulent flux components, in particular evaporation, Q ev ( Figure S9b). The seasonal variability of Q ev is obvious with smaller spatial variability occurring during spring. However, a year-to-year variability is also seen in this pattern reflecting changing large-scale weather conditions. During most years, the convective heat flux, Q co , has a low spatial variability except from late fall to early spring, when it may reach spatial variability values similar to those of Q ev . Individual events of strong spatial SurHF variability such as in early 2012 and late 2013 are almost exclusively caused by the variability of the turbulent heat flux components. The P1-99 distribution (not shown) is similar to the IQR distribution but spread out over a wider range. This analysis demonstrates that the spatial variability of turbulent heat fluxes significantly contributes to the spatial variability of SurHF over Lake Geneva.

Effects of Spatial Variability of Meteorological Forcing on SurHF
In order to determine the major meteorological parameters controlling the observed spatial SurHF variability patterns, we computed the correlation coefficients of hourly spatial patterns (over 7 years) of each meteorological parameter and the net SurHF at the corresponding time. The hourly (spatial) correlation coefficients thus obtained were first smoothed with a 30-day running mean, and then the correlation coefficients for a given hour in the year were averaged over the 6 or 7 years considered (for leap years, the additional day in February was removed). The results are shown in Figure 8. Spatial LSWT (T w ) and wind speed (U 10 ) are negatively correlated with the spatial SurHF pattern. Except for spring (March to May), wind speed has a correlation of >0.6 and is the dominant meteorological parameter controlling net SurHF spatial variability. Global radiation, Q sc , reaches a maximum correlation coefficient of around 0.5 during March and April, the months with the lowest spatial variability in net SurHF ( Figure 5).
The air-water temperature difference (T w − T a ) is close to zero in early spring over a large portion of the lake (March and April in Figures 2c and S4). The wind speed is also relatively low during that period, in particular during April (Figures 2a and S2). As a result, the spatial SurHF variability due to evaporation and convection is lowest in early spring, as indicated by the low spatial correlation coefficients during this time (Figure 7). During the March to May period, when wind speeds are low (Figures 2a and S2), LSWT and the air-water temperature difference are the most important parameters, with correlation coefficients of approximately 0.6 ( Figure 8). During this period, the convective net cooling effect due to the air-water temperature difference (T w − T a ) is larger in the eastern part of Lake Geneva than in its western part (not shown). Thus, the wind sheltering effect in the eastern part of the lake is partially compensated by the air-water temperature

10.1029/2019WR024954
Water Resources Research difference during spring, and consequently, there is a less noticeable difference in SurHF between the eastern and western parts of the Grand Lac ( Figure 5). Schneider and Mauser (1991) suggested the same reason for the estimated spatial variation of SurHF over Lake Constance during May 1989. Using point data over three small lakes, Granger and Hedstrom (2011) found that wind speed was the most significant factor causing turbulent SurHF variations, followed by land-water temperature difference and then land-water vapor pressure. They also reported a minimum in turbulent cooling during spring, as was found in the present study. On the other hand, an investigation in Lake Superior, a much larger lake than Lake Geneva, showed that spatial patterns of evaporation follow the synoptic-scale air masses traveling over the lake (Spence et al., 2011).
The patterns in Figures 2c and S4 demonstrate that the convective net heating effect due to the air-water temperature difference is smaller in the eastern part of Lake Geneva during March to June compared to the western part of the lake. On average, LSWT (Figures 2b and  S3), and hence longwave cooling, is lower in the Petit Lac compared to the Grand Lac. Furthermore, the effect of strong winds on convective cooling is partially compensated for by a higher (T a − T w ) in the Petit Lac. This explains the higher SurHF in the Petit Lac compared to the western Grand Lac (Figures 5 and 6).

Spatial Variability of Atmospheric Boundary Layer Stability
It was shown above that evaporation and convection, given by equation (S4) in Table S2, are the most important heat flux components affecting the spatial variability of net SurHF. The main variable in these terms is the Monin-Obukhov atmospheric boundary layer (ABL) stability parameter, ζ, which characterizes the relative contribution of buoyancy and wind shear to the turbulence generation over the lake surface (Yusup & Liu, 2016). When ζ < 0, the ABL is unstable, and turbulent cooling of the surface is higher. In contrast, under stable ABL conditions (ζ > 0), the turbulent heat exchange is less important. ABL stability is determined by several meteorological parameters, such as wind speed (U 10 ), LSWT (T w ), air-water temperature difference (T a − T w ), and relative humidity (ϕ rel ). To examine the combined effect of the dominant controlling factors, that is, U 10 , T w , and (T w − T a ), the mean temporal variation and the mean spatial pattern of the ABL stability parameter ζ for different months were calculated.
The annual time series of the median and the spatial anomalies of the ABL stability parameter, ζ, are shown in Figure 9. The results (smoothed with a 30-day running mean) demonstrate that the ABL is unstable over Lake Geneva, except for a stable period from early March to early June. Woolway et al. (2017), using data from different lakes, observed that unstable ABL conditions are least common during spring. The transition to a stable period in spring coincides with the strong change in the values of the various correlation coefficients observed before in Figures 7 and 8. Evaporation, as discussed above, was at its minimum during the stable period (ζ > 0). Our results indicate that on an annual basis, the ABL is unstable~74% of the time over Lake Geneva and that the unstable period accounts for~90% of total evaporation. Woolway et al. (2017) reported unstable ABL conditions~72% of the time (annual basis). Persistent unstable conditions were observed over tropical lakes (Verburg & Antenucci, 2010;Woolway et al., 2017).
The range of spatial anomalies of ζ as indicated by IQR and P1-99 is approximately constant in time except in late February and early March, when the stability changes rapidly from unstable to stable and the

10.1029/2019WR024954
Water Resources Research spatial variability is reduced (shaded areas in Figure 9). Since the ABL is more stable during spring (March to May), the air temperature is slightly higher than the LSWT over most of Lake Geneva (Figures 2c and S4), resulting in reduced turbulent heat fluxes, Q ev and Q co . In early spring, this stability is due to intense solar radiation, colder surface layer water temperatures, and a weaker stratification than in summertime. Comparing T w and T a spatial patterns (not shown) suggests that the spatial distribution of air-water temperature differences over Lake Geneva, and hence the spatial variability of the atmospheric thermal boundary layer, is predominantly controlled by LSWT, not air temperatures.
Maps of the mean ABL stability parameter, ζ, based on the 7 years of analysis were created for each month (Figure 10), as described above for SurHF ( Figure 5). The spatial variability of ζ is large compared to its mean throughout the year (Figure 10) for all months. The narrowest range of spatial variability is found during the period of stable stratification (March to May). This corresponds to the narrowest range of spatial variability of SurHF observed in Figure 5. During the unstable period (ζ < 0) from July to February, the ABL is less unstable in the western Grand Lac where the SurHF monthly mean value is at its lowest ( Figure 5). Thus, in this area, the less unstable the ABL is, the greater the heat loss is from the lake. This agrees with Yusup and Liu (2016), who found that the maxima of turbulent heat fluxes did not occur under the most

10.1029/2019WR024954
Water Resources Research unstable conditions but rather at low values of ζ. On average, the ABL is more unstable in the eastern Grand Lac where wind speeds are generally lower than in the western Grand Lac. In the eastern part of Lake Geneva, gradients from the less unstable southern part to the more unstable northern part are evident for some months. Yusup and Liu (2016) reported that low wind speeds usually correspond to very unstable ABL conditions as seen in the eastern Grand Lac. They also mentioned that due to mixing in the near surface water layers caused by wind stirring, the feedback mechanism between ABL stability and turbulent heat fluxes may at times make the interpretation of the relationship between these parameters difficult. For the Petit Lac, it is more difficult to establish a clear link between ζ and SurHF.
The ABL stability concept is part of the Monin-Obukhov similarity theory (equation (S4) in Table S2). Yusup and Liu (2016) suggested that the Monin-Obukhov similarity theory may have limited validity beyond the stability range of |ζ| < 1. In a recent study, Gao et al. (2018) reported that the temperaturehumidity similarity assumption is more uncertain under low wind shear conditions. Studying the effect of thermal inertia and water advection on the Monin-Obukhov scaling, Assouline et al. (2008) also concluded that for the case of small advection and large thermal inertia, as seen in Lake Geneva in the spring, the transport efficiency of humidity is higher than that of temperature, whereas they are usually assumed to be equal in the similarity theory (equations (S4g) and (S5b) in Table S2). This can be a source of error in the presented results.

Spatial Clustering Using k-Means
To further explore similarities in the SurHF spatial patterns and those of meteorological parameters, the kmeans clustering method (Jain, 2010) was applied to the seasonal data. In the present analysis, this method is used to spatially partition SurHF and the meteorological parameters that affect SurHF spatial variability into clusters. The number of clusters used (k ≥ 1) is arbitrary for a given data set. However, the overall variance of the data set decreases with increasing k, and the optimal value of k is usually taken at the point where the variance reduction rate becomes relatively small, because the variance first reduces rapidly, then more slowly (e.g., Tibshirani et al., 2001). The method was applied to the hourly data sets of the net SurHF (Q N ) and the three meteorological parameters that affect its spatial variability the most, that is, LSWT (T w ), wind speed (U 10 ), and global radiation (Q sc ). Results for spring and fall ( Figure 11; results for the other two seasons are presented in Figure S9) show that the net SurHF splits into four clusters (Figures 11a and 11b), thus indicating four lake areas within which SurHF values vary little from the mean of the cluster. These compare well with the pattern seen in Figure 5. There are only minor differences between the patterns in the different seasons (Figures 11a,11b,S9a,and S9b). LSWT and wind speed also split into four clusters in spring (Figures 11c and 11e, respectively). The actual pattern of LSWT more closely resembles that of the clustering of the net SurHF.
Due to convective surface cooling, and therefore strong mixing at the surface, only two LSWT clusters (one in the Grand Lac and one mainly in the Petit Lac) were obtained for the fall season (Figure 11d), while there is more consistency between the SurHF clustering ( Figure 11b) and that of the wind speed (Figure 11f) during that season. Therefore, wind speed, which is a key parameter in the evaporative and convective SurHF terms (equation (S4)), is significant in controlling SurHF spatial variability during fall. Air-water temperature differences (T w − T a ) are also higher in fall compared to spring (Figures 2c and S4) when lower values of (T w − T a ) weaken the effect of the convective and evaporative heat fluxes. As a result, the contribution of the spatial variability of the longwave emission to the observed variability in the total SHF becomes more significant during spring. The similarity between the net SurHF clusters (Figure 11a) and that of LSWT ( Figure 11c) confirms this.
In all seasons, the global radiation, Q sc , splits into three clusters (Figures 11g, 11h, S9g, and S9h) that hardly change in shape or extent. The mean centroid values of the three clusters increase from west to east by~5, 5.5, 1.5, and 1 Wm -2 during spring, summer, fall, and winter, respectively. However, only the spatial variability of springtime is comparable with the monthly standard deviation of net SurHF patterns ( Figure 5 and Table S3). The pattern of three clusters probably reflects the effect of the mountainous topography surrounding the lake and cloud cover. In general, the results of the seasonal k-means analysis agree with the temporal correlation coefficient trends presented in Figure 8 and coincide with spatial variability in the SurHF patterns ( Figure 5). For most of the year, spatial variability is organized in clearly distinguishable zones 10.1029/2019WR024954

Water Resources Research
aligned along the lake axis from the west to the east, thus emphasizing that spatial variability in Lake Geneva is statistically significant.
A source of error in the above results (Figures 4-11) may be related to the horizontal heterogeneity due to the relatively large pixel size of LSWT (~1 km) and meteorological data (~2.2 km). The area-averaged SurHF may be overestimated or underestimated, particularly when air-water temperature differences are small (Mahrt & Khelif, 2010). The change of the ABL regime from unstable to stable and vice versa results in a drastic change of bulk transfer coefficients (Deardorff, 1968). Therefore, we expect higher uncertainties during near-neutral conditions at the end of February and the beginning of March ( Figure 9). The errors mentioned above become more significant when the absolute value of the SurHF is at its minimum in the springtime. However, regardless of these uncertainties, the results of this study, including the spatiotemporal variation of SurHF (Figures 4-6) and the change of the dominant controlling regimes (Figures 7-9), are not affected. Better quantification of the associated uncertainties requires higher-resolution maps of LSWT and meteorological data, which are not available for Lake Geneva.

Spatial Variability of Lake Surface Thermal Energy and Lake Heat Content
The heat content variation in the water column over the full lake depth is due to the net energy flux into it over time. In order to determine whether and how the effect of the spatial variability of the SurHF affects the heat content of the lake, the temporal variation of the water column thermal energy is taken as (Fink et al., 2014;Nussboim et al., 2017;Van Emmerik et al., 2013) Here, the net input energy, ΔG m [Jm -2 ] into the lake, is calculated integrating the net SurHF, Q N [Wm -2 ], for a given period (equation (1)). 3-D numerical simulations of Lake Geneva show that equation (2) is valid in regions far from the shore and the main river inflow and for the spatial and temporal scales of the present study Rahaghi et al., 2018). In Rahaghi et al. (2018), this equation, representing the energy balance in the water column, was applied to the data collected at two stations (SHL2 and GE3 in Figure 1) in order to optimize and calibrate the net SurHF model for Lake Geneva.
To further investigate the validity of the obtained SurHF spatial variability ( Figures 5 and 6), the heat content variation was calculated for two regions (shaded areas in Figure 12a), where long-term temperature profiles were taken in previous studies (as discussed below). One region (denoted EG) is located in the western Grand Lac where SurHF is below the lake-wide mean and the other one (PG) in the eastern Grand Lac where SurHF is above the lake-wide mean ( Figure 6). Time series of Q N were averaged over a small area in these regions, instead of using single points, in order to obtain more representative estimates for these regions. For the 2009 to 2014 period, the spatiotemporal average net SurHF, Q N , is 3.2 Wm -2 (above the mean value of −1 Wm -2 obtained for the whole lake and for the same period) in the EG region and −6.7 Wm -2 (below the mean value) in the PG region. Time series of the heat content calculated with equation (2) for the regions EG and PG (Figure 12b) indicate that the heat content changes in the two areas follow the same seasonal pattern. It should be noted that while there is an increase in heat content in 2013 and 2014 in the EG region, heat is being lost during the entire study period in the PG region. However, in the EG region, winter losses are always larger than summer gains. Thus, spatial variability of SurHF results in the spatial variability of the heat content in the lake. Unfortunately, no measured water temperature profiles in these regions are available for the period 2008 to 2014 to compare these results.
Monthly temperature profiles in these areas were taken for the period July 1986 to June 1995 at several stations in the lake. From these measured temperature profiles, the total heat content of a water column in the lake, G o [Jm -2 ], is obtained by Figure 12. (a) Location of long-term monitoring stations P1 to P4 in Lake Geneva; see legend for details, (b) heat content evolution in the lake based on SurHF estimation for the period 4 July 2008 to 31 December 2014. The shaded regions in (a), denoted by EG and PG, represent the areas used for the net SurHF, Q N , estimation (see the text for more details), and (c) temporal evolution of heat content estimated by approximately monthly temperature profiles at stations P1-P4 for the period 4 July 1986 to 14 June 1995. In order to allow a comparison of the curves in panels (b) and (c), the horizontal axes of the two plots are made for the same number of years, with both starting on 4 July. For a better comparison, all curves in (b) and (c) are plotted relative to their heat content value on the initial date.
where ρ w and C p,w are the density [kgm -3 ] and the specific heat capacity of water at constant pressure [Jkg -1 K -1 ], respectively, and T(z,t) represents vertical temperature profile [°K] at time t down to water column depth z = H [m]. The change in heat content from time t 1 to time t 2 can then be quantified by between two consecutive profiles. These calculations were carried out for stations P1 and P2 in the EG area and P3 and P4 in the PG area. It can be seen (Figure 12c) that the seasonal variation and the amplitude of the variation are similar at all four stations and are comparable to those in Figure 12b for the net input energy, ΔG m [Jm -2 ] into the lake. The temporal change in heat content at P1 and P2 in EG is systematically above that of stations P3 and P4 in the PG area. This result is consistent with the difference in the time variation of the heat content in the EG and PG areas computed from the SurHF (Figure 12b) and confirms that the heat content in the eastern Grand Lac is higher than that in the western Grand Lac. Furthermore, it shows that the SurHF (Q N ) for Lake Geneva is the dominant source of heat input into the lake and that other processes such as advective heat flux play only a minor role, as was demonstrated previously (Rahaghi et al., 2018). Water mass movement in the lake is not sufficient to eliminate these horizontal gradients. This also indicates that in a large lake the effect of spatial variability of heat flux on the heat content may be important and that estimating heat content dynamics of a lake from a single station may not be representative for the whole lake. It is suggested that the trend toward higher heat content values and the variability between the eastern and western Grand Lac regions over the study periods (Figures 12b and 12c) is related to the effects of climate change (Lemmin & Amouroux, 2013). The difference between the rates of heat content variation obtained from model outputs during 2008-2014 ( Figure 12b) and from in situ temperature profiles during 1986-1995 ( Figure 12c) can mainly be due to shifting climate patterns as seen in the long-term lake temperature profile data at SHL2 and GE3 (CIPEL 2018, last accessed 31 October 2019): 1986 to 1995 was a period of overall continuous lake warming, whereas several strong cooling events occurred during the period 2008 to 2014. Uncertainties in the calibrated model, or errors in the input data, that is, COSMO-2 and satellite, could also play a role.

Summary and Conclusions
The primary objective of this study was to determine whether systematic spatial variability of the SurHF exists between different parts of Lake Geneva, the largest lake in Western Europe. If so, can it be quantified? Does it follow a seasonal pattern? And, finally, is it significant? In this study, we addressed these questions by calculating the spatiotemporal distribution of net SurHF and its five main components over a 7-year period for Lake Geneva, using the calibrated SurHF model (Rahaghi et al., 2018) based on observations at two locations. This permitted better identification and quantification of the systematic spatial variability of SurHF, the causes of this variability, and the seasonal pattern of both.
A key finding of our study is a consistent monthly spatial variability of up to approximately ±40 Wm -2 ( Figure 5) between different parts of the lake, which is due to significant spatial patterns of meteorological forcing and LSWT. From the analysis, it became obvious that the observed spatial variability of SurHF is mainly determined by the spatial variability of the latent and sensible heat fluxes (Figure 7). These fluxes are controlled by the spatial variability of wind patterns, LSWT, and air-water temperature difference, all of which affect the ABL stability and hence the air-water heat exchange dynamics.
The spatiotemporal variation of these meteorological forces follows a seasonal pattern. There is a noticeable change of the main controlling forcing in spring compared to the rest of the year. During springtime, the LSWT variability is more important (Figure 8), whereas from the summer to the winter period, the spatial variability of SurHF values predominantly reflects the seasonal variability of wind patterns, which are affected by the topographic sheltering of parts of the lake.
A seasonal change of the controlling regime was evident in the ABL stability conditions curve (Figure 9). Such a regime change is explained by the higher potential of heat absorption during the 10.1029/2019WR024954

Water Resources Research
onset of stratification (March to April) in spring and seasonal variation of the air-water temperature difference. The results also indicate that over Lake Geneva, the ABL is on average only stable during the springtime; during the rest of the year (74% of the time), it is statically unstable. The rate of evaporation is lower under stable ABL conditions in spring, because only about 10% of evaporation occurs during that period.
By carrying out the analysis over the relatively long period of 7 years, it was possible to confirm the long-term consistency and the statistical significance of the obtained monthly and seasonal spatial SurHF patterns. It was also found that there was little year-to-year change in the systematic spatial variability between different parts of the lake. The observed systematic spatial variability of SurHF and heat content in the lake can have consequences for the dynamics of the lake ecosystem.
The spatial resolution of the area-averaged (pixel-wise) analyses was restricted by the LSWT satellite pixel resolution (~1 km). However, our results are dependable in terms of showing for Lake Geneva the range and distribution of the systematic spatiotemporal variability of SurHF and surface thermal energy (heat content) between different parts of this large lake. A similar systematic spatial SurHF variability can be expected over other large lakes if a systematic spatial variability of atmospheric forcing exists.
In order to further investigate the spatial variability of energy exchange dynamics at a higher spatial resolution, more detailed spatiotemporal in situ measurements of the controlling parameters in both air and water would be helpful. However, carrying out such a campaign over large lakes or other large water bodies would be difficult. Therefore, a detailed analysis of the importance of spatial variability, such as the one presented here, may guide the planning/optimization of a field sampling strategy.