A Simplistic Approach for Assessing Hydroclimatic Vulnerability of Lakes and Reservoirs with Regulated Superﬁcial Outﬂow

: This study proposes a simplistic model for assessing the hydroclimatic vulnerability of lakes / reservoirs (LRs) that preserve their steady-state conditions based on regulated superﬁcial discharge ( Q d ) out of the LR drainage basin. The model is a modiﬁcation of the Bracht-Flyr et al. method that was initially proposed for natural lakes in closed basins with no superﬁcial discharge outside the basin ( Q d = 0) and under water-limited environmental conditions {mean annual ratio of potential / reference evapotranspiration ( ET o ) versus rainfall ( P ) greater than 1}. In the proposed modiﬁed approach, an additional Q d function is included. The modiﬁed model is applied using as a case study the Oreastiada Lake, which is located inside the Kastoria basin in Greece. Six years of observed data of P , ET o , Q d , and lake topography were used to calibrate the modiﬁed model based on the current conditions. The calibrated model was also used to assess the future lake conditions based on the future climatic projections (mean conditions of 2061-2080) derived by 19 general circulation models (GCMs) for three cases of climate change (three cases of Representative Concentration Pathways: RCP2.6, RCP4.5 and RCP8.5). The modiﬁed method can be used as a diagnostic tool in water-limited environments for analyzing the superﬁcial discharge changes of LRs under di ﬀ erent climatic conditions and to support the design of new management strategies for mitigating the impact of climate change on (a) ﬂooding conditions, (b) hydroelectric production, (c) irrigation / industrial / domestic use and (d) minimum ecological ﬂows to downstream rivers.


Introduction
Lakes and reservoirs are important assets due to the high value of ecosystem services that they provide [1][2][3][4]. If the worst predictions of general circulation models (GCMs) about climate change become true, then lakes/reservoirs (LRs) will hardly manage to maintain their current conditions not only in terms of quantity but also of quality [5][6][7][8][9]. Under such conditions, the local authorities will face great challenges for setting strict measures able to preserve and restore the conditions of LRs [10][11][12][13] especially in regions which have been identified as climate change hot spots such as the Mediterranean Basin [14].
The hydrological models are important tools for assessing the water balance components that regulate the water volume of LRs and can assist the design of water management strategies. Depending on the modelling purposes and the specific cases, different types of models and different levels of model complexities can be selected, starting from extremely data demanding and complex models such as MIKE SHE [15,16], GSFLOW [17], WATLAC [18], SMS [19], QGIS-SWAT [20], or using simpler methods such as those provided by Bracht-Flyr et al. [21], Yang et al. [22] and others.
A simple method for analyzing the effect of climate conditions on LRs is the simplistic hydroclimatic lake classification model provided by Bracht-Flyr et al. [21]. The method is applied for water-limited environments (mean annual potential evapotranspiration greater than precipitation) under steady-state conditions and takes into account the mean annual precipitation, mean annual potential evapotranspiration and an empirical factor that describes the effect of the ratio between the terrestrial area of the basin that drains to the lake A B versus the area of the lake A L (A ratio = A B /A L ). The method also assumes that the total lake basin A T (A T = A B +A L ) follows the Budyko hypothesis [23,24], which postulates that the mean annual actual evapotranspiration of A T asymptotically approaches the mean annual precipitation P as the climate gets drier. A basic advantage of the method is that has low data requirements, but its main limitation is that it can be used only for closed lake basins. This is a basic problem since in many cases LRs present regulated superficial discharge of water outside the drainage basin, and for this reason the aforementioned method requires adaptations.
The aim of the study is to propose a modification of Bracht-Flyr et al. [21] hydroclimatic model that considers regulated superficial discharge of water outside the drainage basin of the LR. The new modified model was applied for a case study in Greece (Orestiada Lake in Kastoria basin) for analyzing the current conditions, while it was also used for analyzing the hydroclimatic vulnerability of the lake under climate change considering the predictions of 19 general circulation models for three climate change scenarios RCP2.6, RCP4.5 and RCP8.5.

Study Site
The study area is the Orestiada Lake of the Kastoria hydrologic basin in northwestern Greece (40 • 30 45 N; 21 • 18 10 E) (Figure 1a). The maximum allowed lake surface elevation (LSE) is at 628.19 m above sea level (a.s.l.) with a maximum lake depth of~9.1 m. Values above the maximum LSE can flood the residential areas around the lake, especially the city of Kastoria, which is at the eastern coast of the lake. For this reason, at the southern coast of the lake, water is discharged through a sluice gate to Gioli channel ( Figure 1a). The bed and the upper edge of the sluice gate is at 624.95 and 627.45 m (a.s.l.), respectively. The end of Gioli channel is the final outlet of Kastoria basin and it is connected to Aliakmon river basin ( Figure 1a) [25,26]. The local lake management authorities follow a program of water discharge in order to maintain a secure lake surface elevation at~627 m (a.s.l.), which is accomplished by opening for a few days per month (~2 days per month on average) the sluice gate with a mean approximate rate of discharge~2-3 m 3 /s. The water discharge from the lake satisfies also other purposes such as: • improvement of lake water quality by using the inflow of cooler and denser water for flushing out the warmer lake water [27][28][29]. In this way, a reduction of primary production rates is accomplished, which is necessary because the lake is currently classified as eutrophic with a tendency to hypertrophication [27,30,31]. • partial water contribution to Aliakmon river ( Figure 1a) for preserving its minimum ecological flow and for the production of hydropower by hydroelectric dams at its downstream sections.
The spatial variation of mean annual precipitation and temperature of Kastoria basin is given in Figure 1b,c, respectively [32]. The basin climate is classified as Cfb for the lowlands and Dfb for the uplands at the northern area according to the revised Köppen Geiger classification system of Peel et al. [33]. The coverage of land uses in the study area according to Corine Land Cover 2012 [34] database is approximately 36.4% agricultural lands (AGRI), 17.5% forests (FOR), 9.7% water bodies (WAT), 34.1% scrub and/or herbaceous vegetation associations (SHVEG) and 2.3% artificial surfaces and urban fabric (ARTS) (Figure 1d). The coverage of soil texture classes in the study area consists of Loamy Sand (3.9%), Sandy Loam (47.3%), (Loam 9.6%), Sandy Clay Loam (4.7%), (Clay 34.9%) according to the ESDB-ESDAC soil database [35,36] with sand increasing with altitude and the clayey soils dominating to the lowlands areas around the lake.
Hydrology 2019, 6, x 3 of 24 according to the ESDB-ESDAC soil database [35,36] with sand increasing with altitude and the clayey soils dominating to the lowlands areas around the lake.

