Hydrogeological Study of the Glacial — Fluvioglacial Territory of Grandate ( Como , Italy ) and Stochastical Modeling of Groundwater Rising

On November 2014, the Municipality of Grandate, near Lake Como, had to deal with a great emergency that was caused by the flooding of factory undergrounds. The authors realized a hydrogeological study to understand the causes of groundwater flooding and to prepare a pre-feasibility study concerning possible actions for groundwater control. The hydrogeological structure is rather complex and required time-consuming reconstruction of the conceptual site model. A transient numerical model was developed to analyse the system behaviour in different scenarios. The flow model was calibrated in a steady and unsteady-state using the automatic calibration code Model-Independent Parameter Estimation (PEST). The study demonstrated that the reason for floods was mainly due to the concurrence of three causes: (1) the hydrogeological structure of the area was recognized as a stagnation zone, (2) groundwater rising, and (3) extremely heavy rainfall in 2014. Through the PEST RandPar function, 100 random rainfall scenarios were generated starting from rainfall data for the last 20 years. The model was used to run 100 1-year long simulations considering the probability distribution of recharge related to the 100 randomly generated rainfall scenarios. Through collecting the piezometric heads that resulted from the simulations, monthly probability curves of groundwater exceeding a threshold level were obtained. The results provided an occurrence probability of groundwater level exceeding the underground structures level between 12% and 15%.


