Characteristic time scales of evaporation from a subarctic reservoir

Water bodies such as lakes and reservoirs affect the regional climate by acting as heat sinks and sources through the evaporation of substantial quantities of water over several months of the year. Unfortunately, energy exchange observations between deep reservoirs and the atmosphere remain rare in northeastern North America, which has one of the highest densities of water bodies in the world. This study characterizes the dynamics of turbulent heat fluxes by analysing in‐situ observations of a compact and dimictic reservoir (50.69° N, 63.24° W) located in a subarctic environment. The reservoir is characterized by a mean depth of 44 m and a surface area of 85 km2. Two eddy covariance (EC) systems, one on a raft and one onshore, were deployed from 27 June 2018 to 12 June 2022. The thermal regime of the reservoir was monitored using two vertical chains of thermistors. Results indicate a mean annual evaporation rate of 590 ± 66 mm, which is equivalent to ≈51% of the annual precipitation, with 84% of the evaporation occurring at a high rate from August to freeze‐up in late December through episodic pulses. It was difficult to close the energy balance because of the magnitude and the large time lag of the heat storage term. To circumvent this problem, we opted to perform calculations for a year that started from the first of March, as the heat storage in the water column was at its lowest at that point and could thus be ignored. From June to December, monthly Bowen ratios increased from near‐zero negative values to about 1.5. After September, due to smaller vapour pressure deficits, latent heat fluxes steadily decreased until the reservoir had a complete ice cover. Opposite diurnal cycles of sensible and latent heat fluxes were revealed during the open water period, with sensible heat fluxes peaking at night and latent heat fluxes peaking in the afternoon.


| INTRODUCTION
Reservoirs are subject to significant water level fluctuations due to evaporation (latent heat flux, LE), which is not only a key component of reservoir mass and energy balances, but can affect some vital functions such as freshwater supply, irrigation, hydropower and navigation (Friedrich et al., 2018). In some arid regions of the world, structural measures are put in place to limit evaporation, such as floating balls or lattices (Assouline et al., 2011;Assouline & Narkis, 2021). Evaporation is intangible and therefore a challenging hydrological flux to measure, making it difficult to fully understand its magnitude and controlling mechanisms.
Several studies have focused on reservoir evaporation, but few in cold regions. Therefore, there is a need for additional in situ observations to further improve our understanding of the processes involved and document their characteristic time scales. Tanny et al. (2008) reported an average evaporation rate of 5.5 mm day À1 from July to September in the small Eshkol reservoir in Israel (33 N), which has a hot and arid climate. Further north, in the shallow Eastmain-1 reservoir in Canada (52 N), Strachan et al. (2016) found the evaporation rate to be 3.1 mm day À1 between August and October. These studies indicate that even in cold regions, reservoir evaporation can be substantial. Additional studies over reservoirs are nevertheless necessary to improve the parameterization of water-atmosphere exchanges (Mackay et al., 2017).
In cold climates, dimictic reservoirs undergo two turnover periods per year. Their thermal regime typically evolves into three successive phases (Cole & Weihe, 2016). During the ice cover phase, ice acts as a lid over the water body, preventing direct interactions between the atmosphere and the water column. Latent heat fluxes tend to remain low during this period . After the ice break-up period, the spring turnover takes place, resulting in an overall vertical mixing of the entire column. The heat storage period then starts and lasts until the end of summer. Latent heat fluxes remain low during this phase, with frequent stable atmospheric stratification. The third and final phase corresponds to the heat release period. This is characterized by a decline in water temperature due to a substantial release of energy into the atmosphere through turbulent heat fluxes that are high and sustained day and night (Blanken et al., 2011).
While evaporation varies seasonally in response to the three thermal phases, it also fluctuates on smaller time scales in response to meteorological forcing. For instance, incoming shortwave radiation causes latent heat fluxes to peak during the day (Lensky et al., 2018).
The atmospheric demand for water vapour, driven by wind speed and vapour pressure deficit, is also known to modulate evaporation in water bodies (Pérez et al., 2020). Evaporative demand can vary within a single day. For instance, changing wind direction can lead to a reduced or enhanced sheltering effect, increasing or decreasing evaporation rates (Markfort et al., 2010;Venäläinen et al., 1998). Evaporation can also vary over the course of a few days, due to passing synoptic systems that can generate sustained evaporation (Laird & Kristovich, 2002;Spence et al., 2013). Blanken et al. (2000) found that 50% of annual evaporation over the Great Slave Lake occurred over only 25% of the year through episodic evaporation water losses.
Moreover, thermocline depth and intensity, which depend in part on the reservoir morphometry (Gorham, 1964), influence turbulent heat fluxes by limiting or enhancing the energy available in the upper water layers. Indeed, Piccolroaz et al. (2015) identified a positive feedback between the lake surface temperature and the stratification dynamics of Lake Superior, Canada. The timing of evaporation occurs at different scales and remains poorly documented or correlated with physical drivers (Beck et al., 2018) for deep reservoirs in the boreal zone of northeastern Canada.
Northeastern America is one of the densest regions of lakes and reservoirs around the world (Downing et al., 2006). These lakes and reservoirs are considered to be climate sentinels (Adrian et al., 2009;Williamson, Saros, & Schindler, 2009) as well as integrators and regulators of climate change (Williamson, Saros, Vincent, & Smol, 2009). Wang et al. (2018) showed that modifications in surface energy allocation under warmer climate conditions will accelerate global lake evaporation. Moreover, in-situ evaporation observations are needed to develop and improve lake models (McJannet et al., 2017) for future climate estimates, particularly in remote areas.
There is a lack of direct in-situ measurements of turbulent heat fluxes over reservoirs in remote northern regions. The central goal of this study is to characterize the temporal dynamics of turbulent heat fluxes over a compact (i.e. a low surface to volume ratio) and dimictic hydropower reservoir located in a subarctic environment. To achieve this goal, two flux towers, one mounted on a floating raft deployed during the ice-free period and one on a nearby shoreline, were deployed to measure surface fluxes using the eddy covariance method from June 2018 to June 2022. The specific objectives were to quantify turbulent heat fluxes at daily, monthly, and annual time scales, and to identify the key processes and surface energy budget terms that govern LE at each time scale. The paper is organized as follows. We first introduce the study site and measurement methods.
Then, we describe the meteorological conditions over the whole study period and the driving factors of turbulent heat fluxes for each time scale. Finally, uncertainties underlying the flux calculations are discussed, including a description of the closure of the energy balance.

