Future Hydrology of the Cryospheric Driven Lake Como Catchment in Italy under Climate Change Scenarios

Abstract

We present an assessment of climate change impact on the hydrology of the Lago di Como lake catchment of Italy. On one side, the lake provides water for irrigation of the Po valley during summer, and on the other side its regulation is crucial to prevent flood risk, especially in fall and winter. The dynamics of lake Como are linked to the complex cryospheric hydrology of its Alpine contributing catchment, which is in turn expected to change radically under prospective global warming. The Poli-Hydro model is used here to simulate the cryospheric processes affecting the hydrology of this high-altitude catchment. We demonstrated the model’s accuracy against historical hydrological observations, available during 2002–2018. We then used four Representative Concentration Pathways scenarios, provided by three Global Circulation Models under the AR6 of IPCC, to project potential climate change until 2100. We thereby derived daily series of rainfall and temperature, to be used as inputs for hydrological simulations. The climate projections here highlight a substantial increase in temperature at the end of the century, between +0.61° and +5.96°, which would lead to a decrease in the total ice volume in the catchment, by −50% to −77%. Moreover, there would be a decrease in the contribution of snow melt to the annual lake inflow, and an increase in ice melt under the worst-case scenarios. Overall, the annual Lake inflows would increase during autumn and winter and would decrease in summer. Our study may provide a tool to help policy makers to henceforth evaluate adaptation strategies in the area.

1. Introduction

Transient climate change largely affects the water cycle worldwide. Particularly, mountain regions are highly vulnerable, given that they are already subjected to a wide range of natural hazards, and anthropogenic, and environmental pressures [1]. Given the social and economic value of mountain areas, understanding of the hydrological regime thereby is essential. In the European Alps, temperature increased by ca. +2 °C since 1880, i.e., twice as much as the global average, and the trend has accelerated since 1970 or so, and more so at the highest altitudes [2]. The seasonal snowline therein increased by ca. +150 m for +1 °C in temperature, and less snow accumulates at the lowest altitudes [3,4,5,6,7,8,9,10,11].

Less solid precipitation in winter, and earlier snow melting in spring is causing a shift of peak runoff season, and in the last century Alpine rivers underwent significant changes in hydrology [12]. An increase in winter discharge above 1800 m asl (due to a larger share of rainfall vs. snowfall), and a decrease in spring and summer flows (due to reduction in snow cover, and increase in evapotranspiration under higher temperatures) already occurred in the last century [4,8,13].

Even though precipitation in the last few decades does not follow a clear pattern, the increase in temperature, and the subsequent increase in rainfall share against solid precipitation (snow), plays an important role in the timing of runoff, and on seasonally snow covered areas, that in turn affect the ecosystems, and glaciers’ feeding, and dynamics [4,5,6,8,14].

Indeed, an evident effect of recent global warming in the Alps is the retreat of glaciers, with a change of −30% in area in the last six decades, and an acceleration of area reduction after 1999, for glaciers larger than 0.1 km² or so. The recent increase in the number of glaciers in the Italian Alps is indeed an evidence of the fragmentation of large glaciers into smaller ones [4,5,8,9,10,11,14,15]. These glaciers, with a surface area smaller than 1 km² ca., will disappear at mid altitudes in the next few decades. However, such glaciers represent ca. 80% of the Alpine ice cover, and give a significant contribution to runoff [5,9]. Furthermore, a mass loss about 2 Gt of ice per year is expected in the Alps [14]. Glaciers with median elevation closer to their maximum (accumulation) altitude, and glaciers at lower altitudes are losing more area [11].

Such hydrological changes affect downstream water resource exploitation for socio-economic purposes, e.g., through the regulation of large lakes [1]. Our work here evaluates the effect of climate change upon a high-altitude Alpine river catchment of Italy, watering a very important lake in Northern of Italy, i.e., Lago di Como. This lake collects water, used for irrigation in the Po valley downstream, while avoiding possible flooding along the lakeshores. Accordingly, management of this lake is extremely important, under a multipurpose perspective. Hydrological variations recently are already affecting the operation of the Lake, with negative impacts on agriculture, and hydropower production. Such a situation may worsen in the future, with further increase in global and local temperature [16,17].

Some recent studies dealt with the effects of climate change on the hydrological cycle in the Alpine region [1,18,19,20,21], and some works focused upon management of water resources in the Lake Como, adopting simple models (e.g., HBV, data driven) describing the hydrology of the area [16,22,23,24]. Such models are widely adopted, in that they assess most relevant parameters usable to drive management strategies, e.g., lake inflows. However, recent studies [25,26,27] highlighted the importance of considering additional variables (snow-cover, ground water levels, precipitation, etc.), to improve management policies, also within a climate change context. Thus, the use of distributed physical models represents a useful tool, to provide a more comprehensive understanding of the lake dynamics. The Poli-Hydro, a state-of-the-art, physically based model, is well suited in this sense, since it provides a depiction of most relevant processes in cryospheric driven catchments like the one here.

Accordingly, our main goals here are to

(1) Originally set up the Poli-Hydro model, which is able to represent the main components of the hydrological cycle in the area, and their possible variation in time, providing information to policy makers;

(2) Originally analyze changes of hydrology of the upper Lake Como catchment, under most recent scenarios of AR6 of the IPCC, to provide hydrological scenarios in support to future lake management under climate change.

