Assessing Lake Response to Extreme Climate Change Using the Coupled MIKE SHE/MIKE 11 Model: Case Study of Lake Zazari in Greece

Abstract

Lakes, either artificial or natural, are greatly important as a component in their catchments’ hydrology, but also as ecosystem service providers. However, due to climate change, they have begun to face numerous problems with their water quality and quantity. Furthermore, general circulation models (GCMs) show future climate change with a reduction in rainfall and increase in temperature. The aim of the current study is to present an application where GCMs and state-of-the-art hydrological modelling system MIKE SHE/MIKE 11 are combined for assessing the response of a Greek lake in terms of its water balance and water level under climate change. Four general circulation models (GCMs; GFDL-CM3, MIROC-ESM-CHEM, MIROC-ESM, IPSL-CM5A-LR) for the extreme climate change scenario of RCP8.5 were used in the basin of Lake Zazari in Greece as a case study. Results showed that, by keeping the irrigated demands (the main water user) unchanged in the future, the lake exhibited a lower water level for all GCMs, fluctuating from −0.70 to −1.8 m for the mean (min) water level and from −0.30 to −1.20 m for the mean (max) water level. Instead of the above and by preserving the amount of withdraw water n from the lake at a certain percentage of inflows, the irrigated area should be reduced from 54.1% to 64.05% depending on the circulation model.

1. Introduction

Natural and artificial lakes are greatly important due to their great contribution to the supply of high-value ecosystem services [1]. The impact of climate change and human activities on their water quality and quantity is thoroughly investigated [2,3,4,5,6,7,8,9]. Among the main reasons is to support management actions related to their ecosystem protection, water supply, flood, and water shortage phenomena [10]. Such actions are expected to face great challenges for setting appropriate measures to preserve or restore lake conditions [11,12,13,14], especially in identified climate change hot spots such as the Mediterranean basin [15] where Greece is located, with more than 50 lakes or reservoirs in its territory.

A solution for developing management strategies is the application of hydrological models with future climate data based on predictions of general circulation models (GCMs). In this way, the water balance components of lakes and their response under extreme climate conditions can be investigated. Depending on the purposes of modelling and the specific attributes of a lake and its catchment, there are different types of hydrological model structures and levels of model complexities that can be selected. For example, there are more complex and data-demanding models, such as MIKE SHE/MIKE 11 [16], WASP [17], CE-QUAL-W2 [18], GSFLOW [19], WATLAC [20], PCLake [21], DYRESM-CAEDYM [22], or simpler ones [6,10,23,24]. The MIKE SHE modelling system is among the most integrated and it was used for analyzing climate change impact on water resources [25], and more specifically on lakes [26,27]. The MIKE SHE model consists of various compartments for different hydrologic analysis purposes. The most integrated combination is its coupling with MIKE 11, which can be used to perform detailed hydrologic and hydraulic analysis at catchment scale applications based on the hydraulic characteristics and dimensions of rivers and lakes [28].

As far as the Greek territory is concerned, there are few studies that have been conducted for wetlands and specifically for natural lakes. There are studies based on either statistical or mathematical models, but these refer to streams that flow into lakes [29] or analyze the impact of the whole hydrological response of the basin into lakes [30,31]. For example, the authors in [32] studied Lake Karla, one of the most interesting lakes in Greece that was drained in 1962 and then restored by the authorities in 1999 due to its environmental importance but also its vital role in irrigated agriculture. The purpose of the study was the water resources management of the basin. In that direction, the authors in [33] studied the role of adaptation measures on agriculture to water resources under climate change, while the authors in [34] studied the socioeconomic impact that a catchment can have with a lake. There is no clear approach that focuses on lakes apart from an effort conducted in [35] to simulate Lakes Volvi and Koronia using hydrological models such as MIKE SHE and MIKE Hydro. However, apart from the Greek territory, scientists around the word are studying the whole response of a catchment that contains a lake or a reservoir in climate change and not the lake itself and the lake’s behavior [36,37,38].

