Transient Flow and Transport Modelling of an Historical CHC Source in North-West Milano

Legislative Decree 152/2006 requires Public Authorities to identify the subjects who are responsible for soil and groundwater contamination. In highly urbanized areas with a long industrial history and an elevated number of potential contaminant sources, as in N-W Milano Functional Urban Area (FUA), their identification can be difficult. Since the groundwater flow has showed consistent changes in the last 30 years as in Milan, the problem became even more complicate. The Public Authorities put in charge by the law, i.e., Regione Lombardia and Città Metropolitana Milanese, need new methodologies to assist them in finding the source locations and implementing remediation actions. The aim of this study is, coupling unsteady flow with fate and transport model of Chlorinated Hydrocarbons, to reconstruct the potential impact of a former chemical plant on public wells in the N-W area of Milano. The proposed methodology consists in (a) reconstruction of the piezometric trend over time (1980–2018) by means of a transient flow model (MODFLOW2005 + Parameter Estimation PEST) and (b) simulation of transport as a function of the flow variations in time. The obtained results were compared with the previous ones obtained with a quasi-steady model (no changes in time-dependent parameters). Finally, a predictive scenario was performed to assess the potential evolution of tetrachloroethylene (PCE) in groundwater; on this frame, strategies to monitor and remediate the contamination were proposed.


Introduction
Groundwater contamination by Chlorinated Hydrocarbons (CHC) has been detected since the 1960s in the main aquifers (shallow and semi-confined) in Milano Functional Urban Area (FUA), but often the sources (named hot-spots in [1]) are very difficult to locate and quantify over time. Among the possible methods to identify potential sources, there is the so-called "more probable than not", which may be applied considering two main principles as reported in [2]. Since the end of the II World War, the FUA in Lombardy Region is one of the most urbanized areas in Europe ( Figure 1). The Northern area of FUA has been characterized by a dense agglomeration of industrial plants; in and around the city of Milano, automotive, refining, chemical, and steel plants have been established from the 1950's [3]. Because of the groundwater flow direction, the high hydraulic conductivity and the high groundwater withdrawal rate (yearly average public withdrawal is around 250 Gl/y), Milano represents a groundwater drainage area, that attracts contaminated groundwater from the surrounding areas. Persisting critical issues deriving from past industrial activities represent environmental problems that the Municipality of Milano is facing nowadays: on one side, degradation of drinking water quality due to the capture of plumes by supply wells; on the other side, the cost of extracted water cleanup interventions, which, in the absence of an identified polluter, befalls on the citizens. In accordance with European Legislative Decree (ELD), in application of the Polluter Pay Principle (PPP), artt. 242 and 244 of Italian Environmental Code (IEC) [4], establish that in case of soil/groundwater contamination, the chemical characterization of the site and the remediation measures (if necessary) must be imposed by the competent Authority (in the Italian legal framework the Provinces/Districts) only to the liable operator according to his/her fault for negligence, imprudence, in-experiences on the basis of the evidence of a causal link. This frame depends not only on the IEC but also on the Italian Civil Code, and in particular on art. 2043 of Italian Civil Code (IT, 1942 and modifications/integration) on civil liability. The PPP specifies that the identification of the responsible for contamination is a compulsory task to be carried out by the Public Authority. In highly urbanized areas, this legal framework is difficult to apply because of (1) the presence of many industrial plants active from the early '50 (some of them have since ceased all activities and the presence of pollution down-gradient cannot be easily linked with their historical position) and (2) the groundwater flow is influenced by long-term modifications due to variations in groundwater abstraction for drinking and industrial purposes [5]. This latter, which produced high variations in flow directions over the latest decades, is known with a high degree of uncertainty. Previous works focused more on developing transport models rather than on carefully reconstructing the groundwater flow directions over time [6]. As numerical models are very important to support the management of groundwater [7][8][9] and both flow and transport are useful in reconstructing the historical contamination of a hot-spot in hydrogeology, the main aim of the present study is to demonstrate that consider flow direction in time is very important for reconstruct the contamination path and provide stronger tools for Public Authorities to identify the polluters. In order to pursue these goals, the paper (1) simulate a 3D groundwater flow with MODFLOW2005 [10] in transient mode (from 1980 to 2018), (2) simulate the fate of dissolved CHC contaminants released by a historical source with MT3DMS [11] (comparing the results previously obtained by a steady-state model [6]), and (3) forecast the dispersion of CHCs in the coming years, thus providing information to Water Managers and the Public Authority that need actions against the polluter and to plan future groundwater abstraction scenarios.