Data
The data used for the implementation of this study are the daily minimum, maximum, mean temperature and precipitation data from the Kastoria meteorological station (632 m, 40°30'38"N; 21°15'15"E) for the period 1/5/2012-31/4/2018 (Figure 2a,bFigure 2a; Figure 2b) and the variation of lake surface elevation (LSE) (Figure 2c). This period corresponds to 6 full years with starting date the 1st of May for each year. For the specific period, the mean annual superficial discharge Qd of the lake to Gioli channel through the sluice gate was estimated at 5.6 million m 3 /year (based on records of the Environmental and Spatial Planning Agency of Western Macedonia Prefecture, Department of Environment and Hydroeconomy).

Data
The data used for the implementation of this study are the daily minimum, maximum, mean temperature and precipitation data from the Kastoria meteorological station (632 m, 40 • 30 38 N; 21 • 15 15 E) for the period 1/5/2012-31/4/2018 (Figure 2a,b) and the variation of lake surface elevation (LSE) (Figure 2c). This period corresponds to 6 full years with starting date the 1st of May for each year. For the specific period, the mean annual superficial discharge Q d of the lake to Gioli channel through the sluice gate was estimated at 5.6 million m 3 /year (based on records of the Environmental and Spatial Planning Agency of Western Macedonia Prefecture, Department of Environment and Hydroeconomy). Additional GIS datasets and databases used for this study include: • The digital elevation model of the Kastoria basin obtained by the EU-DEM v1.1 database [37] with a spatial resolution of 25 m (Figure 1a). • The bathymetry map of Orestiada lake designed for reference LSE = 628.00 m a.s.l. with isolines of 0.5 m depth [38] (Figure 1e). • The average monthly climate data for minimum, mean, maximum temperature and precipitation of the recent past (1970-2000, WorlClim Version 2 database) and the respective parameters based on 19 general circulation models (GCMs) for three scenarios (RCP2.6, RCP4.5 and RCP8.5 according to IPCC/CMIP5) of future climate conditions (mean conditions of 2061-2080 period) derived from the WorldClim Version 1 database [39] ( Table 1). The datasets are in raster form with spatial resolution of 30 arc-sec (~1 km 2 ). The WorldClim database also provides the results of RCP6.0 scenario from the respective GCMs, but it was not used because it showed very small differences with the RCP4.5 for the study area. • The local revised coefficients of Hargreaves and Samani equation [40] provided by Aschonitis et al. [41] were also used for achieving equivalent estimations of reference crop evapotranspiration, ETo, with the complete formula of ASCE/FAO-56 method for short grass [42,43]. The revised coefficients have been produced based on the WorldClim database and they are provided in raster format with 30 arc-sec (~1 km 2 ) spatial resolution [41]. Additional GIS datasets and databases used for this study include: • The digital elevation model of the Kastoria basin obtained by the EU-DEM v1.1 database [37] with a spatial resolution of 25 m (Figure 1a).

•
The average monthly climate data for minimum, mean, maximum temperature and precipitation of the recent past (1970-2000, WorlClim Version 2 database) and the respective parameters based on 19 general circulation models (GCMs) for three scenarios (RCP2.6, RCP4.5 and RCP8.5 according to IPCC/CMIP5) of future climate conditions (mean conditions of 2061-2080 period) derived from the WorldClim Version 1 database [39] ( Table 1). The datasets are in raster form with spatial resolution of 30 arc-sec (~1 km 2 ). The WorldClim database also provides the results of RCP6.0 scenario from the respective GCMs, but it was not used because it showed very small differences with the RCP4.5 for the study area.

•
The local revised coefficients of Hargreaves and Samani equation [40] provided by Aschonitis et al. [41] were also used for achieving equivalent estimations of reference crop evapotranspiration, ET o , with the complete formula of ASCE/FAO-56 method for short grass [42,43].
The revised coefficients have been produced based on the WorldClim database and they are provided in raster format with 30 arc-sec (~1 km 2 ) spatial resolution [41].  [21] considers a simplified water balance for lakes assuming steady-state conditions. According to this model, the general mean annual water balance of a lake for a multiyear period is described by: where ∆V L is the net annual change in the lake volume (L 3 ), A L is the lake surface area (L 2 ), P L is the mean annual precipitation over the lake (L), E L is the mean annual evaporation from the lake (L), R B is the mean annual runoff/recharge generated by the terrestrial area of the basin that drains to the lake (L), A B is the contributing basin area that drains to the lake (L 2 ), and G in , G out are the mean annual incoming and outgoing groundwater fluxes to and from the lake (L 3 ), respectively. Assuming that an equilibrium condition is attained between all the inflows and outflows from the lake volume (∆V L ≈0) and that a net groundwater transfer is negligible (G in ≈G out ), Equation (1) is simplified according to the following: where A is the ratio of basin surface that drains to the lake versus the lake surface. The R B parameter is described by the following equation [21,44]: where E pB and P B are the mean annual potential evapotranspiration and precipitation (L) of the basin surface that drains to the lake and ω 1 is an empirical factor that represents the effects of soil and land use that both regulate real evapotranspiration. f 1 (Φ,ω 1 ) is an asymptotic function that ranges between 0-1. Bracht-Flyr et al. [21] made the assumption of uniform climate conditions in the whole basin of the lake and provided the following simplified form: where E p and P are the mean annual potential evapotranspiration and precipitation (L) of the whole basin (including the lake), respectively. Equation (4) works only for water-limited conditions (Φ≥1) and not for energy limited ones (Φ < 1). In the second case, the formation of a lake under a stationary climate would be permanent and solely controlled by topography for maintaining constant dimensions [21]. For steady-state conditions, it is easy to estimate the ω 1 parameter using available climatic data of the basin and using the current A B /A L ratio. The function f 1 (Equation (3c)) is a Budyko-type function [23,24] and describes the term ET a /P, where ET a is the mean annual actual evapotranspiration [21,44,45], which is given by the following equation: The ω 1 value can vary depending on the soil type, land use and the method used for estimating potential evapotranspiration [44][45][46]. Zhang et al. [46] noted that ω 1 is a lumped factor that represents the integrated effect of multiple basin processes and land cover conditions on evapotranspiration.
Assuming that the basin features (e.g., land uses, soil properties, etc.) that regulate ω 1 value remain constant (i.e., ω 1 remains constant in the future), then the changes in the term A using climatic data derived from GCMs for different future climate projections can be assessed.
The approach of Bracht-Flyr et al. [21] is interesting and constitutes a useful tool for types of analysis as the aforementioned one, but there are some aspects in the methodology that need to be clarified. For example, the initial assumption that was made for model building was that the system is closed (Equation (1)) while for building the final empirical model of Equations (4) and (5), other additional assumptions were made (e.g., uniform evapotranspiration in the whole A T , Budyko-type of real evapotranspiration described by Equation (5). If we consider that there are no other factors that extract water from the A T basin and that the initial assumption of closed system is still valid, then the final water balance under steady-state conditions for the A T area should simply be described as P-ET a = 0. In this case, and under a water limited environment (Φ>1), the Budyko-type function of ET a (Equation (5)) requires a value of ω 1 →+∞ for achieving P-ET a ≈0. Thus, in a theoretical steady-state and water-limited environment (Φ>1) of a closed lake basin with 100% natural landscape and water losses only by evapotranspiration, it is expected ω 1 →+∞. The above indicate that a lake basin system where the ET a is described by Equation (5) with a regular ω 1 value (e.g., <10) [45,46], cannot support the initial assumption of a closed system because P-ET a will be always >0, which consequently means that there are additional water extraction factors for preserving the system under steady state conditions. Based on the above, if we consider that the initial assumption of "closed system" is violated because (a) it is more probable other factors of water extraction exist and (b) the Budyko-type of real evapotranspiration is valid based on a regular ω 1 value (i.e., not reaching +∞), then the final water balance under steady conditions in the A T basin should be described by: where Q m is the annual amount of water extraction (L) from the system of A T basin by factors not related to ET a (e.g., groundwater recharge outside the A T , net industrial/domestic water consumption inside the A T without return flows etc). Violating the initial assumption of closed system affects also the ω 1 value because the larger Q m is, the smaller ω 1 will be.

Case 2: Outflows Outside the Drainage Basin of the Lake (Modified Method)
In many LRs, as in the case of Orestiada Lake, steady-state conditions may be preserved because water is superficially discharged outside the A T area (e.g., typical case of a reservoir with hydroelectric dam, where water is discharged for energy production or a reservoir used for irrigating downstream areas). In such cases, Equation (4) can be modified as follows using an additional factor f 2 according to the following: where where Q d is the mean annual amount of water (L) that superficially discharges outside the A T area for maintaining steady state conditions of the LR, and ω 2 is the coefficient that regulates the rate of superficial discharge. The f 2 (Φ,ω 2 ) is also an asymptotic function that ranges between 0-1 but with inverted shape compared to f 1 (Φ,ω 1 ), which denotes that the higher is the ET a and Φ, the lower will be the Q d . If the mean annual values of Q d , E p , P and A are known for a relatively large period of years then Equation (8a) can be used for estimating ω 2 and Equation (7) for estimating ω 1 using ω 2 .
We have to note that any abstraction or water outflow from the lake that returns to A T (e.g., using lake water for irrigation inside A T ) is not considered in the Q d . Moreover, an additional modification was made by substituting the term of mean annual potential evapotranspiration E p with the term of mean annual reference evapotranspiration ET o in Φ calculation (i.e., Φ = ET o /P) of Equations (7) and (8). The ET o is calculated by the formula of Hargreaves and Samani [40]: where ET o,d is the reference evapotranspiration (mm d −1 ), R a is the daily extraterrestrial radiation (MJ m −2 ), λ is the latent heat of vaporization (2.45 MJ kg −1 ), T mean is the mean daily temperature estimated as the average value of maximum and minimum daily temperature ( • C), TD is the difference between maximum and minimum daily temperature ( • C) and c hs2 is a recalibrated coefficient that corrects the original Hargreaves-Samani equation in order to be equivalent to ASCE/FAO-56 method for short reference crop [42,43]. The original value of c hs2 is 0.0023, while the local recalibrated c hs2 values for any place of the world can be obtained by Aschonitis et al. [41].
Using the climate conditions of a relatively large period of years where the LR is under steady-state conditions due to superficial outflows outside the A T (meaning that Q d , ET o , P and A are known parameters), then ω 2 is estimated after re-arranging the parameters of Equation (8a) as follows: Using the estimated ω 2 from Equation (10b), ω 1 can also be estimated after re-arranging the parameters of Equation (7) as follows: where Equations (10) and (11) are valid for Φ = ET o /P>1, P 0, X 1 and Y 1. The rule Y 1 denotes that Q d is always 0. For Q d = 0 the original Bracht-Flyr et al. [21] method should be used.
Using the same assumptions of Section 2.3.1. for deriving Equation (6), the final water balance under steady conditions in the A T basin for Q d 0 is described by: where Q m is the annual amount of water extraction (L) from the system of A T basin by other factors not related to Q d and ET a (e.g., groundwater recharge outside the A T , net industrial/domestic water consumption inside the A T without return flows etc).

Estimating LR Conditions under the Effects of Climate Change
The modified method (Equations (7)- (12)) allows evaluating the LR conditions under climate change considering that a) any changes in future soil and land use management or other basin characteristics do not alter the ω 1 value and b) Q d is always 0 and the Q d rates are defined by the future P and Φ considering a stable shape of f 2 function (i.e., ω 2 remains constant). For example, in the case of Orestiada Lake, it is assumed that outflow continues to exist with Q d rates regulated by Equstion 8a despite the fact that LSE may fall below the bed elevation of the sluice gate. Under such conditions, the final A B /A L ratio is equal to A (Equation (7)) considering that lake water is pumped and guided to Gioli channel in order to preserve the outflow for improving lake water quality, the minimum ecological flow and the hydropower production of Aliakmon river as described in the "Study site" section (pumping water from Orestiada Lake to Gioli channel could be a realistic case if we consider a similar effect of climate change on the whole Aliakmon river basin). The aforementioned strategy is denoted as "strategy S1" and the results of A lake volume V L , superficial outflow Q d , real evapotranspiration ET a and other extractions Q m of the final steady state conditions are symbolized as A S1 , V L , S1, Q d,S1 , ET a,S1 , and Q m,S1 , respectively.
Taking into account the results of the aforementioned procedure, any deviation from the current climate conditions could indicate what the final expected steady-state conditions (A S1 , V L , S1, Q d,S1 , ET a,S1 , and Q m,S1 ) will be even before they occur. This allows designing alternative strategies for the volume preservation of LR based on the estimated Q d,S1 under climate change. For example, the estimated Q d,S1 can be returned to LR (in reality Q d,S1 is just calculated but it is not discharged) for restoring LR's volume at any time, even before V L,S1 reaches its final steady-state condition. If we set a minimum threshold value of lake volume (V L,th ) for which S1 strategy should stop and the estimated Q d,S1 should return to/remain in LR (denoted as Strategy S2) then we can estimate a respective minimum time in years that is required for restoring LR volume to the initial conditions as follows: where t S2 is the minimum time (years) that it is required for restoring LR volume to V L,0 from V L,th condition, V L,0 is the LR volume at the initial/current conditions (L 3 ), V L,th is the minimum threshold value of lake volume (L 3 ) for which S1 strategy should stop and Q d,S1 should return to/remain in the lake. Q d,S1 is expressed in volume units (L 3 ) in Equation (13a). The term "minimum time" for t S2 is used because at the initiation of S2 strategy, the lake volume is not at steady state conditions and still declines until reaching V L,S1 (case of Equation (13a)). Thus, the annual Q d,S1 are added to an initial volume which continues declining. Thus, the real time for LR volume recovery is expected to be longer. In this study, as V L,th was set the lake volume that corresponds to LSE at the bed of the sluice gate while V L,0 corresponds to the mean volume of Orestiada Lake during the period 1/5/2012-31/4/2018 (i.e., current conditions).