The current study investigates Zazari Lake, situated in northern Greece, due to its environmental importance. There are few studies that focused on the origin of Zazari. The authors in [39] investigated the effect of using annual updated land use maps derived from satellite images instead of a static coarse land use map in a SWAT model for the estimation of sediments and nitrates that load in Zazari Lake. In [40], the authors investigated the internal trends of surface temperature, precipitation, and wind speed over 18 lakes of Greece using ERA5-Land monthly averaged data from 1981–2020. Moreover, the authors in [29] investigated the impacts of climate change on the hydrology of Lakes Cheimaditida and Kerkini by using monthly data and a conceptual model (MIKE BASIN-NAM). The basin of Zazari Lake is part of the wider basin of Cheimaditida Lake. Lastly, the authors in [41] studied the mean monthly water balance of the wider basin of Cheimaditida Lake using the Thorthwaite and Turc methods for actual evapotranspiration estimation.

The aim of this study is to present an application for assessing the response of a lake to climate change using the coupled MIKE SHE/MIKE 11 modelling system in the catchment of Lake Zazari in Greece, thereby creating a fully appropriate system for wetland hydrological applications. Although the coupled MIKE SHE/MIKE 11 model is not new, the novelty of this work is the fact that this coupled MIKE SHE/MIKE 11 model approach has not been used so far for investigating changes in the lake levels due to climate change in a daily time step in Greek lakes. The study was performed for Lake Zazari in Greece considering the four worst climate change predictions of four respective general circulation models (GCMs) for the climate change scenario of RCP8.5. Inflows and outflows from the lake, and its water level fluctuation were used to assess climate change pressure on the lake’s hydrology.

2. Materials and Methods

2.1. Study Area

The study area is the hydrologic basin of Zazari Lake (40°37′30″ N, 21°32′50″ E; northwestern Greece; Figure 1), which covers an area of 98.4 km² and belongs to the River Basin District of Western Macedonia–Greece. According to the updated world map of the Köppen–Geiger climate classification of [42], the region belongs to the Csa climatic zone (hot summer Mediterranean climate). Annual precipitation during the five-year period of the study ranged from 600 to 800 mm, and mean temperature from 4 °C during winter, and 22.6 °C during summer.

Τhe mountainous area of the basin is mainly composed of low-permeability rocks of a crystalline background, while alluvial deposits cover its lower part. In these quaternary sediments of about 120–130 m thickness [43], two main hydraulically independent aquifer systems were formed, a surface shallow one with aquifer thickness of up to 45 m that grows in coarse-grained surface formations, and a deeper confined aquifer. The shallow aquifer is hydraulically connected with the lake [44,45].

The lake is supplied with water from both its own subcatchment and the Sclithros torrent. The latter, before outflowing into the lake, receives water from the Asprogia and Aetos streams (Figure 1a). Aetos stream, just after it leaves the mountainous area, becomes a drainage ditch of nonpermanent flow.

A sharp crested weir of 17 m width that is located at the outlet of the lake regulates its surface water level. Overspill starts at 599.4 m above sea level (a.s.l.). At this elevation, the lake covers an area of ~1.95 km² with a maximal depth of 6 m.

The main water user in the area is agriculture. The lake supplies with water the irrigation network of the Limnohori area (220 ha, 1.19 × 10⁶ m³ mean yearly abstraction from May to September), which is located out of its basin, and 60 more hectares of individual plots close to its shore. From the agricultural area that is located close to the settlement of Sclithros, 364 ha are irrigated from Sclithros River, Asprogia stream, and pumping wells. Moreover, 480 ha of agricultural land that is located south of the Aetos and Agrapidies settlements are irrigated by groundwater.

2.2. Hydrological Modeling System

In the context of this study, the coupled MIKE SHE/MIKE 11 modelling system was used. MIKE SHE is a physical, deterministic, fully distributed hydrological model [16,46] that can simulate the major processes of water cycles on both surface and subsurface [47]. Evapotranspiration, infiltration, overland flow, unsaturated and saturated flow, irrigation and drainage, and overland–groundwater exchange were simulated in this study using the corresponding MIKE SHE modules.

MIKE 11 is a 1-dimensional hydraulic modelling system for rivers, lakes, and estuaries that is based on the complete dynamic wave formulation of the Saint Venant equations. Furthermore, it can represent a wide range of hydraulic structures including sluice gates, weirs, gates, bridges, and culverts [28,48].

The two systems, MIKE SHE and MIKE 11, were dynamically coupled [49,50], exchanging information after each computational time step [47,48]. The main streams in the basin and Lake Zazari were included in MIKE 11 for the purposes of this study.

2.3. Data