Study Area, Hydrogeology, and Conceptual Model
The study area (also called Pilot Area or PA) is located in the central part of Lombardy Region, and it comprises the North-Western area of Milano City and several surrounding municipalities (coloured in violet in Figure 1). The extension of the Pilot Area is about 395 km 2 (21 km × 18.8 km). The area is densely populated: the 2014 census reports a total of 617,773 inhabitants. The main hydrography is represented by Seveso and Olona rivers and the Villoresi irrigation canal. Stratigraphic data are available from drilling data collected during different studies [6,[12][13][14][15][16] during which the collected sediments were classified and used to interpolate (Leapfrog software [17]) the main aquifer geometry and hydraulic properties. In this phase of the work was followed the Lombardy Region Aquifer classification [18] that considers four different hydrostratigraphic unit called A, B, C, and D, originated thanks to the overlay of plio-pleistocenic alluvial sediments [19]. The main aquifers affected by the contamination are the shallow (Aquifer A, mainly of sandy-gravel composition) and the underlying semi-confined aquifer (Aquifer B), mainly composed by sand interrupted with clay lenses that subdivide it in several silty layers that locally confine this aquifer.
Firstly, by using collected sediments classification, more than ten cross sections were developed (from N to S and from W to E). The structure thus reconstructed was subsequently translated into a MODFLOW2005 3D model. The model only represents the shallow and semi-confined aquifers (called A and B), which are impacted by CHC contamination. The underlying aquifers (Green and Brown), represented in Figure 2, were not represented considering that the thick clay layer (>5 m) between Aquifer B and C is able to confine the Aquifer C and then is hydraulically separated. A North-South and a East-West cross section of the area are depicted in Figure 2 and their location is shown in Figure 1b.
Aquifer A is delimited by the topography surface on top and by a discontinuous clay-silty layer, separating it from aquifer B, on bottom. In Pilot area its average thickness is about 40 m and it is mainly constituted by gravel with presence of sand and a few clay lenses. The clay-silty lens divides hydraulically the two aquifers in the Southern part of the area and it is absent in the North, so that the two bodies constitute an undistinguished aquifer with a saturated thickness of 80-100 m. Aquifer B is separated by the underlying aquifer C by a layer of clay. Its thickness increases towards South from about 40 m to 60 m. The Northern part is mainly constituted by coarse sediments (gravels and coarse sands) while towards the South, an increase in the sandy fraction is observed. The wide presence of clay-silt lenses makes this hydrostratigraphic unit a multilevel aquifer.

Numerical Modeling
The work presents a 3D numerical modeling of transient flow governing the transport of a PCE plume from a historical source using MODFLOW-2005 [10] and MT3DMS [11] implemented through the graphical interface Groundwater Vistas 6. The hydraulic conductivity distribution (K) was extracted from a steady state model (concerning the same domain) calibrated using PEST [20] based on the piezometric heads monitored in 2018, as described in Appendix A. Then, 3 new models were implemented:

1.
Steady flow simulation (1980s) calibrated using PEST [20] (abstraction wells and recharge parameters) in order to have initial condition for transient model.

2.
Transient flow simulation from 1990 to 2018 calibrated using PEST [20] (abstraction wells and recharge parameters) in order to have initial condition for transient model.

3.
The simulation of a suspected source affecting the quality of groundwater in N-W Milano city taking into account the advective-dispersive transport equation.
The model grid has 210 rows and 188 columns (uniform cells with dimension 100 m × 100 m for a total of 355,320 active cells) and 9 layers (shown in Figure A1). The cell size is chosen to be the least in order to maximize the result accuracy, compatibly with the computational capability of the software. The first layer represents Aquifer A. The second layer represents the discontinuous aquitard that locally confines the Aquifer B. The underlying seven layers represent Aquifer B. The detailed representation of the semi-confined aquifer is due to better represent the fate of the contaminant towards deepest aquifer portion that is mainly exploited by public wells.

Temporal Discretization
The simulation time is 29-year (from 1990 to 2018), subdivided into time periods (named "Stress Periods-SP" in Modflow) of different lengths, as shown in Table 1. The calculated heads from a steady-state model (called SP0) representing the 1980's has been used as initial heads for the transient model.

Boundary Conditions
The boundary conditions (BCs) (shown in Figure 3) have been applied as follows:  • Villoresi is an irrigation canal whose bed is made of cement and it is considered almost impermeable (Hydraulic conductivity of 4.32 mm/day was estimated in steady flow model) • The Olona and Seveso rivers paths after the 1960's were modified because the urbanization development and large part of their bed were sealed using cement (Hydraulic conductivity of 8.64 mm/day for river sections presenting bed sediments) (c) Recharge (Neumann condition) applied through 10 zones divided into 3 different uses of soil: • Green areas with water infiltration computed by Thornthwaite-Mather's method [21], using data from a total meteorological station "Milano Lambrate" (it is the unique station that monitors rainfall from 1900 to nowadays) • Agricultural areas considering an uniform irrigation equal to 1.6 mm/day [1,9,22] • Urbanized areas with estimation of the losses of the water network of 0.80 mm/day [1] which is 15% of the circulating water in the aqueduct system (quantified by the Water Municipality Manager of Milan) (d) 1077 (Neumann condition) analytical wells screened in different layers and with different discharge in time simulation.
Concerning the latter, to recover all the data about the withdrawals (well registry info, latitude and longitude coordinates, screened levels and abstraction well) several existing databases were compared (Sistema Informativo Falda and Catasto Utenze Idriche, WebGIS Ambiente Comune and Censimento Pozzi Comune di Milano) and some historical documents were consulted [23][24][25]. The flow model calibration was done in transient conditions. Hydraulic conductivity was not included among the calibration parameters, because they were deemed to be already sufficiently defined in the [19] model (further details are in Appendix A). Other parameters were optimized instead: • CH for each SP was set with a "Trial and Error" procedure using the available head observations • RIV package parameters were kept the same of the steady model because their leakage represents less than 1% of the total budget due to the low hydraulic conductivity of river and canal beds • Recharge and withdrawals have been calibrated in transient state through PEST Calibration of well withdrawals is described in detail in the next Section 2.4. Concerning the effective porosity, it was considered uniform within the homogeneous zone, with a value of 0.2 for the permeable zones and 0.01 for the clay lenses. This parameter was initially involved into the calibration process, but it was observed that there was a very low influence on the objective function during PEST runs. Thus, the porosity was not included in the calibration set.

Abstraction Calibration: WMG Software Coupled with PEST
The automatic calibration was carried out for both recharge and well withdrawals due to the high uncertainty related with some wells (i.e., the flow rate, the active interval time of wells). In particular, during calibration process, the wells were divided in 63 groups according to Table 3: 1.
The surrounding municipalities of Milano where the information of the private/public wells are affected by higher uncertainties, the factors (parameter multiplier) of calibration interval for flow rate are imposet to be more relaxed (from half to double of the starting abstraction) 2.
Milano and public pumping station (within the area of interest) where the wells information present less uncertainties (historical monitoring of Water Municipality Manager of Milan); the factors (parameter multiplier) of calibration interval for flow rate are closer.
Then, each withdrawal was multiplied by the calibrated multiplication factor relative to its belonging group (Equation (1)).
where m i is the calibrated multiplication factor relative to the group to which the i-th well belongs, Q initial,i is the initial withdraw value of the i-th well and Q calibrated,i is the calibrated withdraw of the i-th well. To carry out the calibration, the code PEST [20] was used, coupled with an external software (WMG), that allowed to prepare the files needed for PEST to perform the optimization process. WMG works independently with steady-state or transient models. The program flow-chart is shown in Figure 4.

PCE Transport Modeling
The transport model was based on the transient flow solution described in the previous section. Concerning the setting of hydrodispersive and chemical parameters, the longitudinal dispersivity α L was 20 m, the transverse α T and vertical dispersivity α V were 1 and 0.03 m, respectively (calibrated in previous work [6,9,22] where the transport model was set in a site-specific area and with many available data). The values for distribution coefficient k d and bulk density were 0.000113 m 3 /kg and 1700 kg/m 3 . No degradation reactions were simulated as the aquifers are oxygenated and PCE dehalogenation was not expected. For this reason, a relatively low value of half-life time (10 years) was considered. Considering the transport boundary conditions, the PCE sources were imposed through a Dirichlet boundary condition (i.e., constant concentration) inserting a suspected historical source. The location and PCE concentration of the source were considered taking into account (a) the historical data inside the contaminated site (from 1954 to 2018), (b) the screening position of the monitoring wells between the 1 and 5 layer, and (c) the presence of an active barrier from early 1980 with a known extraction rate. For solving the transport equation a numerical method (TVD-Total Variation Diminishing) was used in order to consider properly time stepping during calibration transport simulation.

Model Settings for PCE Scenarios
In order to simulate a PCE transport forecast Scenario, 3 new SPs are added to the model (representative of the time periods 2018-2021, 2021-2027, and 2027-2040 respectively). The aim of the transport modeling for the Scenario is to understand what relevant impacts will be on sensitive receptors, i.e., water supply wells. The results will help to provide possible remediation strategies and to design an effective monitoring network, since the suspected area lies up-gradient of some public drinking well stations. The transport model for the Scenario considers that is achieved a full remediation of the source zone and simulates the natural dispersion of the plume in saturated zone (10-years half-life time degradation considered). No modification of the flow settings were made with respect of SP5, whereas some changes have been made on the transport model:

1.
Constant concentration of source was set to 0; 2.
The modeled concentrations at 2018 were imported as initial conditions concentrations for the Scenario. Table 4 describes the modeled Scenario divided in three different time-frames.

Transient Model: Calibration and Validation
The simulated heads in the overall area are reported in Figure A3 for SP0, SP2, SP5. The simulated contour map shows that the groundwater flow changed over time (in Milano the flow is strongly perturbed by pumping well stations and locally the flow becomes radial towards Milano City). During the calibration process, weights were applied to head targets. They were calculated as inverse to the standard deviation of each target measurement, according to [27]. Standard deviations were estimated according to [28]. Pumping wells and piezometers close to wells were given higher measurement standard deviations, i.e., lower weights, than those attributed to monitoring wells in undisturbed areas. This choice accounts for the high uncertainty tied to not knowing if the measured level was affected (or not) by pumping. A standard deviation of 2.5 m was attributed to pumping wells; a standard deviation of 1 m was attributed to undisturbed piezometers. An intermediate value of standard deviation and weight was attributed to piezometers close to pumping wells. Overall model performance based on the available head observations was evaluated with several statistics, summarized in Table 5 and graphically presented in Figure A2. The error of the model mass balance is 0.762 percent and the scaled root mean square is less than 3%, highlighting the goodness of the obtained fit.  Figure 5), a decreasing trend is observed and it is linked to a more reliable data set available in the last years. The error among the different SPs and the whole model are below the reference values (<10%, [28][29][30]).
For the model domain, a total groundwater flow rate between 14 and 20 m 3 /s was estimated, whereas the total flux for the area of interest is represented by the 25% of these values. For the whole model domain, Table 6 reports the inflow-outflow components for each SP.

