*Article* **Quantifying Groundwater Infiltrations into Subway Lines and Underground Car Parks Using MODFLOW-USG**

**Davide Sartirana \*, Chiara Zanotti , Marco Rotiroti , Mattia De Amicis, Mariachiara Caschetto, Agnese Redaelli, Letizia Fumagalli and Tullia Bonomi**

> Department of Earth and Environmental Sciences, University of Milano-Bicocca, Piazza della Scienza 1, 20126 Milan, Italy

**\*** Correspondence: d.sartirana1@campus.unimib.it

**Abstract:** Urbanization is a worldwide process that recently has culminated in wider use of the subsurface, determining a significant interaction between groundwater and underground infrastructures. This can result in infiltrations, corrosion, and stability issues for the subsurface elements. Numerical models are the most applied tools to manage these situations. Using MODFLOW-USG and combining the use of Wall (HFB) and DRN packages, this study aimed at simulating underground infrastructures (i.e., subway lines and public car parks) and quantifying their infiltrations. This issue has been deeply investigated to evaluate water inrush during tunnel construction, but problems also occur with regard to the operation of tunnels. The methodology has involved developing a steady-state groundwater flow model, calibrated against a maximum groundwater condition, for the western portion of Milan city (Northern Italy, Lombardy Region). Overall findings pointed out that the most impacted areas are sections of subway tunnels already identified as submerged. This spatial coherence with historical information could act both as validation of the model and a step forward, as infiltrations resulting from an interaction with the water table were quantified. The methodology allowed for the improvement of the urban conceptual model and could support the stakeholders in adopting proper measures to manage the interactions between groundwater and the underground infrastructures.

**Keywords:** urban hydrogeology; rising groundwater levels; shallow aquifer; 3D geodatabase; horizontal flow barrier; Milan; Italy

#### **1. Introduction**

Urban hydrogeology is a specific branch of research [1,2] that has been constantly developed in recent years as a consequence of rapid urbanization phenomena that have been witnessed in most parts of the world [3]. Considering that 70% of the world population is expected to live in urban areas by 2050 [4], urbanization can be defined as a world-wide process [5]. Thus, it is reasonable to think that in the next few years a huge effort will be allocated to research into urban hydrogeology [6].

Overexploitation and deterioration of urban water resources act as the main consequences of this rapid urbanization [7]. To put a brake on urban sprawl, a vertical urban development has occurred, determining an augmented use of urban underground [8–12]. However, the construction of ever-deeper structures [13] can impact groundwater (GW) with regards to flow, quality, and thermal issues [5,14,15].

With respect to GW flow, different cities around the world have observed rising water table levels, as a consequence of the deindustrialization process, that have generated some interference between GW and underground infrastructures (UIs) such as basements, car parks, and subway lines [16–23]. Numerical GW flow modeling was widely adopted as the main tool to evaluate the barrier effect of UIs to flow patterns, GW budget [14], and the possible side effects on the underground elements (i.e., corrosion and stability issues).

**Citation:** Sartirana, D.; Zanotti, C.; Rotiroti, M.; De Amicis, M.; Caschetto, M.; Redaelli, A.; Fumagalli, L.; Bonomi, T. Quantifying Groundwater Infiltrations into Subway Lines and Underground Car Parks Using MODFLOW-USG. *Water* **2022**, *14*, 4130. https://doi.org/10.3390/ w14244130

Academic Editor: Craig Allan

Received: 11 November 2022 Accepted: 15 December 2022 Published: 19 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Concerning engineering issues, GW inflow into tunnels has been predicted in urban areas by adopting analytical solutions [24], synthetic modeling [25], and steady-state numerical modeling on real cases [26,27] to properly design the tunnel drainage system during the construction phase. In fact, water inrush is a challenging issue to face, causing negative impacts on tunnel stability, generating subsidence damage [25,28], heavy financial losses, and losing construction time [26,29].

At the same time, the problem of damages in operating tunnels, as water seepage or lining cracking, requires consideration [30–33]. Despite the lower water amounts penetrating inside the UIs over a long period, GW could determine severe issues, such as temporary unusability, which require waterproofing works and lead to economic losses. Thus, quantifying infiltrations could help to assess proper mitigation strategies [34], supporting the stakeholders in the complex task of urban GW management. To do so, among the different approaches applied in the literature, groundwater infiltrations into subsurface elements have been evaluated by modeling the underground infrastructures by means of the DRN package [34–36]. Recently, a single model layer was developed by Golian et al. [37] to restore groundwater levels after tunneling. In this work, an unsealed and a sealed underground tunnel were modeled using RIV and HFB packages, respectively. The latter has been applied in various fields of groundwater modeling: from coastal areas to model slurry walls containing seawater intrusion [38,39], to geophysical modeling to simulate faults [40,41], to urban contexts in industrial sites [42], or to evaluate the impact of underground infrastructures on groundwater levels [43,44].