In the context of this study, MIKE SHE/MIKE 11 was calibrated on the basis of the conditions of the period of 25 July 2012–31 December 2016, and validated on the conditions of the period of 1 January 2017–31 December 2017. The observed data were the variation in lake surface elevation (LSE), with maximal and minimal water levels of 599.98 and 598.56 m, respectively, and the daily temperature and precipitation data from the Limnochori meteorological station (located nearby the lake, 40°37′20.70″ N, 21°33′19.27″ E) for the aforementioned period. Reference evapotranspiration (ET_o) was calculated on a daily time step by the equation in [51], taking into consideration the local value from the global map of revised coefficients provided by [52] for achieving equivalent estimations of FAO-56/ASCE reference crop evapotranspiration ET_o for short grass (Figure 2) [53,54]. On the basis of ET_o, MIKE SHE estimates the actual evapotranspiration from land surface by taking into account the LAI and RD of the land use, the available ponded water, and the soil moisture content in each cell. In the case of flooded areas such as lakes, MIKE SHE assumes that actual evaporation equals to ET_o, unless a different approach is adopted [47].

Τhe digital elevation model for the study area was obtained by the EU-DEM v1.1 database with a spatial resolution of 25 m, in which lake bathymetry was incorporated [55].

The land uses of the basin according to the Corine Land Cover 2012 database [56] consist of 32.6% croplands, 46.8% forests, 8.7% pastures/grasslands, 8.3% shrubs/sclerophyllous vegetation, 2% lake, 1.3% artificial surfaces, and 0.3% bare soils (Figure 1b). The general soil properties obtained from the ESDB-ESDAC soil database [57,58] showed that the study area is dominated by loamy sand and sandy loam soils at 36.9% and 44.2% coverage, respectively.

For the analysis of climate change, the predictions of four GCMs {GFDL-CM3 (gf), MIROC-ESM-CHEM (mi), MIROC-ESM (mr), IPSL-CM5A-LR (ip)} for the worst scenario of the highest greenhouse-gas emissions RCP8.5 (mean conditions of 2061–2080) were used. These 4 GCMs were selected among 19 GCMs provided by the WorldClim database [59] because they represented the worst conditions in terms of rainfall reduction and temperature increase compared to the current conditions at the position of the lake. The data are in raster form with spatial resolution of 30 arc-s (1 km²).

2.4. Model Setup and Parameterization

The model area included the whole basin of the lake. Two grid sizes, 400 × 400 m and 200 × 200 m were tested [60] for its discretization. The 400 × 400 m size did not allow for an accurate representation of the river network and lake shape, so it was abandoned before any calibration attempt. Instead, the grid size of 200 × 200 m gave a finer and more accurate representation of the river and lake shape, and of the spatial distribution of physical parameters in the basin. This discretization resulted in 2709 cells (70 rows × 85 columns) and was adopted.

The overall maximal time step was set at 6 h, with maximal time step for the sheet overland flow, 0.5 h; water flow in rivers, 30 s, in the unsaturated zone, 2 h, and 24 h in the saturated zone.

Land use spatial extent was further corrected, taking into account Google Earth images. Leaf area index (LAI) and root depth distribution (RDD) over time were acquired from [61] for each class of CLC 2012.

For the simulation of the irrigation in the basin, eight irrigation command areas were adopted (areas that receive water from the same source). Their spatial extent and water sources are given in Figure 3. In all command areas, irrigation starts at wilting point and ends at field capacity during the period for which LAI and root depth are greater than zero. The irrigation method was sprinklers in all cases.

The hydraulic properties of soil types (Figure 4) were estimated using Van Genuchten equations according to the dataset of [62] (

Table 1. Soil type, saturated moisture content (θ_s, cm³/cm³), residual moisture content (θ_r, cm³/cm³), suction at field capacity (pF_fc), suction at wilting point (pF_w), saturated hydraulic conductivity (K_s, m/s), empirical constants in Van Genuchten formulas a (cm⁻¹) and n (-), used in the hydrological model of Lake Zazari basin.

Parameter	Clay	Loam	Silt	Loamy Sand	Sandy Clay Loam	Sandy Loam
θs	0.38	0.43	0.46	0.41	0.39	0.41
θr	0.068	0.078	0.034	0.057	0.1	0.065
pFfc 1	2.8	2.2	2.4	1.9	2.5	2
pFw 1	4.2	4.2	4.2	4.2	4.2	4.2
Ks	5.55 × 10−7	2.89 × 10−6	6.944 × 10−7	4.05 × 10−5	3.64 × 10−6	1.23 × 10−5
a	0.008	0.036	0.016	0.124	0.059	0.075
n	1.09	1.56	1.37	2.28	1.48	1.89

¹ pF: log cm of H₂O column.

). The full Richard’s equation was used with the computations to be carried out in all grid sells (soil columns) of the model. For the functional relationship of moisture–retention curve and effective conductivity, the formulas of Van Genuchten were selected using default values for their parameters (Table 1). The unsaturated zone was considered to be uniform to the vertical direction, and its discretization was chosen to vary from 10 cm for the first 0.5 m (5 nodes), 50 cm from 0.5 to 9.0 m (17 nodes) depth, and 200 cm from 9.0 to 45 m depth (18 nodes).