| Study site
The study site is located at the southern tip of the Romaine-2 hydropower reservoir (50.68 N, 63.25 W; Figure 1a), 80 km north of the city of Havre-Saint-Pierre, which is on the North Shore of the Gulf of St. Lawrence, Quebec, Canada. The Romaine-2 reservoir is a dimictic water body that is ice-free from May to December. The regional climate has a mean annual air temperature and precipitation of 1 C and 1167 mm, respectively. Precipitation was estimated using our measurements and Environment Canada data recorded near our study site. This is typical of the subarctic (Dfc) Köppen-Geiger climate classification type (Beck et al., 2018). The reservoir, which flooded in 2014, drains a 14 351-km 2 area that is mostly covered with a spruce-moss forest. At its southern end, the reservoir is about 1 km-wide, and sits at an elevation of 244 m above mean sea level. When full, it has a surface area of 85.6 km 2 and a maximum depth of 101 m (mean depth of 44 m). It has an elongated north-south shape and steep banks channelling the near-surface winds. The study period extended from 27 June 2018 to 12 June 2022.  (Fournier et al., 2021;Pierre et al., 2022). At both sites, the eddy covariance technique was used to calculate 30-min averaged turbulent heat fluxes from raw 10-Hz turbulence data.

| Raft and shore flux towers
F I G U R E 1 (a) Location of the Romaine-2 reservoir in North America (red square). (b) Satellite image (Sentinel-2A, 26 January 2022) of the Romaine-2 reservoir in winter. The black rectangle indicates the area of the experimental setup (c). (c) Elevation map indicating precisely where the spillway, intake, flux towers and thermistor chains (indicated by TC) are located. The red contour lines represent the flux footprints areas (90% and 80% for the raft and for the shore station respectively, refer to Section 2.2).
F I G U R E 2 The two eddy covariance measurement setups deployed on and around the Romaine-2 reservoir. Raft station (left, photo taken looking west) and shore station (right, photo taken looking north) with the following instruments: (a) net radiometer, (b) air temperature probe, (c) accelerometer, (d) sonic anemometer and infrared gas analyser, and (e) propeller anemometer. Photos were provided by Hydro-Québec and taken on 13 June 2019, at approximately 2 m below the maximum reservoir level.
The raft was anchored between two islands to offer protection ometer measured all terms of the radiation budget, namely incoming and outgoing shortwave and longwave radiation, but since the instrument was installed on the shore, it did not report the radiation emitted/reflected by the reservoir.
Flux footprints were calculated for both eddy-covariance set-up following Kljun et al. (2015). Both flux towers displayed a local flux footprint, i.e. an area of a few hectares ( Figure 1c). Moreover, the footprint of the raft always remained above the water surface, no matter what the wind direction and water level of the reservoir were.
On the other hand, the footprint of the flux tower installed on the shore encompassed the water and the land surfaces depending on the wind direction and the reservoir water level. However, only data with wind direction from the reservoir were considered (Section 2.3).
Therefore, the results reported in this study represent turbulent heat fluxes from the reservoir only.
The atmospheric variables (wind speed, air temperature, vapour pressure deficit, etc.) measured at the raft and at the shore station were merged into a single dataset as follows: the raft data were kept when available, otherwise the shore data were retained.
A complete year-long net radiation time series was created using the following strategy. We used the net radiation that was measured on the raft when it was deployed (June-October). When ice was present on the reservoir, it was typically covered in snow. We therefore assumed that the emitted/reflected radiation measured at the shore station during that time was equivalent to what would have been measured on the reservoir. During the transition periods, we estimated the reflected radiation on the reservoir with a simple albedo formula based on that by Patel and Rix (2019). The longwave radiation emitted during the transition periods was estimated from the Stefan-Boltzmann law. This was done by taking an emissivity of 0.995 and extrapolating the water temperature profile in the top meter of water to estimate the water skin temperature for periods when no direct measurements were available.

