Data-driven decision management of urban underground infrastructure through groundwater-level time-series cluster analysis: the case of Milan (Italy)

The significant increase in urbanization has resulted in greater use of the subsurface in urban planning and, therefore, increased interaction between groundwater and underground infrastructure. Numerical models are the primary tool adopted to manage the resulting problems; however, their construction is time- and cost-consuming. Groundwater-level time-series analysis can be a complementary method, as this data-driven approach does not require an extensive understanding of the geological and boundary conditions, even if providing insights into the hydrogeologic behaviour. Thus, a data-driven approach was adopted to analyse groundwater time-series of the shallow aquifer, occupied by several underground structures, beneath Milan city (Northern Italy). Statistical (Mann-Kendall and Sen’s slope estimator, autocorrelation and cross-correlation, hierarchical cluster analysis) and geospatial techniques were used to detect the potential variables influencing the groundwater levels of 95 monitoring wells, covering the period 2005–2019. A general rising trend of the water table was identified, with local hydrogeologic differences in the western and southernmost areas. Based on time-series analysis results, four management areas have been identified. These areas could act as future geographic units with specific groundwater management strategies. In particular, subsurface public car parks can be classified with respect to groundwater flooding as (1) not submerged, (2) possibly critical, or (3) submerged at different groundwater conditions. According to these outcomes, targeted guidelines for constructing new car parks have been elaborated for each management area. The methodology proved to be efficient in improving the urban conceptual model and helping stakeholders design the planned underground development, considering groundwater aspects.


Introduction
Groundwater represents the most valuable source of freshwater on a global scale (Li et al. 2013a), as it is used for drinking, agricultural and industrial purposes (Filimonau and Barth 2016). As the world population is projected to grow in the next decades, increasing demand for this resource is expected, and this may provoke stress both in terms of groundwater quality and quantity (Pollicino et al. 2021); in particular, groundwater demand is gradually increasing in urban areas, where 70% of the world population is expected to be living by 2050 (Un-Habitat 2012). To limit urban sprawl, the tendency to look for new building spaces in the vertical direction has increased (Li et al. 2013b, c;Koziatek and Dragićević 2017), leading to enhanced use of urban underground (Parriaux et al. 2007;Sterling et al. 2012;Bobylev 2016;Vähäaho 2016). In the last few decades, as a result of the deindustrialization process, several cities worldwide have faced rising groundwater levels, which has generated some interference between underground infrastructure (i.e. basements, car parks, and subway lines) and groundwater (Wilkinson 1985;George 1992;Hernández et al. 1997;Vázquez-Suñé et al. 1997;Mudd et al. 2004;Hayashi et al. 2009;Lamé 2013;Ducci and Sellerino 2015). The groundwater/infrastructure interaction needs to be properly managed to ensure a sustainable use of urban 1 3 underground spaces. Numerical groundwater flow modelling can be used for evaluating the interference between groundwater and underground structures (Vázquez-Suñé et al. 2005;García-Gil et al. 2015) and for identifying possible management strategies. A variety of three-dimensional (3D) numerical models has been realized worldwide (Attard et al. 2015); however, the construction of a 3D numerical model is generally a multistep procedure requiring detailed knowledge of the aquifer system, both from the geological and hydrological features and, thus, such models take a considerable amount of time to develop. As the information is sometimes lacking or scattered among different sources, the conceptual model also needs to be frequently revised, which could be considered costly and time-consuming (Vázquez-Suñé et al. 2005). Within this scope, data-driven techniques can constitute a valid complementary approach, providing preliminary insights into the main drivers affecting the aquifer behaviour (Chae et al. 2010;Obergfell et al. 2013), without requiring an extensive approximation of subsurface parameters and boundary conditions. Generally, adopting data-driven methods is faster and more easy-applicable than numerical models (Bakker and Schaars 2019) and could represent a reasonable compromise between model complexity and costs (Stevenazzi et al. 2017); however, time-series analysis is indispensable to support groundwater model calibration (Bakker and Schaars 2019). In particular, time-series analysis of groundwater levels allows one to detect trends (Patle et al. 2015), and thereby identify the causes of water level variations through the investigation of input-output connections (Obergfell et al. 2013); for these reasons, this approach can be adopted to plan effective groundwater management strategies, thus bringing valuable information to planners and decision-makers (Best et al. 2021). Also, as modelling should make maximum use of existing data, by assimilating its information content (Lotti et al. 2021), data-driven methods could support the construction of a three-dimensional (3D) numerical model. Different time-series analysis techniques have been applied to interpret groundwater systems processes. Statistical tests such as the Mann-Kendall test and Sen's slope estimator have been widely applied to detect long-term trends from groundwater level series, for evaluating their significance and magnitude (Tabari et al. 2012;Singh et al. 2015;Meggiorin et al. 2021). Furthermore, statistical tests have been sometimes combined with geospatial approaches in a geographical information system (GIS) environment to identify those areas of groundwater-level rise or drawdown, determined by the existence of a trend (Chaudhuri and Ale 2014;de Brito Neto et al. 2016;Anand et al. 2020). On the other hand, clustering techniques have been mainly adopted for groundwater quality, but recently different studies have applied these techniques to groundwater level time series (Pathak and Dodamani 2019;Rinderer et al. 2019;Naranjo-Fernández et al. 2020;Sottani et al. 2020).
The present study aimed at proposing, for the city of Milan (Lombardy, Northern Italy), a data-driven approach to support decision making for urban underground infrastructure in relation to groundwater aspects. In this way, the stakeholders could properly manage the future underground development of the city in a more targeted way, also considering hydrogeologic features. To the best of the authors' knowledge, this is the first time that groundwater/ underground infrastructure interactions have been managed through the support of data-driven methods. A set of statistical techniques has been applied to groundwater-level time series of the shallow aquifer in Milan, with two main purposes: (1) to identify different spatial zones with local hydrogeologic peculiarities, possibly acting as future groundwater management territorial units, and to investigate the potential features influencing groundwater levels (i.e. geology, precipitation, withdrawals, land use); and (2) to elaborate targeted guidelines for underground infrastructure management, mainly focusing on the constructions of new public car parks. Underground car parks have been mainly considered, as their interference with the groundwater has been only partially investigated. Nonetheless, this procedure could be applied to every category of underground infrastructure.
The data-driven analysis proposed in this study has been fine-tuned for the Milan metropolitan area; notwithstanding, it could be applied to other urban areas worldwide that have subsurface infrastructure development and a complex shallow aquifer system that requires an integrated urban groundwater management strategy.