According to [44,45], the saturated zone was simulated by adopting a single surface aquifer seated on an impermeable layer that extended across the catchment. In addition, the aquifer is extended beneath the lake, allowing for a dynamic exchange of water between them. Both the depth of the aquifer and hydraulic conductivities were subjected to calibration. Initial values were set to 5 m depth at the mountainous area, and 40 m at the lower part of the catchment. Following [63], almost 90% of the aquifers under exploitation in Greece consist mainly of sand and gravel with hydraulic conductivities between 10⁻⁶ and greater than 10⁻⁴ m/s. Since no further information regarding this parameter was found for the study area, the initial value of 5 × 10⁻⁵ m/s was adopted for both horizontal (K_x) and vertical (K_y) hydraulic conductivity.

The outer boundary condition at the periphery of the aquifer was set as zero flux, except for a section east of the catchment. For the latter, a permanent hydraulic gradient of −0.0025 was adopted that followed the natural slope of the ground surface. This boundary condition drains ground water downstream out of the catchment.

The dense network of streams in the mountainous areas and the artificial drainage network in the agricultural land were simulated by using the drainage module of MIKE SHE. The drain level was set at −1.5 m below ground surface for the whole catchment, while the time constant (leakage factor) was set to be 1 × 10⁻⁷ for the whole catchment except for the valleys, where it was set to be 1 × 10⁻⁶. For drainage routing, grid codes were adopted in order for the drain flow to be routed to the nearest river link. The subcatchment of the lake was excluded from the drainage (drain codes set to 0) due to the absence of significant streams.

For the simulation of water flow in Sklithros River, and the Asprogia and Aetos streams (Figure 1a), MIKE 11 was used. Only the sections of the above water bodies located in the lowlands of the catchment were included in the model (Figure 5). In this way, steep slopes of the mountainous branches and accompanied instabilities for the model were avoided. Lake Zazari was simulated as part of Sklithros River by inserting several cross-sections following the shape of its shores (Figure 5). The whole length of the 11 MIKE branches was coupled with MIKE SHE.

Cross-sections for each branch were specified from field investigations. The Manning coefficient for all branches was set to be n = 0.04 [64], except for the lake section, which was set to be n = 0.03. Leakage coefficients were initially set to be 1 × 10⁻⁷, but they were subjected to calibration for all branches.

As far as the initial conditions of the model are concerned, two conditions were set at the beginning of the simulation, the first was water depth along the branches that was set to be to 0.01 m, and the second was the initial water depth of the lake, which was set to be equal to the observed water level at the beginning of the simulation.

A small discharge of 0.05 m³/s for the first few days that became zero for the rest of the simulation period was set to be the boundary condition at the upper free ends of the branches. At the outlet of the lake, a Q/h boundary condition was set representing the functioning of the crested weir (Figure 6). Further, the water abstraction from the lake for the supply of the Limnohori collective irrigation network was set as a discharge boundary condition.

2.5. Criteria for Evaluation of Lake Zazari Basin Hydrological Model

Model performance was evaluated by comparing the simulated lake surface elevation (LSE) against its observed measurements, since it was the only parameter that is monitored in the catchment almost at a monthly time step since 2012 [65]. The statistical criteria that were used, the range of their expected values, and values that indicated the best fitting of the model results are described in

Table 2. Statistical criteria for evaluating the performance of the Lake Zazari basin hydrological model.