| Turbulent heat fluxes
Sensible (H) and latent heat fluxes (LE) (both in W m À2 ) were calculated using the covariance between the vertical wind speed w (m s À1 ) and air temperature T (K), as well as between w and the specific humidity q (kg kg À1 ), such that: where ρ a is the dry air density (kg m À3 ), c pa is the specific heat of humid air (J kg À1 K À1 ) and L v is the latent heat of vaporization (J kg À1 ). Here, primes denote fluctuations from the 30-min average, indicated by an overbar.
Raw wind velocity data from the raft were first corrected for wave-induced motion, using the accelerometer data, following the method proposed by Miller et al. (2008). In short, the apparent wind measured by the sonic anemometer was corrected to account for Euler angles, angular velocities and linear accelerations monitored by the accelerometer. Then, for the shore and raft flux towers, we processed the turbulence data using the EddyPro ® software, version 7.0 (LI-COR Biosciences, USA). In doing so, we applied time-lag compensation, linear detrending, double rotation approach (Baldocchi et al., 1988;Wilczak et al., 2001), density fluctuation compensation (Webb et al., 1980), spike removal (Papale et al., 2006), and other statistical tests (Vickers & Mahrt, 1997). Poor-quality data were flagged (Mauder & Foken, 2011) and removed. Data from the raft and shore stations were aggregated into one dataset by favouring data with the best quality criteria (Mauder et al., 2013). Note that shore data were retained only when winds originated from the reservoir. To complete the final dataset, gap-filling was implemented based on the method developed by Reichstein et al. (2005).
Over the whole study period (1447 days), 57% of the turbulent flux data had to be gap-filled since the raft was only deployed from June to October and the shore flux tower was frequently exposed to winds from the surrounding land. We assessed flux errors by applying the Finkelstein and Sims (2001) random uncertainty method.

| Water temperature and transparency
The Water transparency was also periodically measured using a Secchi disc. The mean Secchi depth (SD) was 4 ± 0.04 m. Following Koenings and Edmundson (1991), we then assumed that the mean product of SD Â K d was 2.28, thus we estimated the vertical attenuation coefficient K d to be 0.57 m À1 for water of moderate transparency. Consequently, 50% of the incident energy flux density is absorbed in the first 1.2 m of water and 99% over the first 8.1 m of water.