Comparison between Time Series Data Simulated and Measured
In order to better evaluate the model performance nearby the suspected source before its simulation in the transport model, a comparison between simulated and observed heads in time for some observation points was done ( Figure 6). Some points are spread around the model domain (Pz265, Well_3A and Well_3B) while some others are close to the suspected source (FOG57, FOG56, and Pz5), as shown in Figure 7. Piezometers FOG57 (SIF code 0151461496), FOG56 (SIF code 0151461495), Pz265 (SIF code 0151461265), and Well_3A (SIF code 0151540003) are screened in the Aquifer A (shallow) and piezometer Pz5 (SIF code 0151700005) and Well_3B (SIF code 0152040003) are screened in the Aquifer B (semi-confined). As shown in Figure 7, the calculated heads in monitoring points outside the area of interest show some discrepancies with the observed head time series; within the area of interest, adherence to observations is higher, even if the graphs still show differences of some meters in one or more SPs. Considering the complexity of the model, the wide period of time considered, the targeted average head values considered and the uncertainty tied to the observations, the results were considered acceptable in both the shallow and deep aquifer with regard to the main direction of groundwater flow affected by abstraction during the whole simulation period.

Comparison between CHC Plumes Simulated in 2014 and 2018 Model: Differences and Upgrade
The comparison between the two models shows multiple differences between the simulated plumes. First of all, there are several differences in the modeling scheme that must be taken into account, described in Table 7. The upgrade developed within the 2018 model consist in adopting a 64-year transient scheme for flow and transport. The effort made in reconstructing the time-varying condition is necessary to better represent the real contaminant path over time. Generally, the plume modeled in 2018 spreads more widely than the older one, due to the flow variations over time and the different discretization grid dimensions. As shown in Figure 8a,b, the simulated transports in [6] and AMIIGA model in the first stage (1989) differ, meaning that the contaminant path is strongly dependent on flow direction. Looking at [6] last SP, the AMIIGA model well simulates the contamination that affects the area among San Siro, Novara, Chiusabella and Cimabue pumping stations (Figure 8c,d). In [6] model, such contamination has been modeled as coming from two different sources, one located in the suspected source and one unknown source not directly linked to a contaminated site but in a disposal area. Further studies show that the second source does not contribute to the modeled contamination and the AMIIGA model is able to simulate it considering a single suspected point source. The new model also shows a behaviour (close to where the secondary source in [6] is located) that, although indicative of the presence of a secondary source, is instead a consequence of the denser vertical discretization introduced to better represent the local clay lenses in the model. The discontinuities of such lenses modify the transport directions, allowing the contaminant to flow vertically and to deepen within Aquifer B through more permeable areas. In conclusion, the new model indicates that the sensitive receptors of CHC contamination in the time period between 1980 and 2018 are several public wells and the pumping stations of San Siro, Novara, Chiusabella, and Cimabue. The model results were validated comparing the simulated concentrations with the field-measured data (shown in Appendix B, Figure A4).   Figure 9 shows the results of the Scenario transport model respectively for the three different time-frames in shallow aquifer (layer 1) and in the deeper part of the semi-confined aquifer (layer 7), while Figure 10 represents the breakthrough (B-T) curves relative to three monitoring wells (Monitoring Well 07, 99, and 19) along the plume's path (their locations are shown in Figure 6). In the freatic aquifer, the shallow contamination proceeds, in low concentrations, towards the principal flow direction for about 2 km between 2018 and 2040, without influences of pumping wells abstraction and the expected concentrations will be close to drinkable threshold limit (i.e., 12 µg/L in 2040 is close to the 10 µg/L drinkable limit [4]).
As expected, the removal of the PCE source leads to a more rapid and complete remediation of the groundwater plume near the suspected source and in general in the shallow aquifer. In the upper part of the Aquifer B (layers 3 and 4) the B-T curve of the Monitoring Well 19 (down gradient of San Siro pumping station) shows that part of the contaminant passes over San Siro pumping station, with low (around 3 µg/L) but increasing concentrations. Differently in the deeper part of the semi-confined aquifer (layers 5 to 9), the model shows that the contaminant is stopped due to a "capture zone" [31] induced by the superimposition effects of San Siro, Novara and Cimabue pumping stations, with a reduction of the up-gradient plume length of about 3 km between 2018 and 2040; in San Siro, the contamination is still present in 2040 with concentrations approximately close to 100 µg/L. In general, after 22-years simulation, the higher concentrations will mainly affect the deeper part of Aquifer B where most drinking wells are screened, rather than Aquifer A. The impact on Aquifer B derives from the (a) leakage through San Siro wells in layers 3 and 4 due to the low abstractions rate and (b) from a slow release due to the PCE trapped in the aquitard, that appears as a "secondary source" which is difficult to be remediated.