Criterion	Equation	Range	Best Value
Mean error	$M{E}_j=\frac{{{\displaystyle \sum }}_t\left(Ob{s}_{it}-Si{m}_{it}\right)}{n}$	−∞ ~ +∞	0
Mean absolute error	$MA{E}_j=\frac{{{\displaystyle \sum }}_t\left.\left|Ob{s}_{i,t}\right.-Si{m}_{i,t}\right|}{n}$	0 ~ +∞	0
Root mean square error	$RMS{E}_j=\sqrt{\frac{{{\displaystyle \sum }}_t{\left(Ob{s}_{i,t}-Si{m}_{i,t}\right)}^2}{n}}$	0 ~ +∞	0
Standard deviation of residuals	$STDre{s}_j=\sqrt{\frac{{{\displaystyle \sum }}_t{\left(\left(Ob{s}_{i,t}-Si{m}_{i,t}\right)-\left(\frac{{{\displaystyle \sum }}_t\left(Ob{s}_{i,t}-Si{m}_{i,t}\right)}{n}\right)\right)}^2}{n}}$	0 ~ +∞	0
Correlation coefficient	$R_j=\frac{{{\displaystyle \sum }}_t\left(Si{m}_{i,t}-\overline{Si{m}_{i,t}}\right)·\left(Ob{s}_{i,t}-\overline{Ob{s}_{i,t}}\right)}{\sqrt{{{\displaystyle \sum }}_t{\left(Si{m}_{i,t}-\overline{Si{m}_{i,t}}\right)}^2·{{\displaystyle \sum }}_t{\left(Ob{s}_{i,t}-\overline{Ob{s}_{i,t}}\right)}^2}}$	0 ~ 1	1
Nash–Sutcliffe correlation coefficient	$R_j^2=1-\frac{{{\displaystyle \sum }}_t{\left(Ob{s}_{i,t}-Si{m}_{i,t}\right)}^2}{{{\displaystyle \sum }}_t{\left(Ob{s}_{i,t}-\overline{Ob{s}_{i,t}}\right)}^2}$	−∞ ~ 1.0	1

Obs: observed values at a location i and at a time t, Sim: simulated values at a location i and at a time t, n: number of observations, $\overline{Ob{s}_{i,t}}$: means of observations at location i, $\overline{Si{m}_{i,t}}$: means of calculations at location i.

[47,66].

2.6. Climate Change Scenarios

For the analysis of climate change, the predictions of four GCMs {GFDL-CM3 (gf), MIROC-ESM-CHEM (mi), MIROC-ESM (mr), IPSL-CM5A-LR (ip)} for the worst scenario of highest greenhouse-gas emissions RCP8.5 (2061–2080) were used. These 4 GCMs were selected among 19 GCMs provided by the WorldClim database [59] because they represented the worst conditions in terms of precipitation reduction and temperature increase compared to the current conditions at the position of the lake. As the coupled MIKE SHE/MIKE 11 needs daily data for the simulation, the following was assumed in order for the model to be able to simulate in future time: historical data series taken from the weather station of Limnochori were assumed to have the same future pattern. So, the average values of meteorological (temperature) or aggregated (precipitation) parameters in the historical series for each calculated month and the respective values of each month from the Fick and Hijmans dataset were used and compared, and the percentage of increase or decrease of each parameter was calculated. Then, future time series were generated by multiplying each value of the historical dataset with the percentage of increase or decrease that had been calculated in the previous step. After that, the coupled models were again applied using the four new climate datasets that represented the four worst future climate conditions in order to assess the response of the lake.

3. Results and Discussion

3.1. Model Calibration and Validation

The calibration and validation of the model were conducted against observed measurements of the lake surface elevation (LSE) for the period of 25 July 2012–31 December 2016 and 1 January–31 December 2017, respectively.

Simulated vs. observed LSE data for both of the above-mentioned phases are given in Figure 7, which shows that the model satisfactorily followed the seasonal pattern of the lake’s water level fluctuation.

The values of the statistical criteria described in Table 2 were calculated and are given in

Table 3. Values of statistical criteria used for the evaluation of the Lake Zazari basin hydrological model performance.

	ME	MAE	RMSE	STDres	R	R2
Calibration 25 July 2012–31 December 2016	0.0193	0.0719	0.0942	0.0922	0.9720	0.9331
Validation 1 January–31 December 2017	0.0778	0.0982	0.1248	0.0975	0.9899	0.8733