The paper is organized as follows. In Section 2, we describe the case study area, and methods, including data collection, and choice of the hydrological model, hydropower modelling, and hydrological projections. Section 3 contains the results of model tuning exercise, and an analysis of the climate change scenarios. Discussions and conclusions are given in Section 4 and Section 5, respectively.

2. Case study, and Methods

2.1. LakeComo Catchment

Lake Como has an estimated volume of 23.4 km³, and it is the third largest lake of Italy [16,23,25,26,28]. It is in Lombardy region, in the Pre-Alpine, and Alpine region of Northern Italy. Its contributing catchment has an area of 4438 km², nested principally within the Adda river catchment (corresponding to Valtellina valley), and the Mera river catchment (corresponding to Valchiavenna valley). The catchment is divided between Italy and Switzerland, approximately by 90–10%, respectively [22]. The Como lake collects water used for irrigation of a large part of the Po valley (ca. 70 × 10³ km² in area), one of the largest irrigated areas in Europe, with a population density of 355 habkm⁻².

The area of the lake is located within a temperate region, while the catchment has principally a cold climate. Usually, the summer is hot at low altitudes, and becomes warm, and cold at middle and high altitudes, in the Upper Valtellina valley, and in Switzerland. The climate can be classified as Cfa in the area of the lake, and as Cfb–Dfb–Dfc in the alpine part [29]. Overall, on the catchment, temperature varies between −4 °C in December and January, and +15 °C in July, with a yearly average of +5 °C. The total precipitation is around 1300 mm/y on average on the catchment, with peaks in May, and November.

Lake Como is regulated by a dam in Olginate, for a regulation volume of 254 Mm³, to supply water to irrigated agriculture downstream, over ca. 1400 km², and to at least eight run-of-river power plants downstream along the Adda river. Regulation is necessary to avoid floods on the lake shores, especially in autumn due to intense precipitation on the catchment [23,25,26,28]. Additionally, due to global warming in the last few years, conflictive operation between different stakeholders has occurred [25], in particular during summer droughts. Water scarcity requested a large release for irrigation from the lake, with respect to the normal, long term management [16,23]. This condition is in conflict with hydropower producers’ interests, that require storing more water in summer for winter power production [23].

A further main effect of climate change is the reduction in solid precipitation (snow) during winter, which in turn reduces water flows from snowmelt in spring, and leads to earlier flow peaks than normal. This implies that water has to be collected in Lake Como for a longer time period, and in larger amount, due to lack of water in summer, when agriculture requires large water volumes [28,30].

The catchment of Lake Como is Alpine, with altitude between 255 m asl, and 3822 m asl. The median altitude is 1582 m asl, 42% of the area is above 1800 m asl, and 14% over 2500 m asl.

A total of 2.7% of the area is covered by glaciers, i.e., 118 km² (Figure 1), of which 92 km² is in the Sondrio District, and the remaining in the Canton of Grisons, in Switzerland [31]. Different groups can be identified. The Ortles-Cevedale group in the East of Valtellina covers an area of 39 km² (66 km² considering the area in the Autonomous Province of Bozen), of which 10.7 km² belong to the Forni glacier [10,31,32]. The Dosdè-Piazzi group covers an area of 4 km² [11,31,33]. The Livignasco group in the North-East of Valtellina covers an area of 1.8 km² [31]. The Spluga-Val di Lei group in the North of Valchiavenna covers an area of 3.4 km² [31]. Glaciers in the Orobie group, in the center of Valtellina on the South bank, cover an area of 2.6 km² [31]. Last, the Bernina-Disgrazia group covers the remaining area between Italy and Switzerland. In this group, we can identify different important glaciers, namely Scerscen, Fellaria, Ventina, Disgrazia and the Codera-Masino sector [9,31,34].

The Adda river is largely exploited for hydropower production. Twenty-seven hydropower dams are located along the river and its tributaries [23,31], of which 16 are active as reservoirs for hydropower production [35]. More than 50% of the reservoir capacity is managed by A2A S.p.A. (a semi-private company), and Enel (the Italian National Agency for Electricity). In particular, the first one controls the main reservoirs in the Upper Valtellina valley, namely Cancano and San Giacomo, of 12 × 10⁶ and 64 × 10⁶ m³, respectively. The second one controls the main reservoir in Valmalenco (in the central Valtellina valley), e.g., the Alpe Gera reservoir of 68 × 10⁶ m³ [23,31,32,34].

The regulation of the hydropower reservoirs has an influence upon stream flows at the lake inlet, and upon water availability and timing in Lake Como. These plants contribute overall to ca. 13% of the hydropower demand of Italy, and this energy is largely used in Lombardy region [1,2,3,18].

The goal of our study here is therefore to evaluate the potential effects of climate change along the century upon the Lake Como inflows. First, we set up, tuned and validated the hydrological model Poli-hydro. Then, we used three Global Circulation Models to derive potential daily series of rainfall and temperature, to be then used as inputs for the hydrological model, to derive potential hydrological projections until the end of the century.

2.2. Hydrological Modelling

A number of presently available hydrological methods/models can be used for the assessment of high altitude, cryospheric driven hydrology [36,37,38,39]. While a review of such models and methods is beyond the scope of the manuscript here, one can briefly highlight the main issues related to the modelling of high-altitude hydrology. These include (i) representing the main processes contributing to the hydrology of the cryospheric areas, and (ii) tuning cryospheric components by way of proper data, from ground stations (snow depth, new snow precipitation, etc.), remote sensing (snow covered areas), and field assessment (ice ablation on glaciers, ice depth from GPR measurements, etc.).