The existence of 3D geodatabases, gathering information on underground structures [45,46] and frequently scattered over many institutions and stakeholders [1,47–49], could support the adoption of these packages to properly model UIs. In this way, it should be possible to precisely define their relationship with the water levels, thus improving the urban conceptual model.

Based on these assumptions, the aim of this study has been to quantify GW infiltrations into different categories of UIs (i.e., subway lines and underground car parks), considering different UI conditions (i.e., intact, saturated, and leaky walls). The methodology that has been applied involves developing a local 3D GW numerical flow model for the western area of Milan metropolitan city (Lombardy, Northern Italy). Through this model, the most critical portions of the subsurface network suffering from GW infiltrations have been evaluated. Interactions with the water table and possible infiltrations in subway line M4 (to be inaugurated in 2023) and two public car parks that are currently under construction were also analyzed.

By means of this model, the stakeholders would be able to design management solutions to secure the infrastructures from being flooded in the future. The model has been realized as steady-state with MODFLOW-USG [50] and calibrated using a trial and error approach against a GW maximum condition that was defined in a previous work as documented by Sartirana et al. [51]. HFB and DRN packages have been coupled to model the UIs, reproducing their geometries and volumes through the adoption of grid refinement, contributing to the quantification of GW infiltrations into subsurface elements. In particular, the top and the bottom of the UIs were modeled through the HFB package; to the best of the authors' knowledge, this application of the HFB package could represent an improvement in modeling the UIs. In fact, the relation between GW and the UIs along the vertical sides of a model cell could be thus considered. Moreover, as for Milan city, this is the first time that car parks have been considered in a 3D GW numerical flow model, while being studied for the adoption of GW-level time-series clustering to suggest targeted guidelines for the construction of new underground public car parks [51].

The methodology presented here could be implemented for other urban realities, serving as a way of managing a documented interaction between GW and the UIs that may lead to a planned subsurface infrastructure development with possibly great potential for an integrated management strategy.

#### **2. Urban Conceptual Model of the Study Area**

The study area covers 100 km2 inside the Milan metropolitan area (Figure 1). Human activities have always characterized this area, especially through industrial and agricultural activities that are still conducted in the western and southern areas of Milan [52,53]. The city hosts 1.4 million inhabitants [54] and is currently undergoing an important urban transformation [55]. It is located in the middle of the Po Valley, whose hydrogeologic structure has been deeply examined both in the past [56] and recently [57]. Three main hydro structures were identified: a shallow hydro structure (ISS), an intermediate (ISI), and a deep (ISP) hydro structure. Within the model domain, an ISS has a medium thickness of 40 m with a bottom surface ranging from 100 m above sea level (a.s.l.) (to the north) to about 60 m a.s.l. (to the south). It hosts a shallow aquifer (Figure 2) (i.e., Aquifer Group A1, Regione Lombardia and ENI Divisione AGIP 2002 [56]), where all the underground infrastructures are located. This aquifer is not exploited for drinking needs. Sands and gravels mainly characterize this hydro structure. The same lithologies, but with an increasing presence of silty and clayey horizons, constitute the ISI, that mostly corresponds to Aquifer Groups A2 and B of Regione Lombardia and ENI Divisione AGIP 2002 [56]. An ISP, having a more uncertain lithological composition, was not modeled within this study.

**Figure 1.** (**a**) Geographical setting of the study area; (**b**) main hydrogeologic features (lowland springs) Color coding for the subway lines respects the color coding used by the subway managing company. Public car parks have been represented as triangles to differentiate them from wells. (Image readapted from Sartirana et al. [51]).

**Figure 2.** Hydrogeologic schematic cross sections AA' (N–S) of the study area, showing the location of some UIs and their relationship with the groundwater condition of Mar 2015 [51]. For their location on map, please refer to Figure 1.

Industrial needs triggered an extensive groundwater withdrawal since the early 1960s. Consequently, the water table reached its maximum depth of more than 30 m in the northern part of the city around 1975, thus determining the minimum GW levels due to significant water exploitation [58,59]. During the same time frame, some UIs (car parks, subway lines M1 and M2) were built, sometimes without proper lining methods, without consideration for a possible future GW level rise. Subsequently, since the beginning of the 1990s, the decommissioning of many industrial sites, mainly located in the northern sector of the city, generated a rise in GW level, determining flooding episodes for these oldest and shallowest subway lines and for some underground car parks built starting from the middle of the 1980s [60,61]. Consequently, the most recent and deepest subway lines (M3, M4 to be inaugurated in 2023, and M5) have been designed with lining systems. As for underground car parks, 126 public car parks are now listed in the city [51]: 65 out of 126 are located in the model domain. The construction of two new underground car parks (Figure 1b) is currently taking place close to the Gelsomini and Frattini stations of subway line M4. These car parks are named Brasilia (placed just northward of the stations) and Scalabrini (to the south of the stations), respectively; both have been designed to be two floors deep (i.e., 8 m depth as calculated by Sartirana et al, 2020).