, which shows that, in both the calibration and validation phases, ME and MAE were in the order of less than 10 cm when the LSE fluctuated between 598.56 and 599.98 m a.s.l. (in the range of 1.42 m). RMSE and STDres had values close to zero, which is the best value for these criteria. The correlation coefficient (R) and Nash–Sutcliffe model efficiency (R²), also suggested a very good predictive capability of the model, which is also reflected in the validation phase.

Table 4. Final values of calibrated parameters in the coupled MIKE SHE/MIKE 11 hydrological model of Lake Zazari basin.

Parameters	Mountainous Area	Lowlands	Lake
MIKE SHE
Aquifer depth (m)	12	15–25	6
Hydraulic conductivty, Kx, Ky (m/s)	1.5 × 10−5	2 × 10−4	1.5 × 10−5–2 × 10−4
MIKE 11
Leakage coefficient (s−1)		1 × 10−7–1 × 10−6

shows the final values of the parameters that were subjected to calibration. Hence, for the mountainous areas, the depth of the aquifer was 12 m, while for the lowlands, a depth from 15 to 25 m was adopted. Beneath the lake, the depth of the aquifer was 6 m. The final values for hydraulic conductivity after calibration were 1.5 × 10⁻⁵ m/s for the mountainous areas, 2 × 10⁻⁴ m/s for the lowlands, and 1.5 × 10⁻⁵–2 × 10⁻⁴ beneath the lake. The leakage coefficient for MIKE 11 branches was in the range of 1 × 10⁻⁷ to 1 × 10⁻⁶ s⁻¹. Lower values do not allow for the basin to drain, leading to a continuous increase in the water that is stored in the saturated zone, while greater values lead to the continuous drainage of the basin.

3.2. Response of Lake Zazari to Climate Change

For the assessment of Lake Zazari’s response to climate change, the hydrological model of the basin was rerun for each future GCM, and pressures on the lake’s water balance and state were quantified. Data of daily reference evapotranspiration (ET_o) and precipitation (P) that had been generated from each GCM were used in the hydrological model. The same initial conditions were adopted as those at the current period. In addition, the irrigated area in the basin and the abstraction of water to meet the needs of the Limnohori irrigation network were kept the same. In this way, four sets of results, one for each GCM, were produced for the future period.

Parameters that were used to quantify the pressure on the lake was precipitation (mm) and evaporation (ET_o, mm) from its surface, inflow from Sklithros (mean Q_in, m³) and outflow as overspill (mean Q_out, m³) as well as the mean minimal (mean (min)) and maximal (mean (max)) LSE (m), and are given in

Table 5. Comparison of Lake Zazari water balance components assessed by MIKE SHE/MIKE 11 under current climatic conditions and future climatic projections.

	Current Conditions	Future Predictions
	Observed	gf		ip		mi		mr
	mm	mm	Δ (%)	mm	Δ (%)	mm	Δ (%)	mm	Δ (%)
Precipitation	731	405	−44.7	475	−35.1	544	−25.6	520	−28.9
ETo	961	1136	18.1	1052	9.4	1095	13.9	1117	16.2
Discharge	Q (m3)	Q (m3)	Δ (%)	Q (m3)	Δ (%)	Q (m3)	Δ (%)	Q (m3)	Δ (%)
Mean (Qin)	5.73 × 106	2.06 × 106	−64.15	2.48 × 106	−56.70	2.63 × 106	−54.11	2.42 × 106	−57.75
Mean (Qout)	4.02 × 106	6.81 × 104	−98.30	2.60 × 105	−93.54	2.99 × 105	−92.55	2.11 × 105	−94.76
Water level	LSE (m)	LSE (m)	Δ (%)	LSE (m)	Δ (%)	LSE (m)	Δ (%)	LSE (m)	Δ (%)
Mean (max)	599.70	598.50	−0.20	599.20	−0.08	599.40	−0.05	599.10	−0.10
Mean (min)	598.70	596.90	−0.30	597.80	−0.15	598.00	−0.12	597.70	−0.17
Fluctuation	1.00	1.60	−60.00	1.40	−40.00	1.40	−40.00	1.40	−40.00

.

Inflow from Sklithros was estimated using the components of

Table 6. Water balance components of sub-basins (9116 ha) upstream the lake Zazari under current (obs) (1 January 2013–3 December 2017) and future (2061–2080) climatic conditions.