Linking A Parameter to Lake Water Volume
The lake area, A L , and the lake water volume, V L , can be estimated for different lake surface elevations when the bathymetry of LR is known. A B and A parameter (A B /A L ratio) can also be estimated for each LSE value when the total area A T (A T = A B +A L ) of the basin is known. Thus, each value of LSE can be linked to specific values of A L , A B , A and V L .
In this study, the isolines of bathymetry of Orestiada Lake were used to develop a triangulated irregular network (TIN) of the lake bed. The TIN was used to determine the lake surface A L and volume V L for different values of LSE using the tool {SurfaceVolume} from 3D Analyst toolbox of Arc-GIS 10.1 (ESRI). The Arc-Hydro module of ArcSWAT model [47] was also used to delineate the total area of Kastoria basin that drains to Orestiada lake including the lake (i.e., A T ). The known value of A T and the values of A L for different LSEs were used to calculate the respective A B and A values. In this way, the A values were linked to V L for different LSEs.

Using Meteorological Stations as Descriptors of the Whole Basin Climate
There are various methods (e.g., Thiessen polygons, interpolation methods, temperature, rainfall gradients, etc.) for deriving the mean annual P and ET o for basins that have a good number and distribution of meteorological stations. On the other hand, there are basins that may not have an adequate number of stations and especially basins with large portions of non-accessible mountainous areas where meteorological stations usually do not exist. In such cases, the derivation of the mean annual P and ET o of the basin from a few stations may not be representative for the whole basin. One solution to this problem is to use gridded climatic data that have been produced by GCMs and advanced interpolation methods that take into account topographic characteristics for achieving finer grid resolutions (e.g., such gridded data are those provided by WorldClim database, e.g., Figure 1b,c). The average pixels' value of the WorldClim gridded mean monthly data of precipitation, maximum and minimum temperature of recent past  and future climate scenarios (Table 1) were extracted from the polygon that corresponds to A T . The respective values were also extracted from the specific pixel that corresponds to the meteorological station of Kastoria (Figure 1a). Thus, twelve mean monthly values of precipitation, maximum and minimum temperature from the grids of the 51 simulation scenarios of GCMs (Table 1) plus the case of the (1970-2000) (52 cases × 12 months for each parameter) were extracted for the position of the station and for the whole A T basin. These values were used through regression analysis to develop three linear functions (g 1 , g 2 , g 3 ) that calculate the mean monthly values of T max,T , T min,T and P T of the drainage basin (defined by A T ) using the respective parameters from the position of Kastoria station, respectively: where P S , T max,S , T min,S are the mean monthly precipitation, maximum and minimum temperature at the position of the Kastoria meteorological station, and P T , T max,T , T min,T are the mean monthly precipitation, maximum and minimum temperature of the A T basin. The calibrated Equation (14a)-(14c) were finally used to convert the observed data of monthly P S , T max,S , T min,S of the station for the period 1/5/2012-31/4/2018 to the respective mean values of the A T basin. These values were further used to calculate ET o with Equation (9) and Φ of A T for the respective period. We have to note that mean temperatures are calculated as the average of minimum and maximum temperature in all cases.