With a lack of proper modelling of the main physical processes, flow prediction and projection under future scenarios may become undependable, and thereby affect policy making. Additionally, ancillary pieces of information, such as the state of the cryosphere, the dynamics of snow cover in time, and the assessment of water input as related to cryospheric dynamics, may help to refine water policy making. Thus, one needs to use a model that is well suited in this sense, in that it provides a depiction of the most important processes in the cryospheric drive catchments, and allows model tuning based upon proper data. One may refer e.g., to [33] for an in depth discussion about this topic.

Here, we used a weather-driven, semi distributed hydrological model Poli-Hydro, which was demonstrated to be dependable for the exercise of mimicking hydrology of high-altitude catchments, such as those watering the Po valley, and Lake Como here [13,33,34,40]. In Figure 2, we report a schematic view of the model, including the main processes considered, the necessary data, and possible data sources, and model validation metrics. In Figure 2, we also report the conceptual framework of our procedure. Three main blocks in the main frame are highlighted, namely 1 hydrological modelling in present conditions, 2 reservoir operation in present conditions, and 3 hydrological projections (future). In bold, we report the most important outputs for each block. Black arrows indicate the notable input passed from one block to another. Block 1, the main core of Poli-Hydro model, passes observed/modelled discharge, Q(t)/Q_m(t) to block 2, hydropower optimization, in turn delivering reservoir release R_t(t) for downstream hydrological modelling including regulation. Block 1 also delivers properly tuned hydrological modelling to block 3, for hydrological projections under future climate scenarios from GCMs.

Poli-Hydro can simulate water budget, including dynamics of glaciers, snow melt, and evapotranspiration, and subsequently provide routing of overland, and underground flow at any river section of the river network. The model mimics the variation of water content in the ground in two consecutive time steps, i.e., every day, taking as an input the liquid rain R [mmd⁻¹], the snow melt M_s [mmd⁻¹], the ice melt M_i [mmd⁻¹], and assessing evapotranspiration ET [mmd⁻¹], and groundwater discharge Q_g [mmd⁻¹]. Evapotranspiration is evaluated using Hargreaves formula, and groundwater discharge is calculated as a function of the hydraulic conductivity of soil, and of percentage soil content (with respect to a maximum value S_Max). Surface (overland) flow Q_s occurs when the actual soil storage is bigger than the greatest potential one. The Q_s of each cell is routed to catchment outlet via a Nash model [41], with two parallel systems, one for the superficial flow and one for the underground flow. Each system is characterized by a lag time t and a number n of reservoirs (normally taken as n = 3).

Snow and ice melting are modelled by a T-index approach [42,43,44,45], i.e., proportional to the atmospheric temperature (T) through a degree-day factor. Snow melt and ice melt, M_s and M_i respectively, are calculated as

$$M_s=D{D}_s\left(T_0-T_t\right)$$

(1)

$$M_i=D{D}_i\left(T_0-T_t\right)$$

(2)

where DDs is the degree day snow, DDi is the degree day ice, T₀ temperature threshold for melting, here 0 °C, and T_t [°C] temperature at time step t. Ice ablation begins after complete melting of snow cover upon ice [13,34,45,46]. To allow projections of hydrological scenarios in this area, accounting for the dynamics of ice covered area is crucial. Ice dynamic is modeled here, based upon simplified gravity-driven ice flow [33,34,47], where ice flow velocity V_ice,i is calculated as a function of the basal shear stress τ_b:

$$V_{ice,i}=K_d {\tau }_{b,i}^n h_{ice,i}+K_s\frac{{\tau }_{b,i}^n}{h_{ice,i}}$$

(3)

with h_ice,i ice thickness in cell i [m], and K_s [m⁻³yr⁻¹], K_d [m⁻¹yr⁻¹] parameters of basal sliding and internal deformation, respectively, with n exponent of Glen’s flow law [33]. Basal shear stress can be taken as

$${\tau }_{b,i}={\rho }_{i}g{h}_{ice,i}sin{\alpha }_i$$

(4)

with ρ_i ice density [kgm⁻³], g gravity acceleration [ms⁻²], and α_i local slope. Ice flow is evaluated here starting from an initial, spatially constant value of ice thickness, which from the literature can be taken here as h_ice,i = 100 m [33,34], from which the shear stress and the velocity are calculated. Then, the mass balance for each time step is performed, and the new ice thickness is evaluated as the difference between previous ice thickness, and ice melt M_i

$$h_{ice,i}^{t+\Delta t}=h_{ice,i}^t-M_i$$

(5)

2.3. Input Data

The model has a cell resolution of 500 m, and it is implemented here for the historical period 2002–2018. As inputs here, we considered a daily series of precipitation and temperature from 13 AWS stations of the Regional Agency for Protection of the Environment (ARPA Lombardia [48], Figure 1), widespread in the catchment. We spatialized over the catchment daily mean temperature, and precipitation, applying proper vertical lapse rates (at monthly scale, that were found negative for temperature as expected, and positive for precipitation within the altitude range of the catchment). The lapse rate was applied to each cell of the catchment, considering the area of influence given by Thiessen polygons built for the AWS stations.

Other input data are the digital elevation model (DEM) of the catchment (initially at 30 m cell size [49]), the land use maps from CORINE land cover [50], to estimate the maximum soil storage S_Max, and an initial map of glaciers, from the Global Land Ice Measurements from Space (GLIMS [51]).