Water Balance Components	Future Climatic Projections
	obs (mm)	gf (mm)	ip (mm)	mi (mm)	mr (mm)
Precipitation	731.3	404.7	474.6	544.5	520.1
Evapotranspiration	593.2	432.0	460.3	513.3	510.6
Overland runoff to river (ORR)	2.7	0.2	0.3	0.7	0.3
Irrigation (IR)	31.3	40.7	40.3	41.3	40.3
Pumping from groundwater (PG)	29.5	38.9	38.5	39.1	38.4
Saturated zone drain to rivers (DR)	44.2	9.8	12.0	12.8	11.7
Baseflow to river (BR)	17.8	14.5	16.6	17.3	16.4
Baseflow from river (RB)	0.1	0.1	0.1	0.1	0.1
Sklithros River outflow, mm	63	23	27	28	26
Sklithros River outflow, m3	5.73 × 106	2.06 × 106	2.48 × 106	2.59 × 106	2.41 × 106

following the equation:

Q_in = ORR + DR + BR − (IR − PG) − RB

(1)

where ORR is the amount of overland runoff to the river (mm), DR is the amount of water from the saturated zone that drains to rivers (mm), BR is the amount of water from the baseflow that flows to the river (mm), IR is the amount needed for irrigation (mm), PG is the pumping amount of water abstracted from the groundwater (mm), and RB is the amount of water that flows from the river to the baseflow (mm).

Table 6 shows the components of the water balance of the area (sub-basins of Sklithros torrent, Asprogia and Aetos streams, 9116 ha) that contributes with runoff to torrent Sklithros and hence to Lake Zazari (Table 5). In Table 6, the irrigation component includes both pumping water from groundwater and extractions from surface water (Sklithros torrent and Asprogia stream).

The values of Table 5 and Table 6 represent mean yearly values of the five-year period of 1 January 2013–31 December 2017 and future conditions of the period of 2061–2080 of each climatic projection. Pressure was evaluated by comparing the values of the above parameters under the current period against future ones.

In Figure 8, the observed values of LSE for the period of 1 January 2013–31 December 2017 are compared against the estimation from the hydrological model values of LSE for the future period and for each GCM {GFDL-CM3 (gf), MIROC-ESM-CHEM (mi), MIROC-ESM (mr), IPSL-CM5A-LR (ip)}. Under all GCMs, the LSE shifted to lower levels, with that of the gf-model being the lowest. The mean (max) water level is expected to drop from −0.30 to −1.20 m, while the mean (min) from −0.70 to −1.80 m, according to each GCM.

Changes in the lake’s level fluctuation are due to changes caused in the components of its water balance, that is, evaporation from and precipitation on the lake’s surface, water inflows and outflows, and direct abstractions.

Evaporation from the lake is increasing, as is indicated by the values of ET_o (Table 5). Hence, greater amounts of water (from 18.1% to 9.4%) directly leave the lake.

In all GCMs, mean yearly precipitation substantially dropped, from 544 mm (−25.6%) under the mi-model up to 405 mm (−44.7%) under the gf-model, instead of 731 mm under current climatic conditions. This leads to less direct rainfall on the lake’s surface and reduced amount of water inflows into the lake as runoff from its basin.

Table 5 shows that the mean yearly discharge under the four GCMs from Sklithros River into Zazari drops by more than its half (−54.11% to −64.15%) compared to current climatic conditions. This can be seen from Table 6, which shows the water balance components of the Sklithros river. The lake’s water balance and thereby its level fluctuation are greatly affected, since the irrigated area was considered to be the same.

The reduction in lake level results in smaller quantities of water that overspill downstream. Figure 8 shows that, just after the second year of the future period, the LSE drops below the crest of the weir (599.4 m a.s.l.), and water stops overspilling downstream.