Introduction
The development of industry in several major cities in Europe caused groundwater levels to drawdown until they were tens of meters below ground level.As the water tables were so deep, it was often assumed that they were stable and could be safely ignored.Over the last 10-20 years, a reduction in groundwater abstraction has led to a rise in water levels almost everywhere in Europe [1][2][3][4].Then, since the end of the last century, many European cities have been affected by the groundwater rising phenomenon that is the result in many causes, such as the deindustrialization process, urbanization, and the variation of the aquifers recharge rate because of climate change.These different conditions have led to local rises of the water table, ranging from 3 m to 30-40 m with different consequences [5], with the flooding of underground structures among them.The rising groundwater trends are registered especially in particularly exposed areas, such as the depressions that in the recent past have been the site of lakes, swamps, and peat bogs.On November 2014, the Municipality of Grandate, located 10 km southward of Lake Como, had to deal with a great emergency that was caused by the flooding of numerous undergrounds of buildings and factories that were located in that area.Flooding produced several problems in terms of direct and indirect damages to the economy and the territory.
Grandate is located in the Prealps area and in a rather complex geological context.The hydrogeological asset of glacial amphitheatres in the Lombard Prealps has some properties that lead to genesis of stagnation areas [6][7][8], caused by the convergence of groundwater flow from the glacial hills (moraines) towards the morphological depressions lying between the morainic arches, that are filled by fluvioglacial and colluvial sediments.Groundwater stagnation areas represent a zero flux and equilibrium points where applied forces are balanced in all directions [9,10].The convergence of the flow lines towards these areas makes their drainage difficult, lowering the seepage velocity and generating groundwater stagnation.For these reasons, along the morphological depressions, groundwater level often approaches the ground surface, even generating springs, lakes (Sartirana and Alserio lakes are an example in the Prealps area), and swamps.Analysis and studies about stagnation points and their possible practical applications have been conducted by many researchers [7,[11][12][13][14][15], providing that, in order to avoid the risks for existing and projected buildings related to this phenomenon, a good knowledge of the system and its behavior is the first requirement.
Through the hydrogeological study that was carried out during this work (November 2015-November 2016), the Grandate area has been characterized for the first time-no continuous monitoring data existed in this zone before 2015-and has been recognized as an example of stagnation area.The reconstructed hydrogeological conceptual model has been the base for the implementation of a 3D numerical one that has been run for complete knowledge of the system behavior.Thereafter, the model has been used in order to assess the hydrogeological hazard for the underground structures because of the water table rise.Different scenarios of groundwater system evolution have been defined based on stochastically generated rainfall scenarios.The recharge has been considered as probability distributions and Monte Carlo simulations [16] have been carried out to obtain the probability distributions of the groundwater levels.Therefore, the occurrence probability of groundwater exceeding a defined threshold level has been assessed.Uncertainty is inherent in any model-based risk assessment [17] and this uncertainty is particularly high when analyzing future impact [18,19].Regardless, the current paper demonstrates that a hydrogeological analysis coupled with mathematical modelling and stochastic analysis can provide satisfying approximate forecasting of the behavior of the piezometric levels and can be a useful tool for the assessment of risk related to groundwater flooding.

Geology of the Study Area
The Como Province is located in the North-Central area of the Lombardy Region (Figure 1).This part of the region is a recharge zone for the aquifers of the exploited area of Milan and is recharged by the zones in contact with the pre-Alpine rising area.
The hilly and high plain of the Como province, given the geological importance that the glacial deposits-here well preserved-acquire in this territory and the peculiarity of their distribution, has been at the center of the scientific interest of many researchers, Italians, and foreigners since the beginning of the nineteenth century.In particular, Penck and Brückner [20] recognize four glaciations in this area called, from the most recent, Würm, Riss, Mindel, and Günz which are separated by three interglacial ones called: Riss-Würm, Mindel-Riss, and Günz-Mindel.Ugolini and Orombelli [21] and Billard [22] at first emphasize the correlation between morainic and fluvioglacial deposits, examining in detail the characteristics of the main soils.In 1987, Bini [23] reviewed the model of the four glaciations of Penck and Brückner and, through a detailed survey, divided the sedimentary bodies into "glacial complexes".The studies that have been undertaken [23][24][25] show that, since the late Pliocene, the Lake Como area was subjected to 13 glacial episodes during which the glaciers extended to the upper Po Plain.However, given the application purposes of this work, the traditional classification is used in this study.Because of its fluvio-glacial origin, the territory of the Como province discloses singular morphological characteristics and a rather complex structure.In the late Miocene, the pre-Alpine arc was subjected to a strong erosive stage with the consequent formation of canyons at existing prealpine lakes.During the Pliocene highstand sea level stage, clay (Villafranchiano) filled the low areas in the pre-Alps.From late Pliocene-lower Pleistocene there was a low stand sea level with glaciers that moved southwards and occupied the pre-alps zone, and the area was covered by fluvioglacial coarse sediments that subsequently lithified (Ceppo Lombardo).During the Quaternary, some glacial-eustatic phases occurred and the bedrock at the foot of the Como pre-Alps was deeply etched by glacial masses coming from Canton Ticino (glacier-Faloppio, Italy) and Lario (Como side).The enormous amount of sediments that were transported by glaciers and meltwater filled the lowest areas and formed the first morainic amphitheaters.Later, due to the decrease in the glacial phenomena, internal and less extensive moraines were formed.Generally, the most recent sediments of the morainic amphitheaters are located at a lower altitude than the oldest ones.The hydrographic network here is mainly developed in a North-South direction.
From the hydrogeological point of view, the area has marginally been taken into consideration by the Consorzio Acque Potabili (CAP) publication [26] and by some subsequent studies [27] where authors examine the southern limits of the province bordering the Milan one.A summary of the hydrogeological aspects of the Como area is presented by Francani et al. [28,29] and Beretta [30].These studies recognize three main aquifers in the area and the presence of major water resources within paleo channels.
A similar geological and hydrogeological structure is that of the main Lombard morainic arches: Lambro river, Molgora stream in Merate and Sartirana, Oglio river, and river Mincio.In all of these cases, in the valley between a morainic arch and the other, ponds, swamps, and even lakes can be noticed [31,32].The presence of piezometric depressions in correspondence with lacustrine basins is frequently reported in detailed geotechnical studies, like that by CEPAV regarding the glacial deposits of Garda Lake [33] and the study by Arpa Lombardia [34].
The Grandate municipality is located less than 10 km south of Lake Como and its territory presents the characteristics that have already been described.From Figure 2a, in fact, it is possible to recognize two morainic arcs that enclose the less topographically elevated zone of the Grandate plain.Seveso River crosses the plain from North to South.In the area, groundwater flows generally from NNW to SSE.For this study, a bigger area has been considered in order to implement the mathematical model (red square in Figure 2).Because of its fluvio-glacial origin, the territory of the Como province discloses singular morphological characteristics and a rather complex structure.In the late Miocene, the pre-Alpine arc was subjected to a strong erosive stage with the consequent formation of canyons at existing pre-alpine lakes.During the Pliocene highstand sea level stage, clay (Villafranchiano) filled the low areas in the pre-Alps.From late Pliocene-lower Pleistocene there was a low stand sea level with glaciers that moved southwards and occupied the pre-alps zone, and the area was covered by fluvioglacial coarse sediments that subsequently lithified (Ceppo Lombardo).During the Quaternary, some glacial-eustatic phases occurred and the bedrock at the foot of the Como pre-Alps was deeply etched by glacial masses coming from Canton Ticino (glacier-Faloppio, Italy) and Lario (Como side).The enormous amount of sediments that were transported by glaciers and meltwater filled the lowest areas and formed the first morainic amphitheaters.Later, due to the decrease in the glacial phenomena, internal and less extensive moraines were formed.Generally, the most recent sediments of the morainic amphitheaters are located at a lower altitude than the oldest ones.The hydrographic network here is mainly developed in a North-South direction.
From the hydrogeological point of view, the area has marginally been taken into consideration by the Consorzio Acque Potabili (CAP) publication [26] and by some subsequent studies [27] where authors examine the southern limits of the province bordering the Milan one.A summary of the hydrogeological aspects of the Como area is presented by Francani et al. [28,29] and Beretta [30].These studies recognize three main aquifers in the area and the presence of major water resources within paleo channels.
A similar geological and hydrogeological structure is that of the main Lombard morainic arches: Lambro river, Molgora stream in Merate and Sartirana, Oglio river, and river Mincio.In all of these cases, in the valley between a morainic arch and the other, ponds, swamps, and even lakes can be noticed [31,32].The presence of piezometric depressions in correspondence with lacustrine basins is frequently reported in detailed geotechnical studies, like that by CEPAV regarding the glacial deposits of Garda Lake [33] and the study by Arpa Lombardia [34].
The Grandate municipality is located less than 10 km south of Lake Como and its territory presents the characteristics that have already been described.From Figure 2a, in fact, it is possible to recognize two morainic arcs that enclose the less topographically elevated zone of the Grandate plain.Seveso River crosses the plain from North to South.In the area, groundwater flows generally from NNW to SSE.For this study, a bigger area has been considered in order to implement the mathematical model (red square in Figure 2).

Conceptual Site Model Reconstruction
The reconstruction of the conceptual site model started from the geometries of aquifers that were provided by previous works and were improved through the analysis of 110 stratigraphic logs, which have been used to obtain an accurate representation of the hydrogeological structure of the study area (Figure 2b).The collected stratigraphic information has been compared with data from previous studies [30,35], with the topographic data of the Digital Terrain Model (DTM) [36,37] and with the geological map that was elaborated for the CARGproject (Cartografia Geologica e geotematica) by the Istituto Superiore per la Protezione e la Ricerca Ambientale (ISPRA) [32].For this study, eleven cross-sections have been prepared, six with a NS direction (1-6) and five are WE oriented (A-E), as shown in Figure 2b.Section 4 and Section C are shown in Figure 3.

Conceptual Site Model Reconstruction
The reconstruction of the conceptual site model started from the geometries of aquifers that were provided by previous works and were improved through the analysis of 110 stratigraphic logs, which have been used to obtain an accurate representation of the hydrogeological structure of the study area (Figure 2b).The collected stratigraphic information has been compared with data from previous studies [30,35], with the topographic data of the Digital Terrain Model (DTM) [36,37] and with the geological map that was elaborated for the CARGproject (Cartografia Geologica e geotematica) by the Istituto Superiore per la Protezione e la Ricerca Ambientale (ISPRA) [32].For this study, eleven cross-sections have been prepared, six with a NS direction (1-6) and five are WE oriented (A-E), as shown in Figure 2b.Section 4 and Section C are shown in Figure 3.
What is possible to understand from the hydrogeological sections is that the Grandate area is a valley that is carved in the rocky substratum of Würm glaciers into which the glacial and alluvial deposits have accumulated.This formation is generally considered semipermeable because there are very few wells with screens in this unit and they are only in the most superficial and fractured part of the rock.Moving southward, the bedrock deepens, covered by two aquifers: the first, generally phreatic, corresponds to the gravelly-sandy complex, while the second, confined, is composed by the fluvioglacial Riss/Mindel deposits and the Villafranchian ones and is found more than 100 m below ground level.The first aquifer is separated from the deepest one by the glacial Riss/Mindel deposits that are about 30 m below the ground level.In the southern part of Grandate, the rocky substratum and the Villafranchian deposits rise again.The valley in the substratum tends in fact to close in this zone where the separation level between the first and the second aquifer also disappears.The aquifers here are thus in contact and reduce their thickness, making the outflow section southward thinner.From Section 4 and Section C (Figure 3), it is possible to notice what was previously explained: the Grandate plain is like a big bucket, surrounded by morainic hills, where water flows from the topographically more elevated zones, infiltrates in the plain, and hardly discharges southward.In fact, for the uplift of the glacial Riss/Mindel deposits and the approaching of the morainic arcs, the outflow section shrinks in the southern zone of the plain, making it difficult for groundwater to flow downgradient.area (Figure 2b).The collected stratigraphic information has been compared with data from previous studies [30,35], with the topographic data of the Digital Terrain Model (DTM) [36,37] and with the geological map that was elaborated for the CARGproject (Cartografia Geologica e geotematica) by the Istituto Superiore per la Protezione e la Ricerca Ambientale (ISPRA) [32].For this study, eleven cross-sections have been prepared, six with a NS direction (1-6) and five are WE oriented (A-E), as shown in Figure 2b.Section 4 and Section C are shown in Figure 3.In November 2015, a piezometric survey over 31 wells/piezometers and 5 hydrometers along Seveso River was carried out.The data that has been collected are not numerous because unfortunately in the past, the area has never been characterized from the hydrogeological point of view and there are not many monitoring points.The piezometric map (Figure 4) confirms the conceptual site model reconstruction: the first aquifer, though presenting a moderate hydraulic conductivity, hosts a groundwater body that is sustained by finer sediments, which confine the second aquifer.Part of this groundwater flows toward the lake (as happens in some similar systems [38]) and part thereof flows southward through the Seveso River depression aquifer.Most of the groundwater, however, flows from the borders of the morainic hills that surround the plain toward the center of the plain and, with difficulty, discharges southward where the flow section shrinks.The low gradient (<0.5%) that is registered in the plain (Figure 4) is a consequence of what was previously explained, because this hydrogeological structure has been recognized in many cases as a stagnation area.In Figure 5, the groundwater levels that were monitored from March 2016 to present in two piezometers that are located in the flooded zone are represented.One of the two piezometers is monitored continuously through a data logger, while in the other, the groundwater level is manually recorded approximately every 15 days.The graph clearly shows a very fast response of the groundwater to no single rain events-especially when these exceed 20 mm-which is probably due to a very small effective porosity of the sediments hosting the groundwater body, but is undoubtedly caused by the fact that the zone is a stagnation area where the outflow is reduced due to the hydrogeological setting.Single events like the heavy rains in November and December of 2017 do not create great piezometric level variation as in the other cases, however an interruption in the downward trend can be observed.In this study, the piezometric map, which was constructed by interpolating the November 2015 survey data (Figure 4), shows that a stagnation area would probably be present in the Grandate plain.
For this reason, the case of Grandate-which has a hydrogeological setting that is similar to that of an extensive part of the Lombard territory, as indicated in the already mentioned studies [31][32][33][34]has been examined in detail through the implementation of a numerical model.The implementation phase was carried out taking into account the hypothesis of the stagnation area existence which was supported by the piezometric map and the geological reconstruction of the area.The mathematical model then allowed a better understanding of the factors that affect groundwater flow and helped to evaluate the most suitable interventions for avoiding future groundwater flooding.In this study, the piezometric map, which was constructed by interpolating the November 2015 survey data (Figure 4), shows that a stagnation area would probably be present in the Grandate plain.
For this reason, the case of Grandate-which has a hydrogeological setting that is similar to that of an extensive part of the Lombard territory, as indicated in the already mentioned studies [31][32][33][34]has been examined in detail through the implementation of a numerical model.The implementation phase was carried out taking into account the hypothesis of the stagnation area existence which was supported by the piezometric map and the geological reconstruction of the area.The mathematical model then allowed a better understanding of the factors that affect groundwater flow and helped to evaluate the most suitable interventions for avoiding future groundwater flooding.Generally, after finding out and outlining the stagnation zones based on the flow network reconstruction, it is essential to understand what the physical phenomena are that control the origin and presence of zero-velocity areas.This knowledge is useful for correctly interpreting the model results that are related to the behavior of the aquifer and the effects of potential drainage systems on groundwater levels (Figure 6).In fact, in order to achieve optimal results in different pluviometrical and hydrological conditions, the design of drainage interventions in these kinds of areas requires more than the usual preliminary hydrogeological analysis of the most suitable water well location and a correct evaluation of the withdrawal and the water level lowering.

Mathematical Model
The conceptual site model has been implemented in Modflow2005 [39] throughout an 80 × 80 m grid and was refined up to 40 × 40 m in the area where flooding occurred in 2014.In the vertical direction, the model has been extended based on the stratigraphic characteristics and the geometry of the aquifers in the area.Four layers have been identified according to the defined hydrogeological sections (Figure 7): the first represents the shallow aquifer that is constituted by the permeable deposits of the glacial and fluvioglacial Würm phase, the second follows the impermeable layer, where present, the third and fourth layers instead represent the deepest part of the hydrogeological system (fluvioglacial Riss/Mindel and Villafranchian, respectively).The bottom of the model is represented by where the rocky formation where present and was set at a distance of 70 m from the layer above where the identification of the substratum was not possible.The grid thus is composed by 87 rows, 76 columns, and 4 layers (26,448 active cells), representing a 40 km 2 area.General Head Boundary conditions (GHB) have been set over the borders of the model domain based on the data that was collected during the piezometric survey (November 2015) which covered a larger area than the model.GHB have been chosen because of the absence of geological or hydrogeological natural limits in the proximity of the modelled area and in order to make the influence of the boundary conditions on the target zone as negligible as possible.The River condition (RIV package) was used to represent the Seveso River (Figure 7) that is a perennial river originating about 10 km North of Grandate and presents, in this zone, a small river stage.The Seveso River has little influence on groundwater, as is possible to notice from Figure 4.The River condition was set on the base of the measured levels (5 points along the river path) and, for the river bed, considered a hydraulic conductivity value that was similar to that of the surrounding sediments [40].In this study, the piezometric map, which was constructed by interpolating the November 2015 survey data (Figure 4), shows that a stagnation area would probably be present in the Grandate plain.
For this reason, the case of Grandate-which has a hydrogeological setting that is similar to that of an extensive part of the Lombard territory, as indicated in the already mentioned studies [31][32][33][34]-has been examined in detail through the implementation of a numerical model.The implementation phase was carried out taking into account the hypothesis of the stagnation area existence which was supported by the piezometric map and the geological reconstruction of the area.The mathematical model then allowed a better understanding of the factors that affect groundwater flow and helped to evaluate the most suitable interventions for avoiding future groundwater flooding.

Mathematical Model
The conceptual site model has been implemented in Modflow2005 [39] throughout an 80 × 80 m grid and was refined up to 40 × 40 m in the area where flooding occurred in 2014.In the vertical direction, the model has been extended based on the stratigraphic characteristics and the geometry of the aquifers in the area.Four layers have been identified according to the defined hydrogeological sections (Figure 7): the first represents the shallow aquifer that is constituted by the permeable deposits of the glacial and fluvioglacial Würm phase, the second follows the impermeable layer, where present, the third and fourth layers instead represent the deepest part of the hydrogeological system (fluvioglacial Riss/Mindel and Villafranchian, respectively).The bottom of the model is represented by where the rocky formation where present and was set at a distance of 70 m from the layer above where the identification of the substratum was not possible.The grid thus is composed by 87 rows, 76 columns, and 4 layers (26,448 active cells), representing a 40 km 2 area.General Head Boundary conditions (GHB) have been set over the borders of the model domain based on the data that was collected during the piezometric survey (November 2015) which covered a larger area than the model.GHB have been chosen because of the absence of geological or hydrogeological natural limits in the proximity of the modelled area and in order to make the influence of the boundary conditions on the target zone as negligible as possible.The River condition (RIV package) was used to represent the Seveso River (Figure 7) that is a perennial river originating about 10 km North of Grandate and presents, in this zone, a small river stage.The Seveso River has little influence on groundwater, as is possible to notice from Figure 4.The River condition was set on the base of the measured levels (5 points along the river path) and, for the river bed, considered a hydraulic conductivity value that was similar to that of the surrounding sediments [40].
conditions on the target zone as negligible as possible.The River condition (RIV package) was used to represent the Seveso River (Figure 7) that is a perennial river originating about 10 km North of Grandate and presents, in this zone, a small river stage.The Seveso River has little influence on groundwater, as is possible to notice from Figure 4.The River condition was set on the base of the measured levels (5 points along the river path) and, for the river bed, considered a hydraulic conductivity value that was similar to that of the surrounding sediments [40].The morainic arc that is located in the right part of the Grandate plain (Figure 8) has been modelled with a lower value of hydraulic conductivity compared to the surrounding zones and a recharge value equal to half of the one that was assigned to the green areas.The correct representation of this zone is fundamental for modelling the system behavior given that part of the water that flows toward the plain arrives from the surrounding hills.Figure 9 shows the monthly rainfall that was registered in the 2010-2016 period in Como, from which the infiltration in non-urbanized areas (8•10 −9 m/s) has been evaluated using the Thornthwaite method and assuming a 30% runoff.For the urbanized areas, the recharge has been set over each municipality as 15% of the public withdrawal, then taking into account the network losses [41].Where the withdrawal data was not available (Luisago and Casnate con Bernate municipalities), a usage of 160 L/year per person has been considered and then 15% of this value was used for the evaluation of recharge.Urbanized and non-urbanized areas are shown in Figure 8.The morainic arc that is located in the right part of the Grandate plain (Figure 8) has been modelled with a lower value of hydraulic conductivity compared to the surrounding zones and a recharge value equal to half of the one that was assigned to the green areas.The correct representation of this zone is fundamental for modelling the system behavior given that part of the water that flows toward the plain arrives from the surrounding hills.Figure 9 shows the monthly rainfall that was registered in the 2010-2016 period in Como, from which the infiltration in non-urbanized areas (8•10 −9 m/s) has been evaluated using the Thornthwaite method and assuming a 30% runoff.For the urbanized areas, the recharge has been set over each municipality as 15% of the public withdrawal, then taking into account the network losses [41].Where the withdrawal data was not available (Luisago and Casnate con Bernate municipalities), a usage of 160 L/year per person has been considered and then 15% of this value was used for the evaluation of recharge.Urbanized and nonurbanized areas are shown in Figure 8.The calibration was initially performed in steady state using the inverse procedure through PEST [42] with reference to the piezometric survey that was carried out in November 2015.The calibration parameter was the hydraulic conductivity, the initial value of which was derived from a few hydraulic tests that were performed in some of the wells that were located in the model domain during the drilling phase.The available data have been interpolated in order to generate the initial hydraulic conductivity distribution for the four layers.The spatially distributed parameter was then the object of the automatic calibration process.
Once the model was calibrated, the attention was focused on the area where flooding occurred in 2014 for the evaluation of the feasibility of groundwater level control through pumping wells.Some wells were then implemented in the model and were optimized, first using PEST to minimize their number and pumping rate, which was then done manually in order to set their position.The goal of the optimization was to maintain the groundwater depth around 9 m below the ground (level constraint), in steady state condition.This level was defined considering a safety distance from the underground structures, which in this zone generally reaches a depth of 4.5-5 m from the ground.A global objective function was constructed in order to take into account the level constraint and the minimization of the pumping rate.This second objective was obtained allowing PEST to set zero as the minimum possible flow rate per well and thereby, at the same time, to optimize the number of wells (which was constrained to a maximum of 10).The position of the wells was then manually arranged according to the area that was effectively available for their installation.The mathematical model was then calibrated in transient mode as well, using the hydraulic heads that were registered in the monitoring points, and was finally used to forecast the system behaviour in three different scenarios, namely Scenario1 (2014 rainfall and no interventions), Scenario2 (2014 rainfall and continuous extraction), and Scenario3 (2014 rainfall and three months extraction).The unsteady calibration was performed taking into account 1 year piezometric data that was collected during the project development (November 2015-November 2016) and has been carried out to define the effective porosity value and the recharge term for the first layer.During this second calibration phase, the still calibrated hydraulic conductivity distribution was not modified.

Stochastic Analysis
The year 2014 has been recognized everywhere in Italy as a very exceptional year in terms of rainfall (100 years return time), as is possible to notice from Figure 8. Anyway, during this study, the numerical model was also used to perform a stochastic analysis with the aim of assessing the probability that flooding, such as what occurred in 2014, could happen again in the future.For this analysis, the available rainfall data in a meteorological station that was located close to Grandate was collected (20 years).It is possible to see in Figure 10a that the years 2014 and 2010 represent two peaks in the graph, with an annual cumulative precipitation that exceeds 2000 mm.It is also interesting to note the great difference between these years and the following ones (2011 and 2015, respectively), which are among the lowest of the entire twenty-year period in terms of rainfall amount.However, in 2010, the problems that occurred in 2014 did not happen.This could be explained by looking at the average monthly rainfall trend (Figure 9), showing that most of the rainfall that was observed in 2010 was recorded in warm or relatively hot months, such as May, August, and September when the evapotranspiration is high and, consequently, there is low infiltration of rainwater in the ground.Secondly, it is important to consider the process of deindustrialization that the area has faced: in 2010, just two years after the start of the 2008 economic crisis, the number of companies and production activities were higher than in 2014, the extraction from the aquifer was accordingly more intense, and the groundwater was at a greater depth.
The probability distribution of the collected rainfall data has been compared with the theoretical one [43] in order to verify whether, as expected, the data could be represented by a Gaussian distribution (Figure 10b).The slight disagreement between the data and theoretical curve occurs just in the tails of the distribution (representative only of the extreme values) and could be considered acceptable for the purposes of this analysis.Starting from the monthly mean and standard deviation of the normally distributed rainfall data, 100 random rainfall scenarios were generated through the PEST RANDPAR utility [42].This utility allows the generating of a series of random values (rainfall values in this case) that are distributed according to a probability distribution that is characterized by a determined mean and standard deviation.In this case, the 100 rainfall values therefore prove to be independent and normally distributed.The rainfall data that was thus generated have been randomly sampled in order to create 100 annual rainfall scenarios from which the recharge to be assigned to the stress periods of 100 transient model runs has been evaluated.Among the 100 scenarios, there are three hypothetical years that referred to very high rainfall, strong drought, and medium conditions.These extreme scenarios have been created by collecting for each month the 10th, 50th, and 90th percentile over the 100 generated rainfall data.Thereafter, the model was used to run 100 1-year long transient simulations considering the probability distribution of recharge related to the 100 randomly generated rainfall scenarios (including the three extreme scenarios).The transient simulations (1 year) resulted in different piezometric level maps.By collecting the piezometric head distributions that resulted from the simulations, monthly probability curves of groundwater exceeding a threshold level have been obtained.The threshold level is the one of the underground buildings in the flooded area (about 5 m below the ground).The initial head distribution that was used in these simulations is the one that is related to the November 2015 monitoring campaign that can be considered as representative of the present piezometric setting.As it is possible to see from Figure 5, the current average groundwater level in the flooded area is about 299 m a.s.l., which is similar to the value that was measured in November 2015.The analysis was thus intended to assess what the probability could be for groundwater to exceed the defined threshold starting from the current piezometric setting.

Numerical Model
The model was calibrated at the beginning in steady state and the general statistics are reported in Table 1.The model results are shown in Figure 11a where the residual values of hydraulic head are also indicated.Figure 11b shows the scatterplot of the observed vs. simulated heads.Calibration was performed using the piezometric levels of the November 2015 survey assuming that they could be representative of the average piezometric distribution in the area since no other useful data were available.The automatic calibration process with Pilot Point technique (PEST, [44]) was applied to Starting from the monthly mean and standard deviation of the normally distributed rainfall data, 100 random rainfall scenarios were generated through the PEST RANDPAR utility [42].This utility allows the generating of a series of random values (rainfall values in this case) that are distributed according to a probability distribution that is characterized by a determined mean and standard deviation.In this case, the 100 rainfall values therefore prove to be independent and normally distributed.The rainfall data that was thus generated have been randomly sampled in order to create 100 annual rainfall scenarios from which the recharge to be assigned to the stress periods of 100 transient model runs has been evaluated.Among the 100 scenarios, there are three hypothetical years that referred to very high rainfall, strong drought, and medium conditions.These extreme scenarios have been created by collecting for each month the 10th, 50th, and 90th percentile over the 100 generated rainfall data.Thereafter, the model was used to run 100 1-year long transient simulations considering the probability distribution of recharge related to the 100 randomly generated rainfall scenarios (including the three extreme scenarios).The transient simulations (1 year) resulted in different piezometric level maps.By collecting the piezometric head distributions that resulted from the simulations, monthly probability curves of groundwater exceeding a threshold level have been obtained.The threshold level is the one of the underground buildings in the flooded area (about 5 m below the ground).The initial head distribution that was used in these simulations is the one that is related to the November 2015 monitoring campaign that can be considered as representative of the present piezometric setting.As it is possible to see from Figure 5, the current average groundwater level in the flooded area is about 299 m a.s.l., which is similar to the value that was measured in November 2015.The analysis was thus intended to assess what the probability could be for groundwater to exceed the defined threshold starting from the current piezometric setting.

Numerical Model
The model was calibrated at the beginning in steady state and the general statistics are reported in Table 1.The model results are shown in Figure 11a where the residual values of hydraulic head are also indicated.Figure 11b shows the scatterplot of the observed vs. simulated heads.Calibration was performed using the piezometric levels of the November 2015 survey assuming that they could be representative of the average piezometric distribution in the area since no other useful data were available.The automatic calibration process with Pilot Point technique (PEST, [44]) was applied to the hydraulic conductivity parameter and resulted in a hydraulic conductivity distribution for each layer.The hydrogeological study and the model runs demonstrated that the reason for the 2014 flooding was mainly due to the concurrence of three causes: (1) the hydrogeological structure of the area recognized as a stagnation zone, (2) the groundwater rising-which has been happening for several years in this area and generally in Lombardy [43,45]-and (3) the abnormal quantity of rainfall that occurred in 2014 (100-years return period).
Once the model was calibrated, the attention was focused on the area where flooding occurred in 2014 to assess the feasibility of groundwater level control through pumping wells.The optimization process result shows that 8 wells are necessary to achieve the goal, with pumping rates varying between 3 and 6 L/s for a total of about 40 L/s in the area (Figure 12).The action of pumping wells close to the areas that mainly experienced the flooding problem could thus be an effective system to keep groundwater levels under control and to avoid the repetition of such a phenomenon.The hydrogeological study and the model runs demonstrated that the reason for the 2014 flooding was mainly due to the concurrence of three causes: (1) the hydrogeological structure of the area recognized as a stagnation zone, (2) the groundwater rising-which has been happening for several years in this area and generally in Lombardy [43,45]-and (3) the abnormal quantity of rainfall that occurred in 2014 (100-years return period).
Once the model was calibrated, the attention was focused on the area where flooding occurred in 2014 to assess the feasibility of groundwater level control through pumping wells.The optimization process result shows that 8 wells are necessary to achieve the goal, with pumping rates varying between 3 and 6 L/s for a total of about 40 L/s in the area (Figure 12).The action of pumping wells close to the areas that mainly experienced the flooding problem could thus be an effective system to keep groundwater levels under control and to avoid the repetition of such a phenomenon.
rainfall that occurred in 2014 (100-years return period).
Once the model was calibrated, the attention was focused on the area where flooding occurred in 2014 to assess the feasibility of groundwater level control through pumping wells.The optimization process result shows that 8 wells are necessary to achieve the goal, with pumping rates varying between 3 and 6 L/s for a total of about 40 L/s in the area (Figure 12).The action of pumping wells close to the areas that mainly experienced the flooding problem could thus be an effective system to keep groundwater levels under control and to avoid the repetition of such a phenomenon.

Stochastic Analisys
Through collecting the head distributions that resulted from the 100 simulations that were run during the stochastic analysis that was carried out, a monthly probability curve of groundwater exceeding a threshold level was obtained.Figure 13 shows the probability distribution of the hydraulic head for some months of the simulated year compared to the elevation of the underground buildings in order to point out the hazard for the underground buildings (i.e., parking area, cellars, and building foundations) to be flooded.The curves resulting from the stochastic analysis show that, in the area, there is an occurrence probability of exceeding the groundwater level of the underground structures between 12% and 15%, respectively in November and December.This means that underground flooding is likely to occur once approximately every 10 years.Thus, even if 2014 was a very abnormal year in terms of rainfall, the underground buildings are exposed to a not negligible (equal to about 0.12-0.15)occurrence probability to be flooded in the future (considering no changing of groundwater extraction conditions compared to the simulated situation).

Stochastic Analisys
Through collecting the head distributions that resulted from the 100 simulations that were run during the stochastic analysis that was carried out, a monthly probability curve of groundwater exceeding a threshold level was obtained.Figure 13 shows the probability distribution of the hydraulic head for some months of the simulated year compared to the elevation of the underground buildings in order to point out the hazard for the underground buildings (i.e., parking area, cellars, and building foundations) to be flooded.The curves resulting from the stochastic analysis show that, in the area, there is an occurrence probability of exceeding the groundwater level of the underground structures between 12% and 15%, respectively in November and December.This means that underground flooding is likely to occur once approximately every 10 years.Thus, even if 2014 was a very abnormal year in terms of rainfall, the underground buildings are exposed to a not negligible (equal to about 0.12-0.15)occurrence probability to be flooded in the future (considering no changing of groundwater extraction conditions compared to the simulated situation).

Discussion
The flooding event that occurred in Grandate in November 2014 pointed out the need to understand the hydrogeological system in this area and to conceive drainage interventions to lower and control groundwater levels.

Discussion
The flooding event that occurred in Grandate in November 2014 pointed out the need to understand the hydrogeological system in this area and to conceive drainage interventions to lower and control groundwater levels.
The mathematical model, which was calibrated in an unsteady state as well as using the levels that were registered in the piezometer that was located in the flooded area since March 2016, was used to forecast the system behavior in three different scenarios: Scenario 1 (Business As Usual: 2014 rainfall, no interventions are simulated)-Starting from the calibrated situation, the recharge values of 2014 were set in the model and the year was simulated.Thus, the model was used to analyze the behaviour of the system in case a very rainy year, such as the year 2014-listed among the first 10 rainiest years of the last century-should occur again.Figure 14a shows the model results at the end of the month of November in terms of flooded cells at ground level, 4.5 m depth and 8 m depth.The mathematical model, which was calibrated in an unsteady state as well as using the levels that were registered in the piezometer that was located in the flooded area since March 2016, was used to forecast the system behavior in three different scenarios: Scenario 1 (Business As Usual: 2014 rainfall, no interventions are simulated)-Starting from the calibrated situation, the recharge values of 2014 were set in the model and the year was simulated.Thus, the model was used to analyze the behaviour of the system in case a very rainy year, such as the year 2014-listed among the first 10 rainiest years of the last century-should occur again.Figure 14a shows the model results at the end of the month of November in terms of flooded cells at ground level, 4.5 m depth and 8 m depth.Scenario 2 (2014 rainfall and continuous extraction)-The same situation of Scenario 1 is simulated, however now the 8 wells are implemented in the model, pumping out continuously throughout the simulation time.The result is set out in Figure 14b, which shows the flooded cells at ground level, 4.5 m depth and 8 m depth at the end of the month of November.As expected, the well pumping keeps groundwater level below the 4.5 m from ground level and the area is just partially flooded at 8 m depth.It should be recalled that a very peculiar year in terms of rainfall (and then recharge) is simulated.
Scenario 3 (2014 rainfall and 3 months extraction)-Scenario 2 has been replicated, but here the wells do not work continuously: their activation was presumed in August, about three months before the most critical period in terms of rainfall which, in these zones, is generally November.In fact, from the perspective of the groundwater control systems, to optimize the extraction and to carry it out close to the most critical period might prove convenient in terms of energy consumption and extracted groundwater disposal.The result of this last simulation is shown in Figure 14c, setting out the model response at the end of the 11th month of the simulation (November).It is still possible to notice the effect of the pumping, which makes it possible to maintain groundwater level below 4.5 m from ground level.The wells' effect is obviously reduced at an 8 m depth from the ground compared to Scenario 2, given that the wells were active for just three months.However, as has already been said, a very exceptional year in terms of rainfall was simulated.

Conclusions
In 2014, the Grandate plain experienced the problem of groundwater flooding of several underground factory buildings.The groundwater table rise can interfere with the underground structures, with obvious inconveniences and high economic damages.This happens especially in particularly exposed areas, such as the depressions that in the recent past have been at the site of lakes, swamps, and peat bogs.Grandate is an example of such areas: its hydrogeological setting, similar to that of extensive parts of the Lombard territory, makes this zone representative of a large number of areas that are affected by similar problems.The hydrogeological characterization that was carried out in Grandate after the flooding event and the numerical model that was implemented have demonstrated that the reason for the floods was the concurrence of three main causes: (1) the hydrogeological setting of the area was recognized as a stagnation zone, (2) the groundwater rising due to the process of deindustrialization that was happening almost everywhere in Lombardy, and (3) the abnormal quantity of rainfall that occurred in 2014 (100-years return period).The geological setting of the Grandate plain, representative of the Lombard morainic amphitheatres, makes the area particularly susceptible to groundwater flooding.The Grandate area has in fact been recognized as a stagnation area and its conformation and surroundings cause water to accumulate in the subsoil and flow with difficulty southward.Because of this geological conformation, a strong correlation between rainfall events and groundwater level increase has been observed, especially when intense rainfall occurrences happen, to the extent of those that were registered in 2014.Considering the variability of rainfall and, consequently, the infiltration rate influencing the groundwater table, a stochastic analysis was carried out in order to define probability exceeding curves in the critical area.Results confirm that the probability for groundwater to exceed the underground levels in the current piezometric setting is not negligible (12-15%), and the numerical model proves to be essential to quantify the effects of groundwater drainage that are operated by pumping wells in such a complex area.In order to control and forecast the system behavior, the monitoring of groundwater levels in real time by remote sensors is one of the first requirements.Moreover, a WebGis tool could be suitable for keeping the groundwater level variation under control and giving indications about the groundwater approaching an exceeded threshold level, whereupon a pre-alarm could be given and/or the pumping wells system could be activated.This kind of control of groundwater could be adequate for avoiding flooding occurrences in the future.

17 Figure 1 .
Figure 1.Study area (red square) and the study area's location; Grandate municipality is contoured in black.

Figure 1 .
Figure 1.Study area (red square) and the study area's location; Grandate municipality is contoured in black.

17 Figure 2 .
Figure 2. (a) Digital Elevation Model of the study area, principal hydrographic network, and model domain (red square); (b) stratigraphic logs (yellow dots), and realized cross-section traces (in black Section 4 and C).

Figure 2 .
Figure 2. (a) Digital Elevation Model of the study area, principal hydrographic network, and model domain (red square); (b) stratigraphic logs (yellow dots), and realized cross-section traces (in black Section 4 and Section C).

Figure 3 .
Figure 3. Cross Section 4 and C passing through the Grandate plain.Figure 3. Cross Section 4 and Section C passing through the Grandate plain.

Figure 3 .
Figure 3. Cross Section 4 and C passing through the Grandate plain.Figure 3. Cross Section 4 and Section C passing through the Grandate plain.

Figure 4 .
Figure 4. Piezometric map of the first aquifer and measured points with their groundwater level (head expressed in m a.s.l.).Encircled in red is the zone that was flooded in 2014.Piezometric lines 305 m a.s.l. and 295 m a.s.l.mark the low gradient zone (stagnation area).

Figure 5 .
Figure 5. Monitoring data for the two points located in the flooded zone.

Figure 4 .
Figure 4. Piezometric map of the first aquifer and measured points with their groundwater level (head expressed in m a.s.l.).Encircled in red is the zone that was flooded in 2014.Piezometric lines 305 m a.s.l. and 295 m a.s.l.mark the low gradient zone (stagnation area).

Figure 4 .
Figure 4. Piezometric map of the first aquifer and measured points with their groundwater level (head expressed in m a.s.l.).Encircled in red is the zone that was flooded in 2014.Piezometric lines 305 m a.s.l. and 295 m a.s.l.mark the low gradient zone (stagnation area).

Figure 5 .
Figure 5. Monitoring data for the two points located in the flooded zone.

Figure 5 .
Figure 5. Monitoring data for the two points located in the flooded zone.

Figure 6 .
Figure 6.Example of a stagnation area (arrows represent flow lines, the light-blue line represents the phreatic surface, and the black lines represent the vertical head contours).

Figure 7 .
Figure 7. Horizontal (a) and vertical (b) discretization (in the plant, red lines indicate the sections traces, the River condition is blue, while light blue represents GHB conditions contouring the model).

Figure 6 .
Figure 6.Example of a stagnation area (arrows represent flow lines, the light-blue line represents the phreatic surface, and the black lines represent the vertical head contours).

Figure 7 .
Figure 7. Horizontal (a) and vertical (b) discretization (in the plant, red lines indicate the sections traces, the River condition is blue, while light blue represents GHB conditions contouring the model).

Figure 7 .
Figure 7. Horizontal (a) and vertical (b) discretization (in the plant, red lines indicate the sections traces, the River condition is blue, while light blue represents GHB conditions contouring the model).

Figure 8 .
Figure 8. Urbanized (colored zones) and non-urbanized (dotted zones) areas in the model domain.The morainic arc in the right part of Grandate is contoured in black.

Figure 8 .
Figure 8. Urbanized (colored zones) and non-urbanized (dotted zones) areas in the model domain.The morainic arc in the right part of Grandate is contoured in black.

Figure 8 .
Urbanized (colored zones) and non-urbanized (dotted zones) areas in the model domain.The morainic arc in the right part of Grandate is contoured in black.

Figure 9 .
Figure 9. Monthly rainfall that was registered over 7 years in a meteorological station in Como.

Figure 9 .
Figure 9. Monthly rainfall that was registered over 7 years in a meteorological station in Como.

17 Figure 10 .
Figure 10.20 years of available rainfall data (a) and a comparison between the rainfall distribution and the Gaussian one (b).

Figure 10 .
Figure 10.20 years of available rainfall data (a) and a comparison between the rainfall distribution and the Gaussian one (b).

Figure 11 .
Figure 11.(a) Model simulation result and observed-simulated values-zoom on Grandate, (b) observed/simulated scatterplot before and after the calibration procedure.

Figure 11 .
Figure 11.(a) Model simulation result and observed-simulated values-zoom on Grandate, (b) observed/simulated scatterplot before and after the calibration procedure.

Figure 12 .
Figure 12.(a) The area that was subjected to the flooding in 2014 (red grid) and pumping wells (blue circles) inserted into the model; (b) line of 9 m groundwater depth that was created by the action of the 8 wells in the November 2015 recharge condition.The line outlines the internal cells of the grid where the groundwater depth is 9 m or even bigger.

Figure 12 .
Figure 12.(a) The area that was subjected to the flooding in 2014 (red grid) and pumping wells (blue circles) inserted into the model; (b) line of 9 m groundwater depth that was created by the action of the 8 wells in the November 2015 recharge condition.The line outlines the internal cells of the grid where the groundwater depth is 9 m or even bigger.

Figure 13 .
Figure 13.Probability distribution of groundwater levels for some months of the simulated year.The black line represents the buildings' underground level (5 m below the ground) in the flooded zone.

Figure 13 .
Figure 13.Probability distribution of groundwater levels for some months of the simulated year.The black line represents the buildings' underground level (5 m below the ground) in the flooded zone.

Figure 14 .
Figure 14.Results of Scenarios 1 (a), 2 (b), and 3 (c); wells (circles) and flooded zones (grey cells) at different depths.Scenario 2 (2014 rainfall and continuous extraction)-The same situation of Scenario 1 is simulated, however now the 8 wells are implemented in the model, pumping out continuously throughout the simulation time.The result is set out in Figure 14b, which shows the flooded cells at ground level, 4.5 m depth and 8 m depth at the end of the month of November.As expected, the well

Table 1 .
General statistics resulting from the model.

Table 1 .
General statistics resulting from the model.