For model tuning and validation, daily discharges and snow measurements were also needed. We took daily discharges from three available stations. The first one, on the river Adda, is located in Fuentes (SO), with data available during 2003–2018 (ARPA station). The second one is a hydrometer of Consorzio dell’Adda (a consortium managing agricultural water needs) at Samolaco, closing the Mera river (SO, for the period 2009–2018). The third station is a virtual one, i.e., we used the calculated inflows at Lake Como, obtained by the Consorzio dell’Adda, applying inverse flow-routing based on the regulated outflows, for the period 2002–2018.

We collected data of gridded snow cover, by the Moderate Resolution Imaging Spectroradiometer (MODIS [52]) on board the Terra satellite, to evaluate the snow-covered area [6,53,54]. In particular, we downloaded 42 images during 2005–2010, approx. every 15 days from 14 March to 18 July, i.e., the seasonal snow melt period, for the purpose of aiding model tuning. We then collected 28 images during 2015–2018 for model validation. For local calibration of the snow melt model, we used data of snow depth in 4 (of the 13) AWS stations of ARPA Lombardia for the period 2002–2018.

2.4. Hydropower Reservoir Modelling

Poli-Hydro simulations of high-altitude hydrology were also pursued by accounting for hydropower reservoir operation (Figure 2). Hydropower operation was modelled using an operating policy, optimized via Stochastic Dynamic Programming (SDP). The revenue from the sale of electricity was the objective function to be maximized [55,56]. In the simple case where one reservoir feeds one hydropower plant, the revenue is computed by multiplying the electricity price by the plant production, which in turn depends on the reservoir pool level, and release. In the optimization problem, the reservoir inflow is modeled as a 1-year daily cyclo-stationary stochastic process, by means of a series of lognormal probability distribution functions. The electricity price is instead modeled as a 1-year daily cyclo-stationary deterministic trajectory. The resulting optimized operating policy consists of a cyclo-stationary function, where the daily reservoir release is determined based on the reservoir storage ([35], fully depicted in [57]). This function is used to simulate hydropower release, then coupled with the hydrological model.

2.5. Hydrological Projections

To evaluate future hydrological scenarios, we performed a stochastic disaggregation [58] of three Global circulation models (GCMs), of the Assessment Report 6 (AR6) of the IPCC, e.g., EC-EARTH3 [59], CESM2 [60], ECHAM6.3 [61]. These models were driven based upon a new set of emissions and land use produced by integrated assessment models (IAMs), based on coupled shared socio-economic pathways (SSPs) and representative concentration pathways (RCPs). The SSPs represent different evolutions of future society, considering both investments in education and health, and energy development. SSPs 1 and 5 provide optimistic human development trends, but the latter at the expense of a fossil-fuel economy, whereas the former with more sustainable economy. SSP 2 envisions a future trend which continues the historical course. SSPs 3 and 4 provide a pessimistic human evolution, the first with a regional security prioritization, while the second is characterized by large inequalities. Here, we used four SSP-based scenarios as continuations of the RCP2.6, RCP4.5, and RCP8.5 forcing levels, SSP1-2.6, SSP2-4.5, SSP5-8.5, respectively, and an additional unmitigated forcing scenario (SSP3-7.0), with high aerosol emissions and land use change [62]. The model was adjusted for each GCM, based on the available daily series of precipitation and temperature for the control period 2002–2018, to obtain for each model GCM and scenario SSP a series of daily precipitation and temperature until 2100, for each AWS. The EC-EARTH, the CESM2, and the ECHAM6.3 have a spatial resolution of 0.7° × 0.7°, 1.25° × 0.9°, and 1.1° × 0.9°, respectively, thus a downscaling procedure is needed to adapt them to hydrological model resolution.

In former studies, the authors here have investigated the accuracy of GCMs models for hydrological projections. They found that previous versions of the considered GCMs, that were used for future climate simulation within the AR4/5 of the IPCC, were well suited for the purpose, and accordingly they developed proper downscaling schemes (see [21,32,40,58]). Here, the new versions of these GCMs, used under the umbrella of the AR6, were adopted.

Precipitation is downscaled with a stochastic space random cascade model [58], while temperature is downscaled using a correction with a mean monthly ΔT approach [21].

In particular, we obtained 12 scenarios usable by the Poli-Hydro to perform the simulation during 2019–2100 (Figure 2), focusing on two reference decades 2051–2060, and 2091–2100.

3. Results

3.1. Poli-Hydro Set up, and Performance

Table 1. Parameters in calibration of Poli-Hydro.

Parameter	Unit	Description	Value	Method
DDS	[mmd−1 °C−1]	Degree Day Snow	3.4	Nivometer, MODIS
DDI	[mmd−1 °C−1]	Degree Day Ice	5	Literature [32]
tg, ts	[d]	Lag times, ground/surface	150, 140	Hydrograph
ng, ns	[-]	Reservoirs, ground/surface	3, 3	Literature [32,33,34]
K	[mmd−1]	Saturated conductivity	3	Calibration
k	[-]	Ground flow exponent	1	Calibration
θw, θs	[-]	Water content, wilting, field capacity	0.15, 0.45	Literature [32,33,34]
Ks	[m−3y−1]	Ice flow basal sliding coeff.	1.5e−21	Literature [33]
Kd	[m−3y−1]	Ice flow internal deformation coeff.	1.9e−25	Literature [33]