Bathymetry, Lake Volume, and A Parameter for Different Lake Surface Elevations (LSEs)
The TIN surface of Orestiada lake bed was created taking into account the 0.5 m isolines of bathymetry (Figure 1e). The TIN was used to determine the lake surface A L and volume V L for different values of LSE (Table 2). Table 2. Lake surface A L , lake volume V L , basin area that drains to the lake A B , A B /A L and mean depth of the lake for different values of lake surface elevation (LSE) considering A T = 272.48 km 2 (grey zone plus the lake in Figure 3).

Estimating ω1 and ω2 for the Current Period (1/5/2012-31/4/2018)
The WorldClim gridded mean monthly data of precipitation, maximum and minimum temperature of recent past  and future climate scenarios (Table 1) Table 2). Using the mean annual The Arc-Hydro module delineated the total area of Kastoria basin (308.6 km 2 ) and created 26 subbasins from which 21 drain to the lake while 5 drain directly to Gioli channel at the southern part of the lake towards the outlet of the Kastoria basin ( Figure 3). The total area A T of the 21 subbasins that drain to the lake (including the area of the lake) was estimated A T = 272.48 km 2 . Thus, using the value of A T , the values of A B and A were estimated based on the respective values of A L for different LSEs (Table 2).

Estimating ω 1 and ω 2 for the Current Period (1/5/2012-31/4/2018)
The WorldClim gridded mean monthly data of precipitation, maximum and minimum temperature of recent past  and future climate scenarios (Table 1) were extracted (a) from the total basin of A T = 272.48 km 2 (mean value of pixels in the basin) and (b) from the specific pixel that corresponds to the meteorological station of Kastoria (Figure 1a). Thus, 52 × 12 pairs of (T max,T , T max,S ), (T min,T , T min,S ) and (P T , P S ) were developed in order to calibrate the functions g 1 , g 2 and g 3 (Equations (14a)-(14c)) through linear regression. The g 1 , g 2 and g 3 functions were the following T max,T = 0.99T max,S -1.128 with R 2 = 0.99, T min,T = 0.988T min,S -1.063 with R 2 = 0.99 and P T = 0.928P S + 3.515 with R 2 = 0.99.
The monthly T max,T , T min,T and P T of the A T basin for the study period 1/5/2012-31/4/2018 were estimated based on the calibrated g 1 , g 2 and g 3 functions and using the respective monthly observed data of T max,S , T min,S and P S of Kastoria station. The monthly T max,T and T min,T of the A T basin were also used to estimate the monthly Hargreaves-Samani ET o values (Equation (9)) of the basin using c hs2 = 0.002087 (average c hs2 value of pixels included in A T ) from the database of Aschonitis et al. [41]. Finally, the mean annual P and ET o of the A T for the 2012-2018 period were estimated at 701 and 968 mm per year, respectively. The mean observed superficial discharge to Gioli channel of the respective period was 5.6 million m 3 /year, which corresponds to Q d = 20.6 mm per year. The mean LSE of the study period was 626.99 m, which corresponds to A = 9.06 and V L = 84.61 million m 3 (estimated based on linear regressions between sequential values of Table 2). Using the mean annual values of P = 701 mm, ET o = 968 mm, Φ = ET o /P = 1.381, Q d = 20.6 mm and A = 9.06, the ω 2 value was found equal to 16.56 according to Equation (10c). Using this ω 2 value, the ω 1 value was found equal to 6.09 according to Equation (11c). Using Equation (5), the mean annual real evapotranspiration ET a of the studied period was found equal to 651 mm while Q m was estimated equal to 29.4 mm per year according to Equation (12b) (equivalent to 8.1 million m 3 per year).