Advantages of Modeled Scenario
The modeled Scenario considers the total source removal. Based on its results, a monitoring network is needed to check the contamination spread in the deeper part of semi-confined aquifer. In Milano city a dismissed monitoring network is present and in order to monitor the attenuation of contamination in shallow aquifer it could be restored. In addition, monitoring wells could be drilled in both aquifers (for example coupling shallow and deeper piezometers), as shown in Figure 11, that will be able to monitor the concentration of the plumes in time.
The advantages for Authority and Water Managers taken from the modeled Scenario results are mainly: (a) to detect possible leakage of the contamination towards the other pumping stations not yet affected and (b) to monitor the contamination among the affected pumping stations. Furthermore, since the model shows that in 2040 (SP8) the contamination will be captured in Aquifer B close to San Siro, Novara and Cimabue pumping stations (with a peak value higher than 100 µg/L in San Siro), some remediation actions can be suggested and designed to reduce the contamination that could affect the water supply wells and consequently the higher cost of water treatments: • Increase the abstraction rate in some of the wells of pumping stations that will act as hydraulic barrier, alleviating the PCE load that have to be removed by Water Managers • Evaluate bioremediation actions in the area directly involved by the plume propagation • Reactivate some existing deeper wells along the plumes propagation to use them as hydraulic barrier.
(a) (b) Figure 11. Proposed monitoring networks for (a) Aquifer A and (b) Aquifer B.
Regarding Aquifer A, the model shows low solute concentrations (below 12 µg/L), therefore no remediation actions are strictly required as the Scenario represents "worst-case".