Hydrogeological conceptual model of the study area
The study area (182 km 2 ) is the city of Milan (Fig. 1). It hosts 1.4 million people (Istat 2011) and it is characterized by intense industrial and agricultural development (Bonomi 2009). The land use in the western and southern areas is still mainly agricultural, while the remaining areas of the city are mainly urbanized. A wide hydrographic network characterizes the area (Fig. 1b): among the main rivers, Olona and Seveso are culverted within the city, while the Naviglio Grande and Naviglio Pavese are man-made surface channels mainly used for irrigation. Milan is located in the middle of the Po Plain (Fig. 1a), in proximity to the lowland springs belt and lies on a sedimentary aquifer system whose hydrogeologic structure has been deeply studied in past years (Regione Lombardia and ENI 2002). A recent classification (Regione Lombardia 2016) has identified three main hydro structures: a shallow hydro structure (ISS), an intermediate hydro structure (ISI) and a deep hydro structure (ISP). In the study area, ISS has a medium thickness of 50 m, with a bottom surface ranging from 100 m above sea level (asl; to the north) down to around 50 m asl (to the south). It hosts a shallow aquifer- Fig. 2; Fig. S1 of the electronic supplementary material (ESM-not exploited for drinking needs. This shallow aquifer is the focus of this study and is mainly composed of sands and gravels. The same lithology, but with an increasing presence of silty and clayey layers, characterizes the ISI. Finally, because of a limited amount of data, the ISP has an uncertain lithological composition. Groundwater has been extensively exploited for industrial use since the early 1960s; around 1975 there was a water-table depth of more than 30 m in the northern sector. Subsequently, since the early 1990s, the dismantling of the biggest factories, mainly located in the northern sector of the town, led the water table to rise once again (i.e., with a maximum rise of about 10-15 m in the northern area), causing flooding episodes for some nonwaterproofed underground structures (Cavallin and Bonomi 1997;Beretta et al. 2004;Gattinoni and Scesi 2017;Colombo et al. 2018), including car parks and the oldest, shallower subway lines (Fig. 1b). Thus, the most recent subway lines (Fig. 1b), also due to their greater depth, were designed with waterproof systems. Public car parks were built in sequential phases starting from the middle of the 1980s. Now, 126 underground public car parks are listed in the territory ( Fig. 1b; Table S1 of the ESM). The information (i.e., bottom elevation (m asl), period of construction, location) on subway lines, and public and private car parks is stored in a 3D geodatabase (Sartirana et al. 2020), allowing for proper management of the underground infrastructure in the framework of groundwater management (Fig. 2). Fig. 1 a Geographical setting of the study area; b main hydrogeologic features (lowland springs and surface-water bodies), nonwaterproofed and waterproofed subway lines, public car parks, groundwater heat pumps (GWHPs) and wells crossed by section BB′; line AA′ points to the location of the cross section that is visible in Fig. 2; line BB′ points to the location of the geological profile that is shown in Fig. S1 of the ESM. To obtain a clear representation, public car parks have been represented as points instead of polygons Moreover, the number of low-enthalpy geothermal systems in the shallow aquifer increased in the last decades due to regional laws in the scope of European goals to reduce greenhouse gas emissions. Particularly, open-loop groundwater heat pumps (GWHPs) serve large building systems, mainly in the downtown area (Previati and Crosta 2021a) (Fig. 1b). The average depth of GWHP filters ranges between 17 and 34 m (Regione Lombardia 2021a). In some cases, GWHPs discharge exploited water to surface-water bodies to limit the groundwater rise.
The adoption of the Plan of Government for the Territory (PGT) (Milan Metropolitan City 2019) in 2019 will favour further subsurface occupation by developing new subway line sections and reducing soil consumption with consequent changes in the subsurface. Hence, the need for coordinated management among the stakeholders arises, supported by evidence from shallow-aquifer studies.

Materials and methods
The methodology proposed in this study is composed of 7 steps: (1) dataset pre-processing, (2) groundwater timeseries trend analysis, including the Mann-Kendall test, Sen's slope estimator and a time-series reversal points assessment, (3) correlation analysis, through autocorrelation function (ACF) and cross-correlation function (CCF) with precipitation data; (4) groundwater time-series clustering, applying hierarchical cluster analysis (HCA); (5) water-table reconstruction and groundwater volume calculation under specific groundwater conditions identified at point 2; based on the outcomes of points 1-5, the subsequent steps were (6) identification of groundwater management areas (MAs); and (7) definition of guidelines for underground infrastructure management. A graphical representation of the adopted workflow is given in Fig. 3.

Available data and preprocessing
Groundwater level data were provided by Metropolitana Milanese Spa, the local water supplier. The original dataset spans from 1950 to 2019, but data are scattered over ca. 300 monitoring wells (MWs) changing over time due to the monitoring network evolution (Fig. 1b). Excluding MWs not having a useful data resolution, a subset of the most consistent time series was identified for the time-series analysis, consisting of 95 MWs (Fig.1b)  Two measures per year are available (March and September; with 30 observations in total). Once the final dataset had been defined, it was checked for possible errors or outliers (Azzellino et al. 2019). Outlier analysis was performed through the nonparametric interquartile range test, since the Shapiro-Wilk test highlighted that the data distribution was not normal. Errors were corrected manually only if clearly identifiable in the original data (i.e. ground level values), after verifying that they could offer insights into some hydrogeologic process that should not be neglected (Peterson et al. 2018).   As for precipitation, the data of Paderno Dugnano rain gauge (located just northward of the city of Milan) were considered, monitored by the regional environmental protection agency (ARPA Lombardia 2021).

Mann-Kendall test and Sen's slope estimator
The Mann-Kendall test (Mann 1945;Kendall 1948) is a rank-based, nonparametric test widely used for trend detection (i.e. upward or downward trend) in hydrological variables since it is nonparametric and robust to extreme data Anand et al. 2020;Meggiorin et al. 2021). The null hypothesis (H 0 ) is the absence of a trend, while the alternative hypothesis (H 1 ) is the presence of a significant trend. A confidence level of 95% was used as the threshold value.
The Sen's slope estimator (Sen 1968) was applied to estimate the magnitude of the analysed trend. It is a nonparametric unbiased estimator of trends, more robust to outliers than a linear regression estimator (Hirsch et al. 1982). A positive value of Sen's slope indicates an increasing trend, whereas a negative one indicates a decreasing trend.