Trends of Climate in the AT of Orestiada Lake according to GCMs for Different Climate Change Scenarios
The box-whisker plots of mean annual values of P, T, Tmax, Tmin, TD, ETo and Φ = ETo/P of GCMs for each climate change scenario (RCP2.6, RCP4.5, RCP8.5) are given in Figure 5a Table 3). The individual values for each GCM and climatic scenario that were used for creating the box-whisker plots of Figure 5a-g are given in the Supplementary Materials (Table S1).
The worst-case scenario of climate change was the one provided by the GFDL-CM3 (no.6 code GF in Table 1) (Equation (3c)) and f 2 (Equation (8b)) values versus Φ using the ω 1 and ω 2 of this study (the other f 1 cases are provided for comparative purposes and they have been derived from [45,46]).

Trends of Climate in the A T of Orestiada Lake according to GCMs for Different Climate Change Scenarios
The box-whisker plots of mean annual values of P, T, T max , T min , TD, ET o and Φ = ET o /P of GCMs for each climate change scenario (RCP2.6, RCP4.5, RCP8.5) are given in Figure 5a Table 3). The individual values for each GCM and climatic scenario that were used for creating the box-whisker plots of Figure 5a-g are given in the Supplementary Materials (Table S1).
The worst-case scenario of climate change was the one provided by the GFDL-CM3 (no.6 code GF in Table 1)