The water table rise occurred differently among different areas of the town, with a maximum rise of about 10–15 m in the north, and a more dampened effect in the other sectors [51]. Particularly, a low significant rising trend was evidenced in the west and south, respectively, due to local geological conditions and the hydraulic gradient that constrains the water table close to the ground level, thus reducing the water table oscillations.

In the downtown area, an increasing presence of open-loop groundwater heat pumps (GWHPs) for geothermal needs (Figure 1b), together with the presence of a huge number of UIs, could induce an anthropogenic control on water table rising; due to extraction and injection wells systems, the water withdrawn is usually returned to the shallow aquifer, thus determining a non-consumptive use of the resource [62]. These systems sometimes discharge exploited water to surface water bodies to control the GW rise.

Public-supply well fields withdraw water used for drinking needs, and have screens to tap the semi-confined and confined aquifer units. A total of 261 wells, belonging to 13 well fields, are located inside the considered domain.

The construction of new underground car parks takes place in the framework of the adoption of the Plan of Government for the Territory (PGT) [63], that regulates further subsurface occupation as a measure against excessive soil consumption. In this context, numerical modeling, possibly combined with the application of other techniques aimed at better understanding the urban conceptual model [51], could represent a valid tool to coordinate urban underground development, thus supporting stakeholders in their decision-making process.

#### **3. Materials and Methods**

The numerical model was built considering an already-existing urban conceptual model [51], integrating its contents, when possible, with Open Data information [64]. The core of the methodology was the modeling of the UIs (see Underground Infrastructures Modeling) to evaluate GW infiltrations. Different scenarios of conductance were realized to quantify infiltrations simulating different wall conditions; the results have been examined in order to discuss possible strategies to manage GW/UI interactions.

#### *3.1. Numerical Model*

A steady-state numerical flow model was developed using MODFLOW-USG [50], and Groundwater Vistas 8 [65] was used as the graphical user interface.

#### 3.1.1. Grid Discretization

The model grid (Figure 3) was composed of 1,668,348 cells and was horizontally structured by applying a quadtree refinement: cell dimension ranges were from 100 m in the peripheral areas, up to 12.5 m around subway lines and public car parks (i.e., fourth level of refinement); in proximity to public car parks currently under construction, a fifth quadtree level of refinement was applied (i.e., 6.25 m) (Figure 3a). The grid was rotated by 35◦ from the offset (X = 1,509,407, Y = 5,026,235, Monte Mario Italy 1; ESPG: 3003) to be perpendicular to the general NW–SE groundwater flow direction of the domain [58,66].

The vertical discretization (Figure 3c) consisted of 18 layers. The first 8 layers, with an average thickness of three meters, included all the UIs lying in the shallow aquifer (layers 1–10, ISS/Aquifer Group A1); layers 11 to 14 had a medium thickness of seven meters to model the first portion of the ISI (Aquifer Group A2). Layers 15 to 17 (with a medium thickness of 6 m) were adopted to represent the aquitard (AQ), while the last layer, with a medium thickness of 20 m, aimed at modeling the final portion of the ISS (Aquifer Group B).

#### 3.1.2. Boundary Conditions

Boundary conditions (Figure 4), used to outline the hydrogeologic system, were represented through Neumann and Cauchy conditions:


their well discharge was mostly unknown, a discharge value of -432 m3/d was initially attributed to these wells.

• Recharge (RCH): 5 zones, based on land use, were identified from the geographic database Dusaf 6.0 [67]; their values were calculated as the contribution of precipitations, irrigations, and runoff. The initial values for each zone were calculated starting from the precipitation data of Paderno Dugnano rain gauge (located just northward of the city of Milan), monitored by the regional environmental protection agency [68]. Precipitations amounted to 1496.2 mm/yr for the twelve months before Mar15, the period chosen for model calibration. Absence of infiltration was considered for urban areas and for surface water elements (i.e., quarries), while 20% of infiltration was attributed to the other recharge areas; moreover, an additional contribution from recharge infiltration was attributed to irrigational areas.

**Figure 3.** (**a**) Grid horizontal discretization; the red rectangle points to the sector area represented in (**b**); (**b**) example of quadtree refinement close to Lotto exchange station (see Figure 2); (**c**) grid vertical discretization. Please note that for (**b**), the same color coding of Figure 1b has been used for subway lines.

#### Underground Infrastructures Modeling

The underground railway (Figure 1b) occupies only a small portion of the northwestern area; thus, it was not considered within the study. All the UIs (Figure 4) were conceptualized and modeled by coupling the Wall (HFB) [69] and the DRN [70] packages. The capabilities of both packages were combined to properly simulate and evaluate the exchange between the UIs and the surrounding aquifer. HFB offers the ability to isolate individual components to consider how water is passed between an engineered element such as a subway line and the aquifer. On the other hand, the DRN package enables the modeler to assign a head inside the engineered structures. In this case, the DRN package was adopted to simulate a fictitious water-collection system within the UIs. Combining the capabilities of these two packages is the core of the proposed methodology. To support and validate the adopted methodology, the model domain was discretized into 187 zones: the aquifer of interest, and all the UIs' elements. Thus, water exchanged between neighboring zones, based on the MODFLOW solution [71,72], was quantified.