reports Poli-Hydro model calibration. Therein are reported the model parameters, with the corresponding values, and calibration/setup method. The DDS was calibrated both in terms of snow covered area (SCA), using MODIS images, and of snow water equivalent (SWE) on the ground, using nivometric data, as reported. Figure 3 shows the estimated (%) SCA from the model, versus MODIS estimates, at the time of satellite passage (8 days bundle, max snow cover). Goodness of fit reaches R² = 0.73 during the calibration period 2005–2010 (Figure 3a), and 0.78 in the validation period 2015–2018 (Figure 3b). In terms of snow depth, the comparison between the modelled and the observed snow depth in 4 AWS stations during the entire period 2002–2018 gave R² = 0.89 (Figure 4).

The Poli-Hydro model can simulate stream flows at any (chosen) river section within the catchment. Therefore, we could simulate stream flows at the sections corresponding to the hydrometric stations of Fuentes, Samolaco, and at the lake inlet section, i.e., considering the whole lake catchments above Malgrate station, where input discharges were evaluated using inverse flow-routing by Consorzio dell’Adda. The lake inflows were used for calibration, resulting in a mean error Bias = +2.15%, and a monthly/daily Nash-Sutcliff efficiency, NSE = 0.77/0.64, respectively, during 2002–2018. In Figure 5, we report the comparison between the modelled mean monthly discharges and the observed ones at the Lake Como inlet section.

For model validation, discharge values at Fuentes, and Samolaco were benchmarked vs. observations, during 2003–2018, and 2009–2018, respectively, based on data availability (

Table 2. Goodness of fit of modelled discharges.

	Lake Inlet	Fuentes	Samolaco
Area [km2]	4438	2598	757
	Calibration 2002–2018	Validation 2003–2018	Validation 2009–2018
Bias [%]	+2.15%	−2.11%	+8.54%
NSE monthly [-]	0.77	0.69	0.80
NSE daily [-]	0.64	0.53	0.32

). We obtained a Bias = −2.11%/+8.54%, and a monthly/daily NSE = 0.69/0.53, and NSE = 0.80/0.32, respectively at Fuentes, and Samolaco.

Concerning Samolaco station, one has to report that this station closes the Mera river, flowing in the Valchiavenna valley. This area is largely regulated for hydropower production, from several reservoirs in the high altitudes. However, it was not possible to us to retrieve any information concerning the regulation of this area. Analysis of the hydrographs, and goodness of fit statistics as reported here, indicate overall agreement of the modelled discharge over the monthly/yearly scale, with worsening at the daily scale, as given by regulation. The Samolaco sub-basin covers however a somewhat low share of the full catchment (ca. 17%), so the regulation along this river, and as a consequence its Bias, should not be largely influent overall.

3.2. Simulation with Hydropower Regulation

Poli-Hydro was coupled with a hydropower regulation model in the high altitudes. We simulated daily discharges in the catchment, considering the regulation of the two main hydropower systems of the valley [32], Lake San Giacomo/Lake Cancano (Upper Valtellina), and Lake Alpe Gera/Lake Campo Moro (Valmalenco), where information was available.

Since energy price was available only during 2008–2012, we evaluated Bias, and monthly NSE between modelled, and observed flows at the inlet section of Lake Como during that period, resulting in Bias = 0.18%, and NSE = 0.75, respectively. As reported above, without regulation, we obtained Bias = +2.15%, and NSE = 0.77, during 2002–2018, apparently with no large loss of accuracy when neglecting hydropower regulation.

Since the aim of this work was to evaluate the impacts given by climate change scenarios upon the natural/undisturbed inflow at the lake, regardless of influence from hydropower regulation, we sought to develop hydrological projections without accounting for hydropower regulation in this phase, leaving assessment of the effect of regulation across the whole catchment for further refinement.

Here also, we want to assess the effect of climate change upon lake inflows over large time scales (i.e., until the end of century), so we expect that small scale (i.e., daily) changes are less relevant, and we assume that the capacity of the model to depict decently seasonal/annual flows as reported, is somewhat representative for our purpose.

3.3. Flow Components in the Lake Como Upstream Catchment

It is interesting here to explicitly assess the different flow components contributing to the lake overall inflow, which we reported in Figure 6. During the control run (CR) period 2002–2018, ca. 22% of the total discharge derived, according to Poli-Hydro simulation, from snow melting, occurring mainly in April and May, with ca. 129.3 m³ s⁻¹, and 96.3 m³ s⁻¹, respectively. Ice melt contributes less to lake inflow, and an average contribution of 0.20 m³ s⁻¹, and 0.22 m³ s⁻¹ is estimated in July, and August. The remaining contribution is made by precipitation. Evapotranspiration also plays an import role. In July, on average 104.8 m³s⁻¹ of water evaporates, and since the temperature is expected to increase, in the future this component would increase too, especially during summer months [21].

3.4. Analysis of Climate Change Scenarios

We analyzed the potential impacts of climate change, on temperature T and precipitation P in the area of interest. We projected an increase in mean annual temperature on the overall catchment, between +0.61 °C and +5.96 °C (Figure 7). For the 2.6 scenario, a larger increase (between +1.09 °C and +1.78 °C) is expected for the decade 2051–2060, with respect to 2091–2100 (between +0.61 °C and 1.12 °C), coherently with the overshoot pathway hypothesized under this scenario. According to the other scenarios, the increase in temperature is clearly larger at the end of the century, in particular for the 8.5 pathway.