Effects of S1 and S2 Strategy on Orestiada Lake Conditions under Different Climate Change Scenarios
The final mean annual A S1 and Q d,S1 of Orestiada lake and the ET a,S1 and Q m,S1 of A T basin under the effects of S1 strategy and climate change were estimated by Equations (5), (7), (8a) and (12b), respectively, taking into account ω 1 = 6.09 and ω 2 = 16.56, and the future conditions defined by the GCMs of each future climate change scenario (Figure 5a-g, Table S1). The Q d,S1 (mm) and Q m,S1 (mm) values were converted to (m 3 × 10 6 per year) considering A T = 272.48 km 2 . The final A S1 values were also expressed as lake volumes V L,S1 (m 3 × 10 6 ) using Table 2 (through linear regression between 0.5 m intervals). The box-whisker plots of A S1 , V L,S1 , Q d,S1 , ET a,S1 and Q m,S1 estimations according to S1 strategy for the climate projections of GCMs based on the three RCP scenarios are given in Figure 6a-e, respectively. For comparison purposes, the current mean annual values of the aforementioned parameters for the period 1/5/2012-31/4/2018 are also given as red lines in Figure 6a-e. The mean trends of V L,S1 , Q d,S1 , ET a,S1 and Q m,S1 of GCMs for each climate change scenario (red crosses in Figure 6b-e) were compared with the mean annual conditions of the period 1/5/2012-31/4/2018 (red lines) in order to evaluate the % changes from current conditions ( Table 3). The individual values for each GCM and climatic scenario that were used for creating the box-whisker plots of Figure 6a-e are given in the Supplementary Materials (Table S2). Table 3. Maximum negative, average, standard deviation and maximum positive % of changes in mean annual values of P, T, T max , T min , TD, ET o , Φ, V L,S1 , Q d,S1, ET a,S1 and Q m,S1 estimated based on the GCMs for RCP2.6, RCP4.5, RCP8.5 in comparison to the respective mean annual values of the period 1/5/2012-31/4/2018. As regards the S2 strategy, the minimum time t S2 for restoring the volume of LR was estimated taking into account that the mean annual lake volume for the current conditions of the period 1/5/2012-31/4/2018 is V L,0 = 84.61 m 3 × 10 6 and that the lake volume that corresponds to the bed elevation of the sluice gate (LSE = 624.95 m a.s.l.) is V L,th = 33.2 m 3 × 10 6 . The box-whisker plots of t S2 values for each RCP scenario are given in Figure 6f and Table S2.
Some interesting observations from the analysis of S1 and S2 strategies are the following: • there were 2 GCM cases (no.9 HD for RCP2.6 and no.5 CN for RCP4.5) out of the 51 (all cases of the three scenarios) where the final V L,S1 was greater than mean annual value V L,0 of the period 1/5/2012-31/4/2018 (i.e., wetter future climatic conditions compared to the current ones). Thus, for these two cases, the t S2 is 0 according to Equation (13b) for S2 strategy.
• 53.3% of GCM cases for RCP2.6, 94.7% of GCM cases for RCP4.5 and 100% of GCM cases for RCP8.5 showed V L,S1 <V L,th (i.e., final lake volumes smaller than the volume that corresponds to the bed elevation of the sluice gate).

Violation of the Assumption for Negligible Net Groundwater Transfer or Other Types of Outflows
Both the original (Equation 4) and the modified method (Equation 7) are based on the simplistic formula of Equation 2, which was built considering a negligible net groundwater transfer (Gin≈Gout). In both cases, this assumption was finally violated in order to close the water budget using Equations 6 and 12, respectively, considering additional terms for describing outflows outside the lake basin. This element could be considered to be a problem in the theoretical background but in reality it is not because the concept behind the formation of both models is (a) first to derive a simplistic formula as a base (i.e., Equation 2) for building the models and (b) to introduce the terms f1 and f2 that convert Equation 2 to two semi-empirical models (Equation 4 or Equation 7), which in Figure 6. Box-whisker plots of mean annual values of (a) A B /A L = A S1 , (b) V L,S1 (m 3 × 10 6 ), (c) Q d,S1 (m 3 × 10 6 ), (d) ET a,S1 (mm), (e) Q m,S1 (m 3 × 10 6 ) according to S1 strategy, and (f) values of t S2 (years) according to S2 strategy. The red lines correspond to the mean annual values of the period 1/5/2012-31/4/2018.