To properly model the depth of the UIs, information regarding the UIs' bottom and the diameter of the subway tunnels has been obtained from an already-existing 3D GDB of the subsurface elements for the study area [46]. Subsequently, the following rule was adopted as the main constraint to model the UIs: if an UI occupied a layer of more than 50% of its thickness, the UI was then represented inside that layer; otherwise, if this constraint was not respected, the UI was then modeled in the overlying layer.

The wall usually goes along any of the four horizontal sides of each cell, but in MODFLOW-2005 there is no option to specify a vertical barrier. Notwithstanding, the adoption of MODFLOW-USG allowed for wider flexibility in using the HFB package, as the barrier could be aligned along any face of the unstructured grid [50]; thus, HFB cells could also be placed at the intersection between two nodes sharing the same X and Y coordinates, in contiguous layers. This enabled the reproduction not only of the lateral sides, but also the top and the bottom of all the subsurface elements. To do so, the initial information about the lateral sides of the UIs was integrated by "manually" compiling the HFB package, adding the position of the top/bottom of the UIs.

The drainage network was placed inside the UI and positioned at the bottom layer of each section of the UI. The drain head (i.e., drain elevation) was fixed as equal to the bottom elevations of the UI. In this way, the possible groundwater inflow into the UIs could be drained, quantifying the amount of water to be withdrawn to dry the infrastructure. The conceptual model of the adopted approach to simulate the underground infrastructures network is represented in Figure 5.

**Figure 4.** Model boundary conditions. GHBs' distance from the model area has been indicated. Please note that color coding of the infrastructural elements (subway lines and underground public car parks) refers to the HFB package color in Groundwater Vistas 8. Public car parks have been represented as triangles to differentiate them from wells.

**Figure 5.** (**a**) Traditional application scheme for HFBs cell; insertion mask taken from Groundwater Vistas 8; (**b**) conceptual model of the adopted approach to model all the UIs.

At the exchange stations, lines were positioned at their real depth, thus properly separating the deepest and more recent lines (i.e., M3, M4, and M5) from the shallowest and older ones (M1 and M2).

Conductance, which can be defined as the ratio between hydraulic conductivity (K) and the wall thickness (d), is the single parameter that controls the ability of the wall to transmit water. The absence/presence of lining systems was represented through different conductance values. With regard to the wall thickness, a value of 1 m was considered representative of all the modeled UIs.

The drainage system was assumed to provide no resistance to GW flow, imposing a value of conductance higher than the wall conductance and the aquifer conductivity [25,35].

#### 3.1.3. Further Modeling Aspects

The hydraulic conductivity parametrization was readapted from a previous project on the study area developed within the same research group [73], where the lithological information, stored within the Tangram database [74] in the form of stratigraphic data and pumping tests, was numerically coded and interpolated into GOCAD software using the kriging method [75]. Initial values to the continuous distribution of hydraulic conductivity were assigned from Tangram reference tables. A refined investigation was conducted which analyzed 3 cross-sections built along public-supply well fields from Airoldi and Casati [76], to infer the spatial distribution of fine materials (i.e., clay lenses).

With regard to calibration, sensitivity analysis using different multiplying factors (from 0.5 to 1.5) and a "trial and error" method were adopted to calibrate the steady-state model, focusing on GHB values and conductance, aquifer recharge (3 out of 5 zones), hydraulic conductivity, and well discharge. A total of 30 head targets, representing field water table measurements, were considered, showing an uneven distribution over the entire domain, with a limited amount of information for the western sector. The calibration process was conducted against the maximum groundwater condition of Mar15, the highest in the last 30 years [51], evaluating the goodness of the obtained results and analyzing the model statistics (i.e., residual sum of squares, scaled RMSE). In this way, the most critical situation for the UIs should be considered; this is also recommended for UIs currently under construction.

#### *3.2. Decision Management Support*

Different scenarios were analyzed to quantify GW infiltrations into UIs. Further engineering aspects, such as possible subsidence issues due to the drainage effect, or potential negative effects determined by buoyancy as a result of the aquifer pressure (i.e., uplift risks), were not considered within the aims of the project.

The conductance value for waterproofed subway lines (Table 1) was defined from the literature references [43,77,78]. Different conductance values (S1–S3) were tested for subway lines M1 and M2 due to a higher uncertainty; considering the absence of lining systems, the conductance was modified simulating possible deteriorations due to a prolonged interaction with the water table over time. In fact, infiltrations may be regarded as a gradual process, ranging from an unsaturated to a saturated flow induced by groundwater flow [79]. Since for public car parks' conductance no information was available, it was decided to attribute the lowest conductance value to all car parks.