Time-series reversal points
The reversal points analysis for the time series aimed at identifying when an upward (or downward) trend has been reversed. This could be useful in groundwater management, as it allows one to precisely comprehend the dynamics of the water table in a given period.
The Pettitt test (Pettitt 1979) is a widely applied statistical tool to detect an abrupt change in a dataset. Despite its common use, it seems not so adaptable when trying to identify a pure reversal point (i.e. a shift from an upward trend to a downward trend), at least in small datasets (Frollini et al. 2021).
Therefore, reversal points were identified here as the time-series point determining a change in the sign of the difference between successive couples of data, previously smoothed through a loess function to avoid noise effect. This analysis was performed in the R environment, using the imputeTS package (Moritz and Bartz-Beielstein 2017), ts function (Becker et al. 1988) and stl function (Cleveland et al. 1990).
In this way, only the most recurrent groundwater-level minimum and maximum conditions among the different MWs have been investigated. An additional description of their subsequent application is given in section "Water table reconstruction and groundwater volumes calculation". Further details on the adopted approach and the R script are available in the ESM.

Autocorrelation and cross-correlation
Autocorrelation is the correlation of a time series with itself at different time lags. It indicates whether the past values of a time-series influence the current state. Here, autocorrelation analysis was performed for all the MWs through the R Package statistics (R Core Team 2021) to investigate at the study scale the "memory effect" of the groundwater system, which can be considered as the requested amount of time to get rid of the dependency from a previous condition (Tamburini and Menichetti 2020).
The cross-correlation is a measure of similarity between two time series as a function of the displacement of one with respect to the other. To evaluate a possible delayed dependency at different time lags of groundwater levels (1) Water 1 3 (output series) on rainfall data (input series), a cross-correlation analysis between these variables was conducted. Thus, rainfall data were aggregated at a 6-month time step (cumulative rainfall periods: October-March/April-September) to keep equal data length while correlating the two variables and to avoid spurious conclusions that would arise from not using equal data records (Kumar et al. 2009).

Groundwater time-series clustering
Hierarchical cluster analysis (HCA) was performed, considering groundwater level measurements at different times as different variables, using the Ward hierarchical method (Ward Jr 1963).
The squared Euclidean distance was adopted as a measure of similarity between samples (Triki et al. 2014;Bloomfield et al. 2015) to compare the groundwater levels of different MWs at individual seasons, without associating the groundwater levels of different seasons, which could have happened when applying other distances typical of the time-series analysis.
In order to analyse both spatial and temporal patterns over a time-series, not only the initial groundwater measurements, but also seasonal differences between couples of successive observations have been considered. In both cases, to guarantee an equal weight for each considered variable in the calculation of the Euclidean distance matrix, HCA was carried out on standardized data (i.e. autoscaling; mean = 0 and variance = 1; Judd 1980). The HCA dataset, therefore, resulted in 95 rows (i.e. monitoring wells) and 59 variables (columns, i.e. groundwater level measures for March-September from 2005 to 2019 and seasonal variations between successive groundwater measures).
The HCA requires continuous data and no missing values. Therefore, missing data were filled, considering trend similarity among neighbouring MWs. First, a GIS analysis was carried out to identify the MWs nearest to those containing missing values. Missing values were then filled by calculating the average difference between the time series containing missing values and the most similar neighbouring MW; this average was calculated for the previous and the following observations surrounding the missing values, and then applied (summing/subtracting) to the date of the missing value. In this way, possible relevant peaks in the groundwater time series were conserved, avoiding any smoothing effect that could result from the time-series interpolation or other filling methods.