Violation of the Assumption for Negligible Net Groundwater Transfer or Other Types of Outflows
Both the original (Equation (4)) and the modified method (Equation (7)) are based on the simplistic formula of Equation (2), which was built considering a negligible net groundwater transfer (G in ≈G out ). In both cases, this assumption was finally violated in order to close the water budget using Equations (6) and (12), respectively, considering additional terms for describing outflows outside the lake basin. This element could be considered to be a problem in the theoretical background but in reality it is not because the concept behind the formation of both models is (a) first to derive a simplistic formula as a base (i.e., Equation (2)) for building the models and (b) to introduce the terms f 1 and f 2 that convert Equation (2) to two semi-empirical models (Equation (4) or Equation (7)), which in reality are non-linear functions where their coefficients are calibrated. Thus, the two final semi-empirical models (Equations (4) and (7)) indirectly include the groundwater or other types of inflows-outflows (if they exist) because they are calibrated with actual data that may include such effects.

Validity of ω 1 in the Budyko Hypothesis and Justification of Q m based on Local Factors
The validity of the final estimated steady-state water balance (Equations (6a) and (12a)) requires evaluation of the individual parameters (P, ET a , Q d , Q m ) based on expert judgement (especially the case of ET a and Q m parameters). The accuracy of mean P in A T can be evaluated based on the availability, spatial distribution and accuracy of means that measure meteorological parameters, while Q d purely depends on the authorities that manage the superficial outflows. On the other hand, ET a and Q m parameters are the parameters with the higher uncertainty and can be evaluated through indirect observations.
As regards ET a , the estimated ω 1 = 6.09 based on the current conditions (1/5/2012-31/4/2018) may seem quite high but it is not unrealistic since similar values have been already observed in other studies [45,46]. Taking into account Figure 4, it is observed that for Φ = 1.38 (mean annual Φ of 1/5/2012-31/4/2018 for the study area), the f 1 for a sandy soil with sparse vegetation is f 1 (ω 1 = 0.1) = 0.61, for a grass-covered basin is f 1 (ω 1 = 2.55) = 0.86, for a forest-covered basin is f 1 (ω 1 = 2.84) = 0.87 while for f 1 (ω 1 = 6.09) = 0.928. These values indicate that for the same amount of P, the conditions of our study area would provide 34.3% higher ET a from a respective basin covered with sandy soil and sparse vegetation, 7.3% higher ET a from a respective mostly grass-covered basin, and 6.3% higher ET a from a respective mainly forest-covered basin. These differences can be considered realistic and justifiable due to the following combined reasons: • the lake covers~10% of the A T (Figure 3) where evaporation is non-limited and approximates ET o [42]. Thus, the high percentage coverage of the lake in the A T could be a basic factor that leads to a high value of ω 1 . • another reason of the high ω 1 value can also be the high percent (~30%) of agricultural land in the A T area (Figure 1d). Summer irrigated crops and especially those in water limited environments (Φ>1) receive irrigation water for achieving maximum evapotranspiration rates during the summer period of low rainfall. In A T area, the % of irrigated land is~9.8% considering the agriculture statistics of Kastoria, Vitsi, Makethnon and Agioi Anargyroi municipalities [48], which are inside A T . The rest winter non-irrigated crops (~20.2%) increase the water availability for their own evapotranspiration demands because they act as cover crops that can reduce runoff especially in areas of higher slopes [49,50]. Thus, the contribution of agricultural lands and especially of irrigated ones in the A T , could also be a basic factor that leads to a higher value of ω 1 .
Considering the above, the estimated ET a rates by Equation (5) using ω 1 = 6.09 have a realistic base at least for the current conditions. As concern the future conditions, which predict smaller lake areas/volumes due to higher evapotranspiration and lower rainfall, we can assume that any possible reduction of ω 1 that could occur due the reduction of lake evaporative surface can be counterbalanced by the higher lake evaporation rates due to the shallower and warmer waters and by the higher evapotranspiration rates of irrigated agricultural lands due to the increase of temperature. The general reduction of ET a of the whole A T (Figure 5d) due to ET o increase and P decrease can be justified by the fact that the available water for evapotranspiration will be reduced for the natural non-aquatic landscape and non-irrigated agricultural land that cover more than 80% of the A T area.
The domestic water use and groundwater outflow are the two main factors that regulate the existence and rate of Q m parameter. The local domestic water consumption is supplied from both ground and surface waters, while all the domestic wastewater goes to the wastewater treatment plant (WWTP) of Kastoria, which is near the Gioli channel. All the treated water is discharged to Gioli channel, thus the return flow to A T by domestic water consumption is almost null. According to UNEP/MAP [51], the WWTP of Kastoria treated 5640 m 3 /d (equivalent to annual amount of 2.058 m 3 × 10 6 or 7.55 mm based on A T ). Considering the above, it is understood that a minimum amount of 7.55 mm per year is extracted from A T system through domestic water consumption due to null return flows. Considering that the mean annual Q m of the 2012-2018 period was 8.01 m 3 × 10 6 (or 29.4 mm) according to Equation (12b), indicates that the domestic water consumption is approximately 25.7% of the current Q m . The rest part of Q m (5.95 m 3 × 10 6 or 21.85 mm) is assumed to be loss of water through groundwater flow from the A T area (lake and grey areas in Figure 3) towards to the drainage areas of Gioli channel (red areas in Figure 3) and consequently to the watershed of Aliakmon river. If we consider a constant domestic water consumption in the future then all the future estimated changes in Q m,S1 (Figure 6e) are associated with changes in groundwater outflows from the A T system.