**Table 1.** Conductance value for all the considered scenarios.

The most impacted locations of S1–S3 were then analyzed, locally increasing the conductance value of the HFB cells to simulate possible wall fractures. A focus was provided only for subway lines M1, M2, and car parks, as historically they have shown the most revealed interference. To reproduce fractures, wall conductance was only modified close to the infiltration area, increasing the initial value of four order magnitudes, as considered in studies on fractured rocks [80]. The change in the conductance was applied to the minimum model dimension (i.e., one cell). In this way, it was possible to compare the amounts of infiltration of intact and leaky walls.

The identified infiltrations were then analyzed to discuss some management proposals with regard to the design of dewatering systems in the most critical locations of the subsurface network, and also proposing the implementation of monitoring systems to manage possible infiltration issues in advance.

#### **4. Results**

#### *4.1. Model Calibration and Statistics*

The final values of GHBs were 127 m a.s.l. for the northern GHB and 102.2 m a.s.l. in the south, while the western and the eastern boundaries varied from 126 to 103 m a.s.l. and from 124 to 103 m a.s.l. from north to south, respectively. Calibrated values of hydraulic conductivity ranged from 235 to 1.15 × <sup>10</sup>−<sup>3</sup> m/d, as visible in Figure 6.

Final recharge values, and their spatial distribution, are represented in Figure 7. Finally, well discharge was reduced by 25% for the GWHPs and private wells, while for well fields, the reduction, when applied, ranged from 25% up to 50% (for the southernmost well field) of the initial value.

With regard to the calibration, the calibrated model generally provided good statistics (Table 2, Figures 8 and 9) for most of the 30 head targets considered. The most critical targets were located in the western and southernmost portions of the domain, quite far from the subsurface network that was the main focus of the study. Although these values could represent some modeling issues for some local areas of the domain, the scaled RMSE (4.6%) respects the international criteria that indicate the goodness of a solution in a scaled RMSE to be less than 8% [81,82].

**Figure 6.** Hydraulic conductivity values for layers (**a**)1(**b**)4(**c**) 11 (**d**) 14. Please note that subway line tracks are plotted inside all layers to provide refence points, since the grid is not rotated in these images.

**Figure 7.** Areal distribution of the 5 recharge zones; final recharge values are provided in legend.


**Table 2.** Model statistics for the considered head targets. Statistics refer to S1.

**Figure 9.** GW potentiometric map of the study area.

The potentiometric map for the shallow aquifer is represented in Figure 9. From the visualization of the head targets, it is visible that the water table map is generally well represented close to the subsurface network, thus allowing for proper assessment of GW/UI interactions and the consequent infiltrations. In the eastern part of the models, located close to Milan's downtown area, the contour lines' behavior is influenced by the pumping effect of both public well fields and GWHPs (Figure 4).

Model mass balance (Table 3) evidences the importance of well discharge inside the domain, both for the outflows and the inflows; the latter are exclusively due to the injection wells of GWHP systems. The water amounts withdrawn by the drains indicate the GW infiltration into the UIs; despite being a limited amount of water, quantification of the water amounts is important to compare them with the results of the other scenarios in the framework of urban underground management. Model percentage discrepancy is considered to be low (4.59 × <sup>10</sup>−3). A good coherence was detected between the drain outflows and the mass balance with neighboring zones (i.e., surrounding aquifer and the UIs), thus validating the obtained results.


#### *4.2. Modeling Scenarios*

Groundwater inflow for all the UIs was calculated, and results are summarized in Figure 10. As can be seen in Figure 10, an absence of inflow was detected for some subway line branches, as the water table level was lower than the bottom of the UIs [29]. Particularly, these inflow gaps were visible in the north, along subway line M1, and in the central portion of the domain, close to Cadorna exchange station (subway lines M1 and M2). The tunnel sections more exposed to GW inflows are the westmost stretch of M1, towards Bisceglie Station (M1-a), and the stretches close to Uruguay station (M1-b) and between QT8 and Lotto (M1-c) for line M1; and the sections from Porta Genova to Sant'Agostino station (M2-a) and from Lanza to Moscova (M2-b) for subway line M2. Due to their major depth, subway lines M3 and M4 were completely submerged by the water table, which also occurred for subway line M5 (Table 4). With regard to public car parks, 34 out of 67 resulted in infiltration; in the central area, Washington/Piemonte (P-a) and Carducci (P-b) turned out to be among the most impacted infrastructures, while, for example, Betulle Est was impacted in the west (P-c). Critical sections for M1 and M2 were already identified as areas where a historical interaction (i.e., submersion) with the water table was evidenced [46,83]. In particular, Sant'Agostino (M2-a) was impacted for both GW minimum and maximum conditions.