Conclusions
In highly urbanized and industrialized areas, where multiple contaminant sources are present (Milano N-W FUA), a robust modeling approach is needed in order to support groundwater management. Under this aim, the paper wants to demonstrate that, in order to simulate fate and transport of a historical source, the groundwater model has to take into account the flow direction variation during time. In fact, in Milano N-W FUA, the superimposition of both industrial site abstraction in surrounding municipalities and the pumping of water supply wells due to the high water demand in Milano city have a huge influence on the main flow direction and consequently on the track of contamination that originated from the area of interest. Due to the necessities to reconstruct historical sources and the assessment of potential contamination receptors, within the present study (1) a 3D transient groundwater model is presented, built after a detailed geological description of the aquifers interested by the contamination, that is able to simulate the flow in the period between 1980 and 2018; (2) The model has been used to describe the fate of dissolved CHC (PCE) and the results have been compared with a previous model [6], highlighting the improvements in representation of the observed concentrations; (3) Scenario in different time-frames (2021, 2027, and 2040) have been performed to understand the plume evolution.
The model results show that in the deeper part of Aquifer B (layers from 6 to 9), the contamination is captured by the San Siro, Novara, and Cimabue drinking well public stations, while in the shallow aquifer the contamination proceeds following the main flow direction (S-E). Observing the evolution of the plume over time, a new monitoring network has been designed starting from an existing one, that enables the Authorities to assess the concentrations in the area of interest and preserve the quality of abstracted drinking water. The high PCE concentrations currently observed and the long time required to remediate the load of contamination, suggest also to propose some actions aimed at reducing the impact to the San Siro, Novara, and Cimabue pumping stations and to avoid further spreading of the PCE contamination. This work represents an example of the application of flow and transport modeling to support the Authorities aimed at improving the groundwater monitoring network and remediation actions. Moreover, considering a transient transport modeling based on a transient flow model, the work demonstrated that the obtained results are a good upgrade with respect to the previous work. The methodology can be applied to several different suspected sources.