Concerning total annual precipitation on the catchment (Figure 8), different projections occur for each scenario, and each GCM. Cumulated precipitation will vary between −9.53% and +11.07%, with an absolute difference of −128 mm/y to +149 mm/y. Negligible differences can be seen between the period 2051–2060, and the period 2091–2100. The EC-Earth3 model shows a decrease in precipitation for the 7.0 scenario for both periods (−3.42% and −6.23%, respectively) and an increase at the middle of the century (+0.46%), with a following decrease at the end of century (−3.69%) for the 8.5 scenario. Little differences are expected for the 2.6 and the 4.5 scenarios. The CESM2 model shows principally an increase in precipitation for all the scenarios and both the periods, up to +11.07% for the 2.6 scenario for the period 2091–2100. A decrease can be seen for the 7.0 scenario at the end of the century, i.e., −2.21%. On the contrary, the ECHAM6.3 model shows mainly a decrease in total precipitation, down to −9.53% for the 7.0 scenario in the period 2051–2060. Large decreases can be seen for the 8.5 scenario both in the period 2051–2060 (−6.20%) and in the period 2091–2100 (−7.86%).

While temperature will increase homogenously in every month, precipitation shows different patterns with months, and scenarios (not reported for shortness). Indicatively, precipitation will increase in February and November in the middle, and at the end of the century for the majority of scenarios and models. Some models for some scenarios show a decrease in March, April, August, September, and October in the period 2051–2060 and, in addition, May, June and July in the period 2091–2100.

In Figure 9 we display the effects of modified climate patterns upon snowpack dynamics. We report therein changes of the average (catchment wide) volume of SWE_cum (snow water equivalent in the snowpack) in 4 months (April, June, September, December).

These represent the main seasons for snowpack in the area, namely snowpack peak, snow thaw season, end of snow melt, and start of snow accumulation. The percentage variation is given as an average on the scenarios of all GCMs.

Looking at the four graphs in Figure 9, we clearly see different behaviors, depending upon the scenarios and the considered decade. For the 2.6 scenario we expect principally an increase in volume of the average snowpack on the catchment, due to an increase in solid precipitation, thanks to the limited decrease or to the increase in temperature. On the contrary, for the 8.5 scenario we expect a strong decrease in solid precipitation, and consequently of the cumulated snow SWE_cum, during the entire year, and at every altitude range. The same result is obtained for the 4.5 and the 7.0 scenarios at the end of the century. These two scenarios in the period 2051–2060 have a different behavior. We can see a decrease in SWE_cum in April (Figure 9a), due to the increase in temperature at low altitudes, that will cause less accumulation of snow and a faster melting at the beginning of the spring. At high altitudes, the quantity of snow remains more or less constant with respect to the reference period, so the overall decrease is due to a decrease at low altitude. This was verified by looking at the SCA in the same month, that is smaller with respect to the CR period and concentrated to high altitudes (not shown for shortness). In June (Figure 9b) there is not a significant variation of SWE_cum with respect to the CR period. This is explained by a small variation of T in the months of April, May, and June (not reported) and to the maintenance of the snow at high altitudes, with the same behavior as seen in the control run period. The increase in SWE_cum in September (Figure 9c) is a little absolute variation, indeed the SWE_cum in September is the lower value of the year. The little increase in T in summer (July, August, September, not reported for shortness) is limited to low altitudes, so at high altitudes snow continues to be present. Again, we compared the SCA in September for the CR period and for each scenario, and snow is concentrated at higher altitudes. In the end, in December (Figure 9d) SWE_cum presents a decrease, due to reduction in solid precipitation, and snow accumulation at low altitudes.

In Figure 10 we report the projected trend of total glaciers volume, since 2020 to the end of the XXI century. For each decade, the average volume of glaciers on the catchment is calculated for each model and each scenario. The chart reports the average (among GCMs) on the four scenarios. The 2.6, 4.5 and 7.0 scenarios have a similar behavior. An initial strong decrease in ice cover is seen, with a plateau after 2050. The 8.5 scenario would project a slower decrease in the first part of the century, and a more significant decrease after 2050. At the middle of the century all scenarios show approximately the same ice volume, at the end of the century the 8.5. scenario volume is the lowest, due to the large increase in radiative forcing therein. In

Table 3. Percentage variation of ice volume with respect to the CR period.

Variation of Total Glaciers Volume [Δ%]
	2051–2060	2091–2100
earth26	−56.79	−57.27
earth45	−55.43	−60.64
earth70	−57.92	−72.45
earth85	−58.04	−77.62
cesm26	−49.28	−50.35
cesm45	−49.19	−51.68
cesm70	−49.76	−52.93
cesm85	−51.67	−57.12
echam26	−52.43	−51.31
echam45	−53.08	−52.94
echam70	−53.49	−54.35
echam85	−52.04	−56.34

the percentage variation of glaciers volume is reported for each scenario, and each GCM in the 2051–2060 period, and in the 2091–2100 period. The EC-Earth3 model projects the worst conditions for all the scenarios, with a decrease down to −77.62% for the 8.5 scenario at the end of the century.

In Figure 11, we report future projections of mean monthly discharge for each GCM and each scenario, at half century, and at the end of the century, with respect to present condition (CR period). The increase in liquid precipitation, at the expenses of solid precipitation in the winter period (December, January, February, March) would cause an increase in stream flows therein. The increase in discharge in November is due to increase in precipitation. On the contrary, the decrease in discharge between May and September is due to decreased precipitation, especially in the last decade of the century, and to lack of snowmelt. The average annual discharge, Q_y = 150 m³s⁻¹ in the CR period, will mainly decrease in the future, as reported in Figure 12. Therein, we report Q_y for all GCMs for each scenario, at half century (solid colour) and end of century (stripes). In