As summarized in Table 4, groundwater inflows for S1-S3 are limited, with low orders of magnitude. The highest values of inflows (10<sup>−</sup>5/10−<sup>3</sup> order of magnitude) were detected for the oldest subway lines M1 and M2, modeled with higher conductance values to simulate the absence of waterproofing systems and a progressive saturation of the walls over time. As for the deepest lines, such small values are attributable to the low conductance representing lining systems. The spatial distribution of these infiltrations is different, as for shallow lines the infiltrations are detected only at certain spots, as visible in Figure 10, thus evidencing local but more critical situations to manage.

**Figure 10.** Areas showing GW infiltrations into UIs.

**Table 4.** GW inflows into UIs (m3/d) for S1–S3. Please remember that for M3, M4, M5, and parks, K was always set equal to 1.16 <sup>×</sup> <sup>10</sup>−<sup>3</sup> m/d. Percentage below the water table is intended as the sections of UIs where the bottom of the infrastructure is lower than the hydraulic head.


At the most critical points highlighted in S1–S3 (Figure 10) for subway lines M1 and M2, and for some public car parks, locally punctual wall fractures have been simulated to quantify the variation in GW infiltrations. The results of these spots are summarized in Table 5. As is visible, the most critical effects, also considering the features of the UIs (i.e., depth, volume), have been identified for M2-a, around Sant'Agostino station. The infiltration for these points generally increased linearly to one or two orders of magnitude.

As for the two public car parks under construction (Figure 1b), an absence of infiltration was detected in both cases, with respect to the considered groundwater maximum condition, due to a lack of interaction with the water table (Figure 11) which was contrastingly evidenced for the close branches of subway line M4.

**Figure 11.** (**a**) 3D geographical setting of the area close to the car parks currently under construction. The car parks and the subway line M4 are visible below the road network; the names of some roads are indicated to provide more geographic details. Three-dimensional underground reconstruction of (**b**) Brasilia car park, (**d**) Scalabrini car park, (**f**) Lorenteggio 124 intervention point. GW/UIs interaction for (**c**) Brasilia car park, (**e**) Scalabrini car park, (**g**) Lorenteggio 124 intervention point. (**b**,**c**) refer to point 1 in (**a**); (**d**,**e**) refer to point 2 in (**a**); (**f**,**g**) refer to point 3 in (**a**). Transparency has been adopted to represent the volumes submerged by the water table; as visible in (**c**,**e**,**g**) this occurs only for subway line M4, and not for public car parks. The red arrows indicate the viewpoints and the view directions adopted in the 3D visualization of the subsurface elements. Images were realized using ArcGIS Pro.

**Table 5.** Comparison of GW inflows into UIs (m3/d) for the initial scenario (S1–S3, intact walls) and their corresponding final scenario (S4–S6, leaky walls). Please remember that, for car parks, K was always set equal to 1.16 <sup>×</sup> <sup>10</sup>−<sup>13</sup> m/d for S1–S3 and to 1.16 <sup>×</sup> <sup>10</sup>−<sup>9</sup> m/d for wall fractures in S4–S6. S means station, T means tunnel, P means park. Depth (m) has been provided for subway stations and parks, as they are designed from the ground field; as for tunnels, since they are not designed from the ground field, thickness was provided rather than depth.


#### **5. Discussion**

Managing GW/UIs interaction in urban areas is a challenging issue. Different problems can arise regarding GW quality, quantity, and thermal issues [5,6,15], but stability, erosion, and infiltration for UIs are some further topics to consider. With regard to GW infiltration into UIs, the scientific literature deals both with water inrush calculation during the construction of tunnels [26,27] and problems regarding already-operating underground tunnels [31,33]; in this study, a local scale numerical model was developed for the western sector of Milan city, applying a methodology to quantify GW infiltrations into completed and operative UIs.

#### *5.1. Modeling Scenarios*

Model results in terms of calibration were generally acceptable (Table 2, Figures 8 and 9). However, some head targets did not show an optimal result. This happened for a couple of targets in the north-western portion of the domain, and for one target in the south. In the north-west, not far from the critical targets, the behavior of the water table is presumably influenced by a multitude of local situations. The proximity of a group of quarries and a public-supply well field with high discharge (Figure 4), and the presence of clay lenses determining the existence of perched aquifers with seasonal oscillations [51,52] make predictions more uncertain. This response could highlight the presence of local mechanisms, possibly uniformed by the targets, that have been neglected. In this sense, the model provides a guide for future data collection, that could allow the improvement of the appropriateness of the conceptual model [84]. In addition, acquiring further data

could contribute to making more effective predictions, thus improving the model use in supporting management decisions [85].