Funding:
This research received funding from Interreg Central Europe "AMIIGA Project"-Grant CE N • 32.

Acknowledgments:
We thank Regione Lombardia and Environmental Protection Agency for providing the information about contaminated site and about extensive available dataset of PCE concentration.

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

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Hydraulic Conductivity Calibration Process
Calibration of hydraulic conductivity has been carried out against the 2018 [19] head observations (updated CH and updated vertical recharge) through the definitions of 337 Pilot Points (PP) and 4 homogeneous zones, shown in Figure A1. An initial hydraulic conductivity value has been assigned at the 4 zones, relying on the CSM values, then 317 PP have been inserted with a regular grid assigning an initial value equal to the one of the homogeneous to which they belong, while 20 PP are relative to pumping tests found in the AGISCO database (Regione Lombardia and ARPA Lombardia). In Table A1 are shown the upper and lower variability bound of the PP with respect of the homogeneous zone. Table A1. Pilot Points upper and lower variability bounds, relative to the homogeneous zone to which they belong.  Figure A1. (a) Pilot Point (green dots) and aquifer tests (dark and blue dots) used for calibration process with observation and validation head targets. The K-zone represents in brown the discontinuous aquitard whereas the light blue is the Aquifer A in layer 2, (b) cross-section in Figure A1b represents the vertical discretization of the numerical model: in blue the Aquifer A and in green the Aquifer B).