## 2 Materials and Methods
### Mathematical Model  
The conceptual site model has been implemented in Modflow2005 [39] throughout an 80 \(\times\) 80 m grid and was refined up to 40 \(\times\) 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 Wurm 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]. 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). 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\(\cdot\)10\({}^{-9}\) m/s) has been evaluated using the Thorinthwaite 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.  
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 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). 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.