Identifying the areas more exposed to infiltrations is important to predicting future risks due to a more severe water inrush; thus, adopting strategies to ensure these infrastructures are preserved is vital [86,87]. Some of the impacted areas (i.e., M1-a, M1-b, M1-c, M2-a) have already been identified as critical in previous works [46,83]. In this case, only a qualitative GW–UIs interaction was detected through a GIS methodology. This spatial coherence among the results could be considered as a validation of the numerical model. At the same time, model findings could represent a step forward in the definition of the urban conceptual model; through this approach, GW infiltrations resulting from GW/UI interactions could be estimated. As for M2-a, P-a, and P-b, the highest depth of downtown infrastructures (Table 5) plays a key role in influencing GW/UI interactions. This is due both to a high population density, thus requiring more space for subsurface infrastructures [88], and to the adoption of specific construction methods; as an example, Sant'Agostino station was built with two overlapping pipes [61]. As for the western sector, the complex geological situation explained above could be a possible driver of the infiltrations both for subway lines (M1-a) and underground car parks (P-c), despite their limited depth (Table 5). To counteract this situation, in the framework of creating a more sustainable and resilient city [5], some residential constructions have been designed with superficial car parks occupying the first floors of the buildings. With regard to public car parks, the new buildings currently under construction have been designed as two floors deep; at this time, this results in an absence of impact even considering a groundwater maximum condition (Figure 11). However, prolonged monitoring should be useful to cope with the evolution of GW/UI interactions. Finally, in the north, a reduced GW/UIs interaction is attributable to a wide unsaturated thickness of the shallow aquifer [51], with the water table located around 10–12 m from the ground field.

#### *5.2. Considerations of the Adopted Modeling Approach*

The applied modeling strategy aimed to quantitatively evaluate the interaction between the GW system and the subsurface structures. With regard to the calibration process, it is not tied to the prediction of interest; in fact, it is based on head targets whose hydraulic measures are not directly connected to the final goal (i.e., GW infiltrations into UIs). The information content on which the head targets are based is not informative about the degree of connection between the UIs and the water table. In technical terms, the K of the walls is completely in the null space and outside the solution space of the model. This does not mean that the calibration is useless, but it does mean that the model could not be so much a predictive tool as a way to understand a phenomenon (inflow across leaky walls) in general terms. The geometry of the UIs is realistic [46], but one can only make hypotheses about the permeability of the intact and leaky walls that are not in any way informed by the calibration. To limit this uncertainty, a literature analysis was conducted to choose the initial conductance values for subsurface impervious structures [78] and the conductance to simulate isolated fractures [80]. Moreover, an ensemble of scenarios [89] was defined to deal with non-lined systems, testing different conductance values. In this way, stakeholders are enabled to visualize a range of impacts and they could consider them to apply different management options [90,91]. As for S1–S3 (Table 4), GW infiltrations are very limited, especially for waterproofed subway lines; thus, the model allowed for the gaining of insights into the conductance values that are needed to simulate an almost impermeable element.

Anyway, obtaining good calibration results was crucial, since they allow GW/UI interactions to be well represented and, consequently, they allow the obtaining of a more reliable estimate of the infiltrations originated by the relationship between the aquifer and the subsurface infrastructures. As visible in Figures 8 and 9, this is mostly true for this specific case, especially for the targets located in the central part of the domain that lie in proximity of the main UIs' elements.

Using MODFLOW-USG as numerical code allowed the refining of the grid horizontally, therefore properly representing the UIs. Moreover, through the implementation of the unstructured grid, the key numerical computations could be limited within the required bounds [50], making the simulations less computationally intensive. Above all, the adoption of MODFLOW-USG was pivotal to model the UIs, as it allowed the WALL package to be used to represent not only the cells' lateral sides through HFB, but also the top and the bottom of the UIs. In this way, the subsurface elements could be modeled with their real depth and volumes, thus refining previous applications of the HFB package to simulate UI fully penetrating single-layered models [37,43,44]; hence, a precise estimate of further modeling aspects (i.e., evaluation of the barrier effect on groundwater flow paths) should also be guaranteed. In MODFLOW-USG, to reduce numerical instability, desaturated cells (i.e., dry cells) are not inactivated, so there could be a small amount of flow from one cell to another. The adoption of the DRN package helped to solve this possible issue, especially in the upper portion of the domain where unsaturated aquifer was present. As a drain is activated only when the hydraulic head is at least equal to the drain elevation, it was possible to unravel where an effective infiltration was present. The choice of the DRN package also came after its previous applications to quantifying flooding episodes during the construction of tunnels [34,36]. Through the developed methodology, modeling GW/UI interactions could be enhanced. In fact, combining the use of HFB, DRN, and mass balance zones to quantify infiltrations depending on different conductance values is possible, instead of deactivating cells of impervious structures. Thus, a step forward could be taken in the development of the urban conceptual model, supporting previous approaches conducted within the same domain [92,93], or in other areas [94] where different aspects of GW/UI interactions have been investigated but GW infiltrations into subsurface elements were not quantified.

The methodology has been tested on a steady-state numerical model. Future applications on transient numerical models would be possible depending on long-term data collections [95]; this could raise awareness about infiltration issues, supporting a deeper interpretation of GW/UI interactions and making the model a useful management tool to make long-term predictions [84].