Table 4. Projected share of flow components for each GCM and each SSP.

	2051–2060				2091–2100
	Rain [%]	Snow [%]	Ice [%]	Q Mean [m3/s]	Rain [%]	Snow [%]	Ice [%]	Q Mean [m3/s]
earth26	82.74	16.80	0.45	140.43	82.42	17.51	0.07	141.08
earth45	85.66	13.52	0.82	145.55	84.58	15.16	0.26	139.09
earth70	85.56	13.63	0.81	143.24	89.68	9.07	1.25	143.24
earth85	86.35	13.03	0.62	160.10	89.92	8.61	1.47	141.21
cesm26	81.65	17.63	0.72	153.35	82.39	17.55	0.07	166.88
cesm45	84.28	14.78	0.94	140.25	84.59	15.32	0.09	143.91
cesm70	84.57	14.64	0.79	159.59	84.64	15.17	0.19	137.68
cesm85	85.79	13.44	0.76	155.66	90.26	9.00	0.74	142.15
echam26	84.82	14.68	0.49	139.22	82.13	17.83	0.05	148.90
echam45	84.85	14.63	0.52	140.84	85.72	14.18	0.10	133.90
echam70	86.43	13.13	0.44	158.58	88.45	11.37	0.18	144.68
echam85	84.76	14.36	0.88	137.79	90.08	9.33	0.59	129.16

we report Q_y for each scenario and each GCM. The overall trend, especially at the end of the century, is characterized by a decrease in the mean annual discharge. However, during 2051–2060, some scenarios show a slight increase in discharge, mainly due to the increase in precipitation they provide at the middle of the century (Figure 8).

In Figure 12 we can also see the contribution of (liquid) precipitation, snow melting and ice melting to the annual discharge, on average between models. During the CR period the calculated contributions were equal to +53.1%, +34.2% and +2.7%, respectively. As we can see in the Figure 12, and in details in Table 4, for each scenario and each GCM, the contribution of ice melting will decrease significantly. The contribution of ice melt would be somewhat higher for the model EC-Earth3 under the 7.0 and 8.5 scenarios, where glaciers’ volume reduction is largely enhanced, however clearly resulting into ice cover depletion, and lower ice melt later on. At the end of the century, the 8.5 SSP scenario displays higher ice melting contribution with respect to other scenarios. Indeed, the 8.5 SSP projects a further decrease in glaciers’ volume after 2050, differently from other scenarios (Figure 10). This reduction in ice volume would result into a larger contribution of ice melting to stream flows. The contribution of snow melting will decrease as well, as expected. It will contribute between +8% and +18%, according to our scenarios and GCMs. On the contrary, the contribution of precipitation, in particular rain, will increase, up to +80%, +90%.

4. Discussion

Looking at our results, our Poli-Hydro model simulates acceptably well Lake Como inflows, despite some background noise, also given by high altitude hydropower regulation. Accordingly, we thought to use the model to preliminarily analyze the effects of climate change upon the hydrology of the Lake Como catchment. A credible projection of future water availability in the lake is essential for planning of its future management, even considering the different future demands of downstream stakeholders, in particular for irrigation [16].

Proper spatialization of the meteorological data was necessary to consider the effects of orography. Temperature clearly displays negative lapse rates with respect to altitude, but it is horizontally quite uniform in the catchment. On the contrary, we observed a certain variability of precipitation within the different valleys of the catchment. The collection of meteorological data from a large number of AWS, and acceptably widespread in the catchment, allowed us to consider differences in the precipitation distribution, paramount important for proper hydrological modeling in high altitudes [14].

For Poli-Hydro model calibration we used some parameters (namely, for ice flow dynamics, and degree day factor for ice) coming from other studies in the same area [25,26,27], e.g., the upper Valtellina valley, and the area of Dosdè, and the Valmalenco valley. Other parameters (such as the degreed day factor for snow) were calibrated for the specific case study.

Some approximations were made, depending upon data availability. We have hypothesized an average initial ice depth, homogeneous within for the entire catchment, based upon a detailed analysis on the glaciers in the study area. This was necessary with the lack of specific distributed data on the overall area, but consistent with present studies covering glaciers’ dynamics therein [33,34]. Moreover, considering that the ice melt contributes to 2.7% to the total annual discharge in the CR, a potential underestimation or overestimation of the ice depth, may not provide large noise.

Snow and ice melt models are based only upon accumulation of thermal time. The thermal melting factor (DDS and DDI) give an acceptable approximation at the daily scale, while physically based model mimicking energy budget would likely provide a more accurate description of cryospheric melt even at shorter scales (e.g., hourly). Degree-day models, however, clearly require a limited amount of data with respect to the full energy balance model [6], and here we used a daily scale of investigation. In this study the temperature-radiation approach is simplified using only the average daily temperature for each time step t. This hypothesis is acceptable considering the orography of the catchment. Using a spatial grid of 500 m × 500 m, the local differences in orography cannot be considered, just as the effects on the ground radiation. Moreover, in the study area few data of radiation are available from AWS.