Water-table reconstruction and groundwater volumes calculation
Potentiometric maps for the shallow aquifer have been reconstructed, for the identified global minimum and maximum conditions (see section 'Groundwater minimum and maximum'), and for the last measure of the time series (September 2019, Sep19). Further piezometric levels, available for each specific period both inside and outside the study area (i.e. control points), have been employed to properly reconstruct the potentiometric maps.
As in the study area the presence of a spatial piezometric trend (NW-SE) has been detected (Bonomi et al. 1998;Beretta et al. 2004), groundwater heads in MWs have been interpolated using universal kriging (Isaaks and Srivastava 1989;Goovaerts 1997;Kitanidis 1997;Webster and Oliver 2001). Subsequently, the volume of water in the shallow aquifer, determined from the difference between the potentiometric heads in maps of the different groundwater conditions, has been calculated and the content of mobile water contained in that portion of the subsurface has been estimated considering the effective porosity (n e ). Reckoning the quantities of both occupied volumes and mobile water could in fact be deemed as a proper way of assessing the impact of groundwater on different resources, including manmade ones (Bonomi et al. 2010). To do so, geological information on the subsurface was required. The parametrisation of the hydrogeologic properties, such as hydraulic conductivity, has been readapted from a previous work of the research group on the same study area (Bonomi et al. 2010). Classes of effective porosity have been then extracted from the hydraulic conductivity values of the subsurface volumes derived from comparison of the water tables. This step was performed through the use of conversion conductivity-effective porosity tables, previously entered in a reference database (Bonomi et al. 2010(Bonomi et al. , 2014.

Decision management
On the basis of the results that emerged in previous steps, different spatial units (management areas, MAs) have been aggregated as a result of their hydrogeologic features (Chaudhuri and Ale 2013) in the context of a more targeted groundwater management strategy. Moreover, a deeper analysis associated with underground infrastructure, with a focus on public car parks, was elaborated through a GIS methodology previously applied by Sartirana et al. 2020, to evaluate the interaction between public car parks and the water table. This analysis was carried out on the identified groundwater-level minimum and maximum conditions and for Sep19. Based on these results, public car parks have been grouped into three categories with respect to groundwater flooding for each period: not submerged, possibly critical, submerged. A public car park was defined as "not submerged" or "submerged" if the water table was respectively below or above the bottom of the car park structure. The category "possibly critical" included all the public car parks showing a difference between the reference plan and the water table of less than 1 m. Engineering aspects, such as the possible uplift risks generated by buoyancy as a result of the aquifer pressure, are beyond the aims of this work. The identified possibility of submersion for each car park, therefore, facilitated a list of targeted management proposals, including suggested guidelines for the construction of new public car parks, based on the detected hydrogeologic evidence.

Mann-Kendall test and Sen's slope estimator
Trend analysis through the Mann-Kendall test and Sen's slope estimator showed a general increasing trend in the groundwater levels. Results are visible in Fig. 4. Particularly, at 95% confidence level, 86 out of 95 MWs showed an increasing trend throughout the considered period, with a medium slope ranging from 0.20 to 0.60 in the northern sector, from 0.10 to 0.20 in the central portion of the domain, and from 0.00 to 0.10 in the southern area. Only 9 out of 95 MWs did not show significant trends, all located in the westernmost portion of the study area, except for one in the southernmost section.

Groundwater minimum and maximum
Through the adopted approach developed in R, the global groundwater minimum and maximum conditions have been identified. Results are visible in Fig. 5; moreover, for a subset of selected MWs, the identified reversal points are visible, together with their groundwater trends in Fig. 6. The complete list of the identified reversal points for all the MWs is available in Table S2 of the ESM. As observable in Fig. 5, among the considered groundwater level time-series, September 2007 (Sep07) has been mainly identified as the global minimum condition; however, Sep07 is not identifiable as the historical global minimum, which took place at the end of the 1970s, as a consequence of intense industrial exploitation of groundwater (Bonomi 1999). Instead, March 2015 (Mar15) has been detected as the global maximum. As a general description for the study area, Sep07 comes at the end of a decreasing phase that began in December 2002, as described in a previous work about the same area (Sartirana et al. 2020). A rising phase, steeper in the northern sector, as confirmed by the Sen's slope estimator (Fig. 4), has been detected from Sep07 up to Mar15, while a subsequent declining phase took place from Mar15 up to Sep19 (final measurement of the dataset). Some MWs, mostly located in the peripheral areas, showed their groundwater minimum condition earlier than Sep07 (i.e. Mar07); in contrast, in the northernmost areas, the groundwater minimum was sometimes detected later than Sep07 (i.e. Mar08). Moreover, a higher number of reversal points has been identified for the MWs located in the peripheral sectors (eastern and western sides) and in the southern areas, with respect to the MWs placed in the northern and central parts of the domain (Fig. 6). Local minimum and maximum conditions have been also identified (Table S2 of

Autocorrelation and cross-correlation with rain data
Autocorrelation analysis was performed on the MWs considering 10 lags. For each MW, the number of significantly autocorrelated lags was determined and this is represented in Fig. 7. Results (Fig. 7a) show that autocorrelation can be significant for up to 6 lags in the northern sector, which correspond to a 3-year time span. In the central part of the area, significant autocorrelation was detected for a maximum of 4 lags (i.e. 2 years), while a more heterogeneous situation is visible in the southernmost part. A westernmost group of MWs shows no degree of correlation. Correlation of groundwater time series with rain data has been evaluated up to lag -2 (possible response time: 1 year), as these lags have been considered as the most environmentally significant; results are shown in Fig. 8.
At lag 0 (i.e. no temporal shift), the MWs located in the most peripheral areas showed a significant correlation with rain data, indicating a direct and fast response of piezometric level to precipitation (Fig. 8a) mainly in the western and the southernmost sectors, which showed a shorter autocorrelation and therefore a lower memory of the groundwater system. A significant but lower correlation was detected in the eastern sector of the domain. Thus, in these areas, groundwater levels respond rapidly to local precipitation events. At lag -1 (response time: 6 months) a more pronounced correlation between rain data and groundwater levels has been noticed in the northern and the central sectors of the domain (Fig. 8b), thus indicating a less rapid response of groundwater levels to rainfall; whereas, a reduced significant correlation has been identified in the western and southern sectors, that showed a significant correlation at lag 0. Finally, at lag -2 (response time: 1 year), a significant correlation with rain data was identified for 74 out of 95 MWs (Fig. 8c), thus showing a widespread dependency of groundwater-level time series on rain data.

Hierarchical cluster analysis
According to the dendrogram in Fig. 9c, the possible solutions are at k = 2, 3, 4 and 8, where k is the number of clusters. MWs have been grouped into eight clusters (C1-C8; Fig. 5 a Identified reversal points for groundwater (GW) level minimum and maximum conditions. Dates on the X-axis are provided according to format mm-yy. Spatial and temporal distribution of the identified reversal points for b GW minimum c GW maximum. NCI stands for 'not clearly identifiable' Fig. 6 Groundwater-level time series and reversal points for the selected monitoring wells (MWs). a-b MWs time series located in the north; c-d MWs time series located in the south; e-f MWs time series located downtown; g-h MWs time series located in the west. i-j MWs time series located in the east; k location of selected MWs time-series. To better compare the time-series, the X-axis scale has been set with intervals of six observations (i.e. 3 years); the Y-axis scale has been kept the same for all the time-series graphics (15 m) Fig. 9a,c), as this allowed a more detailed exploration and interpretation of the study area's environmental variability (Rinderer et al. 2019). Two clusters (C1 and C4) contain most of the MWs (23 out of 95 are contained in C1, while 34 out of 95 are included in C4); two clusters (C2 and C3) both contain 12 MWs, while the remaining four clusters (C5-C8) contain at most four MWs, thus describing more specific situations. The MWs of C1 are mainly located in the southern sector, while a few are close to the downtown area. These MWs are characterized by a relatively stable groundwater level, with hydraulic head variations of 2-3 m. This characteristic is due to the local geology for the southern MWs, showing groundwater heads close to the ground level, with a constrained possibility of oscillation; close to the downtown area, the nonconsumptive employment of groundwater exerted by GWHPs could contribute to maintaining groundwater levels as relatively stable. Cluster C2 includes 12 MWs located in a transition belt between the northern and the central sector. As illustrated in Fig. 9b, the main differences between C3 and C4 are the average groundwater level of the cluster and the range of hydraulic head variations (i.e. within ~5 m), while the cluster trends are rather similar. Despite this similarity, MWs in C3 are characterized by a steeper rising (around 8-10 m) phase in the period between Sep07 and Mar15 (groundwater minimum and maximum conditions respectively). C4 contains 34 MWs, distributed in the central portion of the domain: the average trend of this group is similar to C1, even if a less oscillating declining phase has been observed between Mar05 and Sep07; also for this cluster, stable groundwater levels could depend on the use of GWHPs. Three MWs located in the northern part of the domain are included in C5: a groundwater level peak for Sep08 was the main driver determining a medium groundwater level trend that differs from C3 and C8, whose MWs are similarly sited northward. C6 and C7 contain four and three MWs respectively: both these clusters group MWs located in the westernmost part of the study area; the average groundwater level of these clusters shows a high fluctuating trend, with no evident rising or declining phases. This may prove a rapid response of the water table to local conditions, as identified also from the results of the autocorrelation and cross-correlation functions. Four MWs, located in the northernmost area of the domain, are grouped in C8: the average trend of this cluster is comparable in its first part to the trend of C3, while the declining period from Mar15 is more oscillating than the linear declining phase of C3; thus, this second part of the trend is more comparable to C5 trend. The average hydraulic head variation is around 7-8 m.

Water table and groundwater volumes calculation
Groundwater potentiometric maps of the shallow aquifer for the identified groundwater level minimum and maximum conditions (Sep07, Mar15) and for Sep19 are shown  (Bloomfield et al. 2015). The 95% significance level in both cases is ±0.3578. For the legend of land use classes, refer to Fig. 4 in Fig. 10a-c and their raster differences are displayed in Fig. 10d,e. As observable in Fig. 10d, in the northern portion of the domain, characterized by a steeper rising trend (Fig. 4) and a wide unsaturated zone (Fig. 2), hydraulic heads (m asl) differed by more than 10 m between the minimum and maximum conditions, while in the NW this difference was limited to around 4 m. Moving south, the difference values oscillated between 3 and 5 m close to the downtown area; in the south and west, the differences in hydraulic heads were between zero and 2 m, confirming a reduced/not significant rising trend. In these areas, despite Sep07 being identified as the groundwater level minimum, in a few locations, the hydraulic head of Sep07 was higher than Mar15, which could be due to the irrigation activities generating a groundwater maximum at the end of the summer season. Moreover, as the unsaturated zone of the southern part of the shallow aquifer is thinner than in the northern part (Fig. 2), groundwater level oscillations are limited, and the groundwater minimum and maximum conditions are less evident. On the other hand, the eastern side of the domain shows different behaviour, with a global groundwater maximum higher than the global minimum. As regards the Mar15-Sep19 investigation (Fig. 10e), hydraulic head differences appeared limited with respect to the Mar15-Sep07 analysis (Fig. 10d). In the north, the hydraulic head differences did not exceed 6 m, with 3 m as the maximum difference observable near the downtown area, while a similarly limited oscillation was detected for the southern sector, if compared to the Mar15-Sep07 analysis. With respect to the western portion, hydraulic heads of Sep19 were generally higher than in Mar15, despite Sep19 coming at the end of a declining trend period. In terms of subsurface volumes impacted by these groundwater oscillations, the rising trend between Sep07 and Mar15 determined an occupation of 704.94 × 10 6 m 3 . In contrast, the declining phase that took place from Mar15 onward determined a reduction of 289.62 × 10 6 m 3 of groundwater in the shallow aquifer. As for the mobile water included in these volumes, eight classes of effective porosity have been introduced for the calculation, deriving their values from the hydraulic conductivity values of the lithologies in these subsurface portions. A detailed description of these classes (n e ) has been provided in Table 1.
Effective porosity values ranging from 0.18-0.19 have been attributed to fine grain-size materials (hydraulic conductivity (K) ranging from 38 to 70 m/s), while higher values have been assigned respectively to medium-grain Based on this classification, for the Sep07-Mar15 period there has been an estimated increase of 149.53 × 10 6 m 3 of mobile water, while for the Mar15-Sep19 period, a decrease of 62.69 × 10 6 m 3 .

Discussion
As planned by the local authorities (Milan Metropolitan City 2019), in the framework of a general urban transformation, urban underground development is occurring for the city of Milan. This implies an integrated management scheme, taking groundwater aspects into account. Thus, a data-driven methodology was adopted in this study to identify groundwater management areas (MAs) and elaborate management proposals for the underground infrastructure. In particular, greater attention has been paid to public car parks. To do so, a set of statistical and geospatial techniques, which are scientifically appropriate, has been adopted to assess spatial-temporal patterns of groundwater levels (Chaudhuri and Ale 2014).

Improved conceptual model
The obtained results showed a general spatial agreement, identifying different hydrogeologic features among various domain sectors. The northern sector showed the highest trend magnitude (Fig. 4): in this area, the greater unsaturated thickness (Fig. 2) allows the water table to oscillate widely. This hydrogeologic condition facilitated the groundwater rebound that started at the beginning of the 1990s, triggered by both the industrial decommissioning and a gradual decrease in the resident population (Bonomi 1999;De Caro et al. 2020). In the same sector, the autocorrelation coefficient turned out to be significant at up to six lags, implying a long-term memory of the system and a significant dependency from previous measures (Dinpashoh et al. 2014;Neto et al. 2016). The downtown sector showed a lower trend magnitude than the northern sector. This lower trend can depend on a possible anthropogenic influence associated with extensive use of GWHPs and a significant presence of underground infrastructure controlling groundwater levels. These areas showed a delayed response to precipitation events, with a significant correlation at lag -1 (response time: 6 months): this could be due to (1) the greatest depth of the water table (Fig. 2), mainly in the north, which determines a delayed groundwater response, and (2) a reduced infiltration rate dependent on an urbanized land use, visible for both sectors (Fig. 4). These factors contributed also to a delayed identification of reversal points, as occurred for the groundwater minimum: in the north this condition was in part identified in Mar08 (Fig. 5b), and not in Sep07 as for the majority of the time series (Fig. 5a).
In the south, depending on its hydraulic gradient, the water table is close to the ground level (Fig. 2), which limits the rising trend. Moreover, as for the western sector, a decrease in grain size from coarse (i.e. gravel and sand) to fine sediments (i.e. silt and clay), moving from N to S, determines a change in the aquifer transmissivity, which leads to the formation of the lowland springs, forcing groundwater both close to the ground level and to outflow, thus constraining water-table oscillations (De Luca et al. 2014). Furthermore, in the west, the presence of clay lenses determines the formation of local perched aquifers, with a medium depth ranging from 6-8 m below the ground level ). Western time series are also marked by the presence of seasonal fluctuations (Fig. 9b), with a groundwater level minimum in March and a maximum in September. Land use could influence this cyclical condition. This area is mainly agricultural (Fig. 4), with maize and rice as primary cultures (Bernini et al. 2019; Pulighe and Lupia 2019): as the irrigation period goes from April to September, a groundwater level maximum is expected at the end of this season, with a water-table response generally delayed by less than 1 month compared to the beginning of the flooding activities (Lasagna et al. 2020). All these factors could explain the reduced significant autocorrelation in the south and the almost complete absence of autocorrelation in the west, highlighting a lack of system memory in these areas and a rapid response to external drivers. The significant cross-correlation with  rain data at lag 0, providing insights into a rapid aquifer recharge by rainfall (Lasagna et al. 2020), is in line with these outcomes. This aspect also appeared from reversal points analysis, as in these areas the identified groundwater level minimum and maximum condition partly occurred in advance from Sep07 (groundwater minimum) and Mar15 (groundwater maximum) (Fig. 5).
The widespread correlation at lag -2 (response time: 1 year) could be attributed, in particular for these peripheral sectors of the domain, to a sort of periodicity of groundwater time-series, mainly attributable to the agricultural activities that led to a cyclical hydraulic maximum at the end of the summer season.
In contrast, different hydrogeologic settings could explain the discrepancy in correlation results between western and eastern sectors (Fig. 7a), as land use distribution is mostly similar. In the east, the paleochannel of the Lambro River ( Fig. 1), with coarse materials and a higher transmissivity, determines a difference from the geological system of the lowland springs and the local perched aquifers previously described for the western sector, influencing correlation results.
These findings highlighted that the degree of autocorrelation influences the slope magnitude of groundwater level trends. In addition, the thickness of the unsaturated zone was revealed to be a contributory factor to MWs autocorrelation (Bloomfield and Marchant 2013). As regards cross-correlation with rain data, results pinpointed that lag times between the groundwater level responses to precipitation events are related to the lithological conditions and the thickness of the unsaturated zone (Kotchoni et al. 2019). Moreover, land use distribution could be considered as an additional factor influencing the infiltration rate (Scanlon et al. 2005).
The final grouping of groundwater time series in eight clusters is consistent with the results previously described. Northern time series are grouped into C3, C5 and C8: these time series both show the highest trend magnitude and the highest degree of autocorrelation; moreover, their correlation with rain data is uniformly delayed (lag -1). In addition, their water table rising, of between 8 and 10 m, was the greatest detected. This trend proves the anthropogenic impact (Bloomfield et al. 2015), representing the ending phase of the groundwater rebound, triggered by the industrial decommissioning that started at the beginning of the 1990s (Bonomi 1999). The same explanation provided for the northern clusters could be provided even for C2, grouping MWs located at the transition belt between the northern clusters and C4. This latter cluster groups all the groundwater time series located in the downtown belt, showing both a dampened magnitude of trend and degree of autocorrelation compared to the northern MWs. C1 includes all the southern MWs, with a part of them stretching up to the downtown area; this spatial distribution could depend both on a structural and anthropogenic factor. The structural factor is due to the constrained possibility of oscillations of the water table in this part of domain; the anthropogenic influence could be attributed to the groundwater exploitation for heat exchange; GWHPs active in the downtown area, mainly inducing a thermal stress, could be considered as nonconsumptive (García-Gil et al. 2016), as by using a system of catch and return wells, the water withdrawn is usually returned to the shallow aquifer. C6 and C7 group all the MWs located in the western area; as discussed previously, the fluctuating but constrained groundwater level trend of these MWs could be related both to the presence of lowland springs and local perched aquifers and to the land use, with irrigation activities conducted during the summer season. Thus, also cluster analysis emphasized that, even within the same aquifer, and in a rather limited domain, significant differences in groundwater behaviour could be identified.
The outcomes of the data-driven approach have also been confirmed by the spatial analysis, applied to reconstruct the potentiometric maps of different groundwater conditions and to calculate the subsurface volumes influenced by these groundwater oscillations. In the north, where the highest trend magnitude (Fig. 4) and degree of correlation were detected (Fig. 7), the major groundwater hydraulic head variations have also been identified using GIS. Subsequently, a progressive dampening in water-table oscillations has been registered moving both southward and westward, where a reduction in trend magnitude (Fig. 4) and degree of correlation (Fig. 7) was observed, while a more rapid correlation with rainfall data was identified (Fig. 8a). Moreover, the content of mobile water included in these subsurface volumes has been estimated. This information could be useful for stakeholders, as in groundwater rising phases part of this amount of water could represent a flooding threat to underground structures; for this reason, this water could be locally withdrawn to secure the infrastructure at risk, exploiting this resource for other needs (i.e. street washing, recreational uses).

Considerations on the data-driven approach
This agreement among the obtained results allowed this study to focus on the effectiveness of a data-driven approach. Groundwater-level time series have been analysed, applying different statistical techniques-each one was employed to investigate a specific hydrogeologic aspect; however, in the end, the final outcomes obtained always showed spatially coherent regions of groundwater behaviour. Thus, timeseries analysis provided a deep insight into the main hydrogeologic features of the shallow groundwater system, highlighting fundamental properties that affect groundwater level dynamics (Bloomfield et al. 2015); applying them to support groundwater system conceptualization could guarantee both time-and cost-effective assistance to groundwater models (Obergfell et al. 2013), providing useful complementary information. Besides, even if for many groundwater questions an answer could be obtained only by structuring a spatial groundwater model (Anderson et al. 2015), timeseries analysis can still have a relevant role (Bakker and Schaars 2019)-for example, in data preprocessing and target definition.
In particular, through the approach developed in R to identify groundwater conditions, the highest number of reversal points was identified for the western and southern time series (Table S2 of the ESM), which proved the absence of a clear significant trend, evidencing a more fluctuating behaviour if compared to the most visible rising and declining phases detectable in the northern and central sectors. This outcome is consistent with the hydrogeologic framework described through the Mann-Kendall test and Sen's slope estimator. Thus, the adopted approach turned out to be reliable in identifying the main groundwater conditions typical of the considered shallow aquifer, also identifying local situations.
As regards cross-correlation, combining groundwaterlevel time series and precipitation data provided a better understanding of at least one natural forcing on groundwater trends (Kumar et al. 2009). The adoption of just one meteorological station was a direct consequence of a lack of data and fragmentation, which is still a limiting factor as regards the potential use of Open Data information (Janssen et al. 2012). However, as the study area is a limited domain, located in the middle of an alluvial plain, a rather homogeneous rainfall trend over the area is presumable, considering that only one weather station was believed to have obtained reliable results. Depending on data availability, future perspectives could consider using a higher number of meteorological stations to strengthen this analysis.
Finally, grouping groundwater time series through cluster analysis techniques, focusing on groundwater levels, constitutes an emerging research field (Naranjo-Fernández et al. 2020). In particular, its application represents an element of novelty for the study area: in fact, statistical techniques in this domain have been mostly applied for groundwater chemicals (Alberti et al. 2016;Azzellino et al. 2019;Pollicino et al. 2019) and thermal issues (Previati and Crosta 2021b), but not for quantitative analysis. Additionally, the adopted approach allowed for the accurate pinpointing of the potential features influencing groundwater levels (i.e. geology, precipitation, withdrawals, land use), thus evidencing more defined local groundwater conditions. Moreover, their application in an urban planning framework is supposed to guarantee better policy outcomes. To the best of the authors' knowledge, this is the first time that data-driven analysis of groundwater levels have been adopted in supporting underground infrastructure management, and this approach could help stakeholders in effectively designing the already planned urban underground development that will take place in the next years for the city of Milan, considering also groundwater aspects.

Decision management
In view of the previous results, which showed a similar spatial zoning, four management areas (MAs) have been identified (Fig. 11). Each area showed distinct local hydrogeologic behaviour, thus all need different management policies regarding urban groundwater management and urban planning. In order to elaborate effective guidelines for the management and construction of new public car parks, an already defined spatial methodology has been previously adopted to evaluate how their interference with the water table in terms of flooding (not submerged, possibly critical, submerged) has varied over time. The methodology has been fine-tuned with respect to a previous work (Sartirana et al. 2020), introducing a threshold value that allows one to identify potentially submerged infrastructures (i.e. possibly critical structures). Given the total number of 126 car parks, 35 car parks are located in MA1; 54 and 30 car parks are respectively sited in MA2 and MA3, while only 7 are in MA4. The results of the spatial relationship between public car parks and the water table have been summed up in Table 2 and are visible in Fig. 11. The MA for each public car park is found in Table S1 of the ESM.
Thus, car park submersion mainly occurred in MA2, followed by MA3. This depends also on the depth of these structures-in the downtown area (MA2), car parks can reach up to six floors (i.e. 20 m of depth); in contrast, in the northern sector (MA1) car parks have a medium depth of two-three floors (i.e. 8-11 m of depth), while one floor (i.e. 5 m) is typical in the western sector (MA4). Finally, car parks located in MA3 range from one to two floors in the southernmost part, and up to five floors in the northernmost part, at the border with MA2 (Sartirana et al. 2020). Thus, MA3 also presented recurring impacted car parks, due to both the number of floors and the proximity of the water table to the ground level (Fig. 2). Finally, MA1 did not show any impact for Sep07 (Fig. 11a) and Sep19 (Fig. 11c), while some parks results were noted as being possibly critical with respect to flooding in Mar15 (Fig. 11b), as a consequence of the steep water-table rise occurring in the northern sector. Even if the unsaturated thickness is high (Fig. 2), adopting proper management measures would be recommended. As for MA4, the worst condition has been detected in Sep19, because of the seasonal water-table fluctuations in this sector (Fig. 9b), determining a local groundwater level maximum for the end of the summer. As a consequence of these results, a list of suggested guidelines regarding the construction of new public car parks has been elaborated. These guidelines, together with the hydrogeologic features and the possibility of submersion for each MA, are summarized in Table 3.
In the past, some underground structures have been designed without waterproofing, neglecting the possibility of a groundwater rebound. Consequently, in the last decades, there have been flooding episodes in some public car parks (Sartirana et al. 2020). For this reason, a general guideline on the adoption of waterproofing techniques for new underground infrastructure has been proposed, acting as a reminder for the administrators. Based on local hydrogeologic conditions, a suggested number of floors has been provided for new public car parks for MA1, MA3, and MA4. A maximum of two floors has been indicated only for MA1, aiming to avoid possible interactions with the water table. A similar proposal has not been developed for MA2, as the majority of the economic activities are located in the city centre. This supposes a high demand for new infrastructure due to socio-economic needs (Bobylev 2009), and considers that the already-in-place structures are developed as multistorey, specific guidelines have been elaborated considering an interaction with the water table.
In particular, pumping systems, including the adoption of GWHPs, could be favoured to control the water-table rise, regulating the withdrawals based on regular monitoring of the interaction between groundwater and the underground infrastructure. Moreover, in order to assess the impact of the underground infrastructure on groundwater (Attard et al. 2017), local 3D numerical models should be developed. In this way, a sustainable and integrated management scheme will be developed, balancing the needs of both a secure groundwater resource and underground infrastructure.

Conclusions and future perspectives
This work aimed at proposing a data-driven approach to support decision-based management for urban planning of new underground infrastructure, focusing on public car parks. To facilitate this, a time-series analysis of groundwater-level time series of the shallow aquifer beneath the city of Milan has been conducted, as interaction with groundwater has been observed in the last decades, resulting in the adoption of targeted management strategies in the urban planning framework becoming necessary. In particular, this study has allowed this paper to: 1. Define a workflow of different statistical and geospatial techniques that allows one to evaluate the spatio-temporal patterns of the groundwater levels for the studied shallow aquifer, proving the effectiveness of a datadriven approach. 2. Identify portions of the study area showing different groundwater conditions, mainly attributable to local hydrogeologic settings. Indeed, even in spatially limited domains, a thorough monitoring network could guarantee a broad insight into groundwater processes. Thus,  different groundwater management strategies should be adopted according to specific hydrogeologic situations. 3. Provide decision-makers a (site-specific) methodology to adequately design and handle the planned underground development of the city of Milan, also considering groundwater aspects. In particular, as a main outcome of this study, four areas of groundwater/underground-infrastructure integrated management options have been identified, whereby different guidelines for managing and designing new underground public car parks are recommended for adoption.
In order to properly design urban groundwater management strategy, this approach could be complementary to groundwater numerical modelling. Indeed, analysing groundwater-level time series could allow one to improve the definition of the conceptual model, and to identify local hydrogeologic conditions to be further studied through the development of a specific numerical model. Future perspectives will consider the realization of different groundwater scenarios, simulating possible climate change effects, to evaluate the dynamics of the interferences between groundwater and the underground structures, thus updating when required the suggested guidelines. The adopted methodology could be extended worldwide to other urban areas that have a well-structured groundwater monitoring network, are impacted by an existing interaction between groundwater and the underground infrastructure, and may face future underground development.

Acknowledgements
The authors are grateful to Metropolitana Milanese S.p.A for providing piezometric data for groundwater-level timeseries reconstructions. The authors would like to thank the associate editor, Hamid Kardam Moghaddam and another anonymous reviewer for their comments, which helped to improve this article.
Funding Open access funding provided by Università degli Studi di Milano -Bicocca within the CRUI-CARE Agreement. This research did not receive external funding.

Conflict of interests The authors declare no conflict of interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.