#### *5.3. Decision Management*

The infiltration issue of UIs in Milan city is historical. Different episodes have been documented over time [46,60,83], leading both to economic and management problems for Metropolitana Milanese Spa, the subway managing company. For example, the section between Piola and Lambrate stations, along subway line M2 (outside the numerical model domain), was closed during summer 2019 to complete lining works because of GW infiltrations, thus forcing the use of surface public transport. Although the water inflow is small with respect to water inrush into subway tunnels during their construction [28,80,96], this situation could trigger further issues over a long time period (i.e., corrosion of foundations), resulting in a decline of the subway system efficiency; thus, this problem should not be underestimated.

To ensure sustainable development of GW/UI interactions, effective engagement of the stakeholders should be of great value [97–99]. Open communication is needed to raise awareness about the importance of data to describe the system and conceptualize and develop a model [91,100] with increased predictive capabilities. For this specific case, monitoring, estimation, and control are essential aspects for tunnel management [96]. Having access to existing infiltration measures, if available, or implementing monitoring of the punctual inflows along the tunnels or for car parks would also improve the calibration process; in this way, model uncertainty would be reduced, thus strengthening the usefulness of hydrogeologic models for decision-making bodies [84,85]. The collection of field data could focus on the most critical sectors highlighted (i.e., M1-a, M2-a) by the model results. Amongst these areas, dewatering solutions could be adopted to manage the issue, thus contributing to preserving the status of the subway network, avoiding the development of

more serious issues as occurred for the surrounding areas of Piola and Lambrate stations. In particular, the historical issues of Sant'Agostino station, also due to the adoption of specific construction methods [61], impose an increased degree of attention for this limited branch of subway line M2.

However, applying these solutions would be a consequence of effective GW infiltrations into UIs. A move away from reacting and correcting measures, focusing on preventive actions [6] to secure the UIs, should be evaluated. In a previous work by Sartirana et al. [51], underground car parks were classified as possibly critical for different GW conditions if the difference between the reference plan (i.e., bottom) of the UI and the water table was less than one meter. To avoid infiltration issues, activating localized pumping when a certain threshold is locally exceeded would be a possible measure [101]. To do so, early warning monitoring solutions, such as integrating GIS, BIM, and GPS techniques [102,103], with continuous online data measurements should be implemented in proximity of the most critical UIs.

Moreover, groundwater is not only an annoyance for its side effects, but it is also a heritage [6] in urban frameworks; therefore, further management strategies could be proposed. For example, as GW is a valuable energy reservoir [15,93], increasing the adoption of GWHP systems, possibly only due to extraction wells, could keep the water table levels controlled close to the UIs, thus not only limiting the infiltration issues but also exploiting the thermal potential of these subsurface elements [104].

Finally, in the framework of the goals of the Plan of Government for the Territory, this local-scale urban model could help the decision makers to understand and manage the relationship between new UIs and water table levels, testing possible urban underground development scenarios.

#### **6. Conclusions**

This work aimed to adopt a methodology to quantify GW infiltrations into UIs (subway lines and public car parks) with the view of assisting urban underground management. In this sense, the realization of a local-scale, urban numerical model allowed the following:


The overall findings of this study could provide a useful tool to the stakeholders to properly design new UIs in the framework of the planned underground development of the city. In this sense, the numerical model could be used to realize different GW scenarios, testing their effects on the designed UIs. Furthermore, modeling their tops and bottoms through the HFB package could improve the evaluation of their barrier effect on groundwater flow paths. For future applications, reasoning the combination of the HFB package with different third-type boundary conditions (i.e., River, GHB) to model other subsurface elements (i.e., sewer systems, buried channels, etc., to evaluate their leakance) could represent a challenging task. The methodology has been tested for the city of Milan—nonetheless it should be worth considering its application to other urban realities to enhance the analysis of GW/UI interactions.

**Author Contributions:** Conceptualization, D.S.; methodology, D.S. and T.B.; validation, M.R., M.D.A., L.F. and T.B; formal analysis, D.S. and C.Z.; data curation, D.S. and M.D.A.; writing—original draft preparation, D.S.; writing—review and editing, D.S., C.Z., M.R., M.D.A., M.C., A.R., L.F. and T.B.; visualization, D.S., C.Z., M.C. and A.R.; supervision, M.R., M.D.A., L.F. and T.B.; project administration, T.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research did not receive any external funding.

**Acknowledgments:** The authors are grateful to Metropolitana Milanese S.p.a for providing both the altimetric profiles of the subway lines and the piezometric data that were used as information for model calibration. Moreover, the authors would like to warmly thank Daniel T. Feinstein of USGS and Gennaro Alberto Stefania for their support and suggestions during the development of the work. The authors would also like to thank the three anonymous reviewers for their comments, which helped to improve this article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