| Energy balance ratio
The energy balance ratio (Feng et al., 2016) allows us to quantify the energy balance closure and therefore to correct underestimated turbulent heat fluxes (see Section 3.3). It is defined as follows: where R n is the net radiation and ΔH S is the heat storage variation along the water column, defined as where ρ w is the water density (kg m À3 ), c pw is the specific heat of water (J kg À1 K À1 ), ΔT w is the water temperature difference between two timesteps (K), Δt is the period (30 min), and h is the water layer thickness (m). To determine the EBR, ΔH S was calculated over the entire depth of the epilimnion, which is the distance between the surface (0 m) and the position of the thermocline (h Moreover, it is well known that eddy covariance data underestimate turbulent heat fluxes (Foken, 2008). We implemented an alternative approach to adjust the turbulent heat fluxes according to the annual EBR ratio: the corrected turbulent heat fluxes were obtained by redistributing the missing energy over a so-called 'energy year', namely from March 1 to February 28/29 of the next year. Since energy storage is typically at its minimum value at this time of the year, this allowed us to discard that term and write R n = H + LE. By doing so, Equation (3) becomes 2.6 | Ice cover Ice cover was monitored using a time-lapse camera (Reconyx, Holmen, WI, USA) that took hourly photos of the southern edge of the reservoir. Water temperatures were also used to identify the solar radiation penetration into the water column and ice formation, as ice and snow cover reduce or prevent radiation from entering the water column. More precisely, the completion of freezing takes place when the water surface is completely covered with ice. The end of the icecover period occurs when the reservoir becomes completely free of ice. In both cases, the identification of the dates was done by inspect-  Condensation episodes (LE <0 W m À2 ) occurred occasionally throughout the study period. During winters, small and short condensation episodes occurred almost every day. However, during the vernal turnover, conditions became very stable, and the atmosphere was much warmer than the water surface, which was nearly constant at 4 C. Consequently, only a few condensation events were recorded between May 1 and July 17, with cumulative amounts of 2.2 mm, 0.6 mm and 3.4 mm in 2019, 2020 and 2021, respectively.

Figure 7 presents examples of condensation events that occurred in
June, temporarily reaching À10 W m À2 or 0.35 mm day À1 . Although condensation was highest in June, it remained low compared to evaporation, which easily exceeded 100 W m À2 throughout the September-December period.
The reservoir was exposed to many episodes of sustained evaporation, defined as consecutive 24-h periods with a daily mean LE ≥ 100 W m À2 (i.e., daily mean LE ≥3.5 mm), ranging from one to several days. Of those, the most modest episode caused 3.5 mm of evaporation in one day, while the largest episode caused 29.3 mm of evaporation in 6 days (≈5 mm day À1 ).   Moreover, during these three evaporation events, the winds were from the SSW, indicating that they were channelled into the main axis of the reservoir. In addition, during these evaporation events, the air was drier than usual with a specific humidity below 5.5 g kg À1 , which was within the 30th quantile over the months of August and September. These episodes were each time due to the occurrence of a high-pressure system in conjunction with the passage of a cold front. The anticyclones brought drier air masses from the north and prevailing clear-sky conditions promoted incoming solar radiation on the water surface, while cold fronts generated strong winds after their transit. The combined effect of these phenomena enhanced the evaporative process. September at 78 W m À2 when the cumulative heat storage reached its maximum (refer to Figure 10) and remained greater than 60 W m À2 all the way to December. 84% of the annual total evaporation took place in the last 5 months of the year. However, in our study, LE decreased gradually due to the declining vapour pressure deficit that followed the decline in air temperature. In winter and spring, LE stayed below 20 W m À2 and reached its minimum in June, during the ice-free period ( Figure 5).

| Monthly scale
Sensible heat fluxes were negative from February to July, with values less than À10 W m À2 and reaching as low as À20 W m À2 when the reservoir was colder than the air. On average, values remained positive from August to January. In June, the reservoir was ice-free while the water column was under vernal transition, with upper layers that were not warm enough to stratify (i.e., below 3.98 C, the temperature of maximum water density). The water surface temperature stayed far below the air temperature, preventing H from becoming positive and enabling LE to become negative. In summer and fall, the surface layers were warmer than the air above the water. H increased steadily until December, reaching 80 W m À2 , and then decreased abruptly in January as the water surface froze.
From May to August, net radiation (R n ) was high, contributing mostly to the storage of heat in the reservoir (ΔH S ), while turbulent heat fluxes were relatively low. Nearly all the energy brought in by R n was used to increase the heat storage term. During the fall and early winter, R n declined rapidly while LE and H increased, fuelled by the energy stored in the reservoir, firmly establishing the heat release period of the reservoir. We observed a 3-month delay between the maximum summer net radiation in June and the maximum latent heat flux in September. The delay was 6 months between that peak of net radiation and the maximum sensible heat flux in December. We also observed different delays between the maximum surface water temperature and the maximum LE and H, which were delayed by 1 and 4 months, respectively.
Heat storage started in May and ended in early September (a 3.5-month period): the heat content of the water body rose by 130 W m À2 , 220 W m À2 and 175 W m À2 over the months of May, June and July, respectively. Heat release followed and lasted about 4 months until ice-on. The heat release rate was lower from September to October (À73 W m À2 ) compared to November to December (À250 W m À2 ). Overall, heat storage exhibited larger dayto-day fluctuations than heat release. and ice-cover periods. Figure 11 explores the relationship between daily LE and atmospheric stability ζ (ζ ¼ z=L ob , where L ob is the Obukhov length) for heat storage and release. It confirms that larger daily evaporation occurs under near-neutral (ζ ≈ 0) and unstable conditions (ζ < 0), and that stable conditions (ζ > 0) are related to low evaporation rates. Note that condensation occurred primarily under near-neutral and stable conditions. This is consistent with the findings from Blanken et al. (2011) for Lake Superior, where 89% of the annual evaporation occurred when there was unstable atmospheric stratification. It is also consistent with the Rouse et al. (2003) study at Great Slave Lake, where 85% to 90% of evaporation occurred during a period when the atmosphere was unstable. When looking at the monthly layer energy balance of the reservoir using Equation (3), we noted that there was a significant nonclosure (see Section 3.3). This is due to the time scale difference between the heat storage in the water column and the remaining terms (turbulent heat fluxes and net radiation). Indeed, due to the high thermal inertia and the high specific heat of water, the net radiation that was received in summer accumulated in the water column and was only released in autumn. This feature makes it difficult to calculate the EBR on a monthly scale. Within this context, correcting turbulent heat fluxes to account for the energy imbalance (Foken, 2008) using the monthly EBR is problematic unless accurate measurements of lateral energy inputs are available. This challenge is developed in Section 3.3 with monthly EBR from 2019 to 2021.

| Annual scale
Calculating annual reservoir evaporation values is challenging. The underestimation of turbulent heat fluxes by the eddy-covariance technic implies that the observed annual values must be corrected using the EBR over an 'energy year' (see Section 2.5). Since R n > H + LE, the missing energy can then be redistributed by preserving the observed Bowen ratio, as discussed in Mauder et al. (2018). By doing so, the EBR (Equation (5)) was 80%, 69%, and 76% for the years 2019-2020, 2020-2021, and 2021-2022, respectively.  Table 1). Note that the inter-annual variabilities of cumulative evaporation were 7.2% and 18% for non-corrected and corrected values, respectively. In fact, 50% of the total measured evaporation occurred over 24%, 27% and 26% of the days for 2019-2020, 2020-2021 and 2021-2022, respectively. This is consistent with Blanken et al. (2000), who found a mean value of 22.5% for the Great Slave Lake.
Few studies have focused on annual scales of lake or reservoir evaporation in the boreal biome. In Canada, Strachan et al. (2016) reported that the evaporation rate at the Eastmain-1 reservoir was 595 mm year À1 , reaching 100 mm month À1 in summer and Horizontal bars at the end of the energy period represent corrected annual evaporation, compensating for the lack of closure of the energy balance based on the annual EBR (see Table 1).
T A B L E 1 Annual characteristics from 2018 to 2021. Non-corrected and corrected total of evaporation, ice cover period, turbine outflow volume, temperature and wind speed anomalies compared to 1991-2020 (ERA5 values) and total net radiation. 3.1 mm day À1 from August to October. Rouse et al. (2003) found the annual evaporation at Great Slave Lake to be between 384 mm and 506 mm year À1 . For Lake Tämnaren in Sweden, Heikinheimo et al. (1999) reported 281 mm of evaporation from May to October. Finally, Blanken et al. (2011) found the annual evaporation at Lake Superior to be up to 645 mm. However, it is important to note that these studies were all based on eddy-covariance measurements and did not take into consideration the non-closure of the energy balance, meaning that they did not correct for cumulative evaporation values. This appears to be one major shortcoming in the annual estimation of turbulent heat fluxes. When examining the inter-annual variation of cumulative evaporation, we noted that the seasonal onset and the break-up of ice cover did not have a large impact on total yearly evaporative losses. The longer ice cover period in 2019 (145 days) than that of 2020 (133 days) resulted in lower cumulative evaporation. This trend is not repeated in 2021, which featured the shortest ice cover period (119 days) and lower cumulative evaporation compared to 2020 ( Figure 12). We also noted that the timing of ice onset was very consistent, as it began around January 1 of each year, while ice break-up was more variable between years. As for duration, the onset of ice lasted between 2 and 3 weeks while ice break-up lasted 5 to 7 weeks prior to the reservoir becoming completely ice-free.