The degree-day factor DDS is here calibrated and validated, using both observed ground data and remote sensing data (MODIS). Modeling of current accumulation and melting is difficult, due to the high spatial variability of the territory, in terms of altitude, and topography. For this reason, satellite data are useful to control the seasonality (timing) and the distribution in space of snow cover area (SCA). It is broadly known that calibration of snow models is improved including satellite data of snow cover [46,53], and accordingly here we could exploit such added value.

Here, we focused upon assessing potential changes in the hydrology of the upper catchment of Lake Como, with reference to natural (i.e., unregulated) stream flows. We demonstrated that the model provides an acceptable representation of the Lake Como inflows (as seen by low yearly Bias_% at inlet, and somewhat high NSE at monthly scale therein, Table 2, and visual assessment of monthly flows in Figure 5b), even without considering the upstream hydropower regulation. No large loss of accuracy occurs against the regulation scenarios (as reported in Section 3.2).

The lack of regulation can be partially seen in the daily NSE, but here one may state that the main interest is linked to seasonal and monthly variation of discharge and its contribution. Looking at Figure 5, we can notice the effects of regulation in the abundance of observed discharge in winter (JFM) and November, and the lower discharge in summer (JAS). Regulation upstream is complex to model, and furthermore it may change in the future, so making projections of flow dynamics is possibly less accurate. However, we demonstrated here that overall regulation does not largely affect the performance of the model, in the simplifying hypothesis that “natural like” discharges reach the lake. Accordingly, we assumed that under the simplifying hypothesis of “natural like” regime, assessment of climate impact upon the flow regime upstream of the lake can be pursued using the Poli-Hydro model, at least preliminarily, and even considering large uncertainty therein.

Looking at the future projections, we can benchmark our results against similar studies in the Valtellina valley [32,34]. As from literature, the expected contribution of snow melt and ice melt in the future would drastically reduce, for the reduction in solid precipitation in winter and the increase in glaciers ablation [12,32,34]. The survival of glaciers at the end of the century is limited to the high altitude part, usually above 3000 m asl [32], and to the largest glaciers of the catchment [34].

Annual discharge Q_y changes with GCMs and RCPs/SSPs, but within the somewhat narrow range of −14% to +10%. In all cases, both at mid-century and at end of century, it will increase in autumn and winter (from November to March), and it will decrease in summer (from May to September). Again this is coherent with other studies developed under the Fifth Assessment Report of the IPCC [12,32,34]. This can be explained with different reasons, namely (i) increase in liquid precipitation vs. solid precipitation in winter, especially at low altitudes, (ii) earlier, and shorter snow thaw in spring due to the increase in temperature, and the lack of snow cover, (iii) changes in precipitation timing seasonally.

Different dynamics of seasonal (monthly) discharge during the year will likely cause a different management of the Lake Como. A large amount of water has to be collected during autumn and winter, to have water available for summer irrigation demand (also likely to increase due to larger evapotranspiration in the Po valley). Such need may conflict with flood risks prevention, i.e., with the need for lowering the lake level to avoid shoreline floods. It is likely that more detailed studies have to be pursued, to verify the chance for correct management of the Lake Como, considering future hydrology, according to our projections.

In the approach we proposed here, the use of Poli-Hydro model may help largely in this direction. Poli-Hydro model was used in several areas of the European Alps, and outside, including in several areas of the Himalayas, and South America, always with good results [46,47]. Poli-Hydro, and similarly conceived, physically based hydrological models, whenever well fed with proper inputs, may provide large flexibility of application, and thereby give large insight to scientists, and policy makers in the area.

Possible draw backs for such models may dwell into the large computational time as required for physically based modelling, and projections. However, Poli-Hydro here is rather fast, and can be used for long term, multi-scenarios simulation as here, with reasonable computational time (few hours of simulation for each IPCC scenario), and hardware investment (all simulations were made using a desktop PC), so allowing the exploitation of a large potential array of scenarios for policy making.

5. Conclusions

Lake Como of Italy is an iconic, worldwide renowned place, for historical, cultural, and touristic reasons. Lake operation has tremendous importance therefore, contributing to the beauty of the area, the feasibility of recreational and sport activity, and most importantly supplying water downstream for irrigation, and hydropower, and mitigating flood risk along the lake shores (and possibly downstream in the Adda river). A great deal of literature exists exploring the optimal operation of the lake under present, and potential future hydrological scenarios. However, the complex, cryospheric driven hydrology of the lake upstream catchment calls for proper modeling of the processes therein. Here, we used the Poli-Hydro model, a state-of-the-art model, fully capable of mimicking the most relevant processes of high altitude, cryospheric hydrology of the catchment, including complex snow and ice melt, and glaciers’ flow dynamics, based upon a large data base from field measurements, data at stations, and remote sensing information.

We further used most recent climate projections from models, and scenarios of the IPCC AR6, the state-of-the-art reference in the field, to provide credible scenarios of climate for hydrological projections in the area. Our results, even within the range of the well known uncertainty when dealing with future climatic, and hydrologic scenarios [32], indicate consistently the expectation of an increase in flows during the wet (flood) seasons, winter and especially fall, and subsequent decrease during the dry (drought) seasons, spring, and especially summer, as due to shifted snow cycle, and decreased ice cover. Such scenarios would largely impact the future dynamics of the lake, and surely would call for modified operation thereby. Our results here may therefore provide a credible tool for initial brainstorming of future lake operation under impending climate change, for scientists, and policy makers therein. Moreover, the proposed methodology is useful to evaluate climate change effects on hydrological components and it can be applied to other case studies, such as for large lakes with multipurpose management in northern Italy, and in worldwide mountain areas.