Practical Meaning of S1 and S2 Strategies
It is possible to argue that there are practical uncertainties regarding the assumptions behind the S1 and S2 strategies. For example in S1, there is always outflow unless the lake dries up, which may indicate a non-sustainable approach for preserving the LR. These two strategies have three practical uses: (a) In some LRs, management protocols require the maintenance of minimum outflows, perhaps to support irrigation, electric production, and preservation of ecological flows for protecting important downstream aquatic habitats (e.g., deltas) that host valuable species diversity. Stopping the outflow may also lead to more severe implications to the downstream environments. For example, if the outflow of LR supplies a river that goes to the sea, then stopping the upstream outflow may lead to disturbance of the transition/mixing zone (brackish water) at the outlet/delta, allowing sea water to intrude in the river and in coastal aquifers [52,53]. Thus, the S1 scenario can be used as a decision support and management tool that allows providing estimations of lake dimensions and outflows that can be used to evaluate the sustainability of the whole system (lake basin and downstream areas). (b) The S1 scenario can also be used even when the local conditions allow cessation of the surface discharge for preserving the LR. In this case, S1 is used to evaluate lake dimensions and surface discharge based on specific climate conditions and then the results are used in S2 scenario. In reality, S2 is a mathematical trick where the surface discharge is estimated by S1 but never goes out of the system either because the surface discharge stops or because an equivalent amount of water returns to the lake. The second case is very important since it can support management plans for preserving the LR and its downstream outflows by returning equivalent amounts of water by another source (e.g., by diverting a river of another basin to lake basin) [54,55]. (c) In S2 strategy, t S2 (Equation (13) has an alternative absolute physical meaning, which is the time needed to save water equal to the amount V L,0 −V L,th by setting Q d = 0.

Conclusions
The overall analysis and the provided results of the modified Bracht-Flyr et al. (2013) methodology showed that it can provide (a) a simple evaluation of the mean annual A B /A L ratio of an LR with regulated superficial outflow and (b) a water balance of the A T area that hosts the LR, considering steady state conditions based on the mean climatic conditions of the current studied period but also for the future climate projections of GCMs. The results of the modified method showed that it can be used as a diagnostic tool for analyzing the response of LRs with superficial regulated discharges and to support the design of new LR management strategies for mitigating the impact of climate change. The use of the modified method is proposed for LRs that discharge water outside their drainage basin for supporting necessary purposes such as flooding prevention, hydroelectric production, irrigation/industrial/domestic use and preservation of minimum ecological flows of downstream rivers. Future studies could focus on improvements of the original and modified methodologies in order to be applicable for energy-limited environments (Φ < 1), while in the case of the modified methodology (where ω 1 is calibrated based on ω 2 for Q d 0) a procedure for direct predictions of A under steady-state conditions based on different climatic conditions assuming Q d = 0 is also required.  The supplementary materials contains an excel file with Tables S1 and S2. Author Contributions: Conceptualization, all authors; methodology, V.A., and K.D.; formal analysis, K.D., and D.P.; resources, D.P. and V.A.; data curation, K.D., and D.P.; writing-original draft preparation, all authors.; writing-review and editing, all authors; supervision, D.P. and D.P.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations Abbreviation Description
Unit A B Terrestrial area of the basin that drains to the lake L 2 A L Area of the lake L 2 A T Total area of the basin including lake area L 2 A Ratio between the terrestrial area of the basin that drains to the lake versus the area of the lake unitless ratio A S1 Ratio between the terrestrial area of the basin that drains to the lake versus the area of the lake according to S1 strategy for a future climate change condition Mean annual amount of regulated superficial discharge out of the drainage basin of a lake or reservoir for maintaining steady state conditions according to S1 strategy for a future climate change condition L Q m Mean annual amount of water extraction from the whole basin by factors not related to actual evapotranspiration and not regulated superficial discharge out of the drainage basin L Q m,S1 Mean annual amount of water extraction from the whole basin by factors not related to actual evapotranspiration and not regulated superficial discharge out of the drainage basin according to S1 strategy for a future climate change condition L R 2 Coefficient of determination unitless R a Daily extraterrestrial radiation MJ/m 2 R B Mean annual runoff/recharge generated by the terrestrial area of the basin that drains to the lake L RCPs Representative concentration pathways of greenhouse gass concentration trajectory -S1 Strategy of lake management where the regulated superficial outflow always exist with a rate which is defined considering constant ω 1 and ω 2 values regardless of climate and lake conditions (except if lake volume becomes 0) -S2 Strategy of lake management where the outflow Q d,S1 should return to the lake when a minimum threshold value of lake volume is reached The minimum time that it is required for restoring LR volume to V L,0 from V L,th condition when S2 strategy is applied years V L Lake volume L 3 V L,0 Initial/current conditions of lake volume L 3 V L,S1 Lake volume according to S1 strategy for a future climate change condition L 3 V L,th Minimum threshold value of lake volume for which S1 strategy should stop and the outflow Q d,S1 should return to the lake L 3 WAT Water bodies in Corine Land Cover map -∆V L Net annual change in the lake volume L 3 λ Latent heat of vaporization MJ/kg Φ Ratio of potential/reference evapotranspiration versus precipitation unitless ratio ω 1 Empirical factor that represents the effects of soil and land use that both regulate real evapotranspiration unitless ω 2 Empirical factor that regulates the rate of regulated superficial discharge out of the drainage basin unitless