| Uncertainties
There are uncertainties that affected the EC measurements that were used in this study. Random sampling uncertainties were calculated following Finkelstein and Sims (2001) and amounted to ≈2% for both sensible and latent heat fluxes over the whole study period. In addition, the error associated with gap-filling was noteworthy for H, LE and R n because of the difficulty of measuring fluxes over an inland water body. The use of a raft results in oscillations that must be considered in the calculation process, while the use of a flux tower on the shore limits the measurement period over the reservoir.
EC measurements were also subject to underestimations linked to the lack of energy balance closure (Foken, 2008), which is best described in terms of the energy balance ratio (EBR, Equation (3)). non-closure. The large overestimation of the January EBR is due to the very small denominator during this period. Indeed, both net radiation and heat storage are negative and close together at this time of year, allowing the January EBR to reach very high values ($150).
Some uncertainties lie in the evaluation of ΔH S . The monthly EBR was calculated using a ΔH S for depths of 0-15 m, 0-30 m or 0-44 m, which correspond to the average depths of the thermocline for January to October, November, and December, respectively. The heat content of water layers between the surface and these depths represent good approximations of the heat storage within the water column. This is because the thermal mixing is spatially constrained by the thermocline. However, there is uncertainty in estimating these thermocline depths, which may lead to errors in the evaluation of ΔH S .
Note that for 2019, there were no data for depths below 15 m.
When we calculated ΔH S for a 30-m and a 44-m depth, the EBR for November and December 2020 was 101% and 138% respectively.
This was less than the 150% and 440% calculated for the water between the surface and 15-m depth.

| CONCLUSION
This study quantified the temporal dynamics of evaporation over a deep, boreal, dimictic hydroelectric reservoir using two eddycovariance setups, one mounted on a raft and one onshore. Data were Results showed an annual evaporation of 590 ± 66 mm year À1 that was quite constant from year to year, with frequent 1-day to 2-day sustained events. Latent heat flux increased earlier than the sensible heat flux but also decreased before the sensible heat flux, resulting in a Bowen ratio that varied from a near-zero negative value in July to 1.5 in December. Vapour pressure-controlled evaporation induced a steady decline from September to December due to decreasing air temperature.
The large time lag and the magnitude of the energy storage within the water column made it difficult to close the energy balance. Therefore, in this study, we have taken that issue into account by correcting the annual cumulative evaporation while preserving the measured Bowen ratio.
Moreover, monthly, and seasonal patterns of evaporation can be related to the energy state of the reservoir. Indeed, depending on the time of year, the reservoir was either under ice cover or in heat storage or heat release conditions, which drives the magnitude of evaporation. Finally, because evaporation is likely to increase in the region due to climate change, the assessment of this energy and the associated hydraulic components remains topical and essential to understanding future trends.