In the above future climatic projections, the area of the irrigation network of Limnohori (220 ha) was considered to be unchanged. Under current climatic conditions, 1.19 × 10⁶ m³ (5400 m³/ha) are directly withdrawn every year from the lake for the irrigation of this area. Further to the above, 60 ha more of individual plots located near the shores of the lake are irrigated directly from the lake (0.324 × 10⁶ m³). Hence, the total amount of water that is yearly abstracted is estimated to be Q_ir = 1.512 × 10⁶ m³ (Table 5). R_ir = 0.26 (or 26%) of the mean yearly inflows mean (Q_in) to the lake under current climatic conditions (Table 5 and

Table 7. Pressure analysis of irrigation under current climatic conditions and future climate projections.

	Current Conditions	Future Climatic Projections
	Observed	gf	ip	mi	mr
Qir (m3)	1,512,000.0	1,512,000.0	1,512,000.0	1,512,000.0	1,512,000.0
Mean (Qin) (m3)	5,730,000.0	2,060,000.0	2,480,000.0	2,630,000.0	2,420,000.0
Rir = Qir/Qin	0.26	0.73	0.61	0.57	0.62
Qir under stable Rir = 26%	1,512,000.0	543,581.2	654,408.4	693,989.5	638,575.9
Irrigated area (ha)	280.00	100.66	121.19	128.52	118.25
Reduced irrigated area (%)		64.05	56.72	54.10	57.77

). However, under future climatic projections, the above percentage (R_ir) increases due to the reduced inflows to the lake, and ranges from 0.73 to 0.57. Keeping the percentage of abstractions (and hence of irrigation pressure) unchanged, and equal to that under current conditions (26%), the available amount of water for irrigation varies from 693,989.5 to 543,581.2 m³ (Table 7). These quantities of water (and under the assumption that irrigation needs per hectare remain the same) allow for the irrigation of 100.66, 121.19, 128.52, and 118.25 ha for each of the future projections instead of 280 ha under current climatic conditions. That means that, in order for irrigation pressure to remain at the same level as it is under current climatic conditions, a reduction in irrigated land from 64.05% to 54.10% (Table 7) following the climatic projections should be considered. In fact, this reduction in irrigated land is expected to be greater, since the rise in temperature leads to higher values of evapotranspiration and hence to greater amounts for irrigation per hectare.

According to the results above, evapotranspiration increases within the basin because temperature is expected to increase, while precipitation and runoff face serious yearly reduction. The above trends can be observed in all of the four climate change models (scenario RCP 8.5) that were used. Similar findings with regard to the trends of the above parameters under different climate change models and scenarios were reported from many other researchers [67,68] for the Eastern Mediterranean region.

Furthermore, it is obvious that Lake Zazari faces a serious future problem as the lake shrinks, and available water volume is decreased. This is something that can be seen intensively in the Mediterranean area [69,70]. Lake Zazari is small compared to its basin (1:50). In addition, more than 25% of its yearly inflow is withdrawn during summer for irrigation purposes. This study showed that a reduction in rainfall of more than 25% in the basin, which is expected under severe climate change scenarios, leads to a substantial decrease in LSE. Lakes located at similar latitudes or climatic zones, especially those with similar hydrology and exploitation patterns, are expected to be subjected to analogous pressure due to climate change. As it regards Greece, there are at least 8 lakes located in the northern part of the country that have similar climatic conditions and exploitation patterns. Irrigation is one of the main aspects in the basin that would be affected, in accordance to other studies as well [71,72,73]. The results of this study can offer a quantifiable framework to the competent authorities for the design of adaptation measures related to irrigated agriculture for mitigating the impact of climate change on lake surface elevation (LSE).

4. Conclusions

This study presented an application for assessing the response of a lake to climate change using the coupled MIKE SHE/MIKE 11 modelling system in the catchment of Lake Zazari in Greece, thereby creating a fully appropriate system for wetland hydrological applications. Overall analysis and the provided results of the Lake Zazari hydrological model showed that MIKE SHE/MIKE 11 accurately simulated hydrology in the basin where both surface and ground water not only coexist, but also interact. The modelling system was also able to incorporate hydraulic structures and water management practices, and simulate irrigation at the basin scale through the combined use of surface and groundwater resources.

The calibrated model provided a detailed evaluation of the lake’s water level and its water balance components, taking into consideration the climatic conditions of the current studied period (1 January 2012–31 December 2017), but also of future climate predictions.

The four climate change GCMs showed that precipitation dropped by almost half, while temperature and thereby evaporation increased. Under these climatic conditions, runoff and hence inflow discharge into the lake decrease.

In case that irrigation demands remain the same, the lake faces diminished water level. In addition, fewer quantities of water overspill from the lake, in turn affecting the hydrology in the downstream basin.

On the other hand, any attempt to maintain the abstracted amounts of water for irrigation at the same level as it is under the current climatic conditions would lead to a reduction in irrigated area of more than half.