Next Article in Journal
A Modified Surface Energy Balance to Estimate Crop Transpiration and Soil Evaporation in Micro-Irrigated Orchards
Next Article in Special Issue
Prediction of Suspended Sediment Load Using Data-Driven Models
Previous Article in Journal
Scour Evolution Downstream of Submerged Weirs in Clear Water Scour Conditions
Previous Article in Special Issue
A New National Water Quality Model to Evaluate the Effectiveness of Catchment Management Measures in England
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Dipartimento di Ingegneria Civile e Ambientale (DICA), Politecnico di Milano, 20133 Milan, Italy
2
Tethys Srl, Viale Lombardia 11, 20131 Milan, Italy
*
Author to whom correspondence should be addressed.
Water 2019, 11(9), 1745; https://doi.org/10.3390/w11091745
Submission received: 1 August 2019 / Revised: 13 August 2019 / Accepted: 16 August 2019 / Published: 22 August 2019

Abstract

:
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.

Graphical Abstract

1. 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 G l / 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.

2. Materials and Methods

2.1. 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:
  • Steady flow simulation (1980s) calibrated using PEST [20] (abstraction wells and recharge parameters) in order to have initial condition for transient model.
  • 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.
  • 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.

2.2. 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.

2.3. Boundary Conditions

The boundary conditions (BCs) (shown in Figure 3) have been applied as follows:
(a)
Constant Head conditions (CH) at the model boundaries. In order to represent in a better way the transient conditions of the model, different CH values were chosen for each SP (time-variant CHs). Table 2 describes the documents used to extrapolate the CH values
(b)
Cauchy boundary conditions (RIV) for simulating the surface bodies Olona and Seveso and the Villoresi canal:
  • 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.

2.4. 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:
  • 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)
  • 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)).
Q c a l i b r a t e d , i = m i × Q i n i t i a l , i
where m i is the calibrated multiplication factor relative to the group to which the i-th well belongs, Q i n i t i a l , i is the initial withdraw value of the i-th well and Q c a l i b r a t e d , 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.
The 1980 (SP0) steady-state model (68 parameters: 5 vertical recharge zones and 63 multiplicative factors) was calibrated on a set of 92 observation, while the transient model (from 1989-SP1 to 2018-SP5) (340 parameters: 25 recharge zones and 315 multiplicative factors) was calibrated on a set of 458 head observations.

2.5. 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.

2.6. 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:
  • Constant concentration of source was set to 0;
  • 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.

3. Results and Discussion

3.1. 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.
The final value of the RMSE (Root Mean Square Error) stands at 2.90 , a value deemed acceptable considering the model complexity. Looking at the RMSE Scaled on the observation range (shown in 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.

3.2. 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.

3.3. 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).

3.4. Scenario in Three Different Time-Frame: 2021, 2027 and 2040

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.

3.5. 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.
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”.

4. 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.

Author Contributions

Conceptualization and methodology, all Authors contributed equally; software, G.F. and P.M.; writing—original draft preparation, L.C. and P.M.; writing—review, editing and supervision, L.A. and G.F.

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:
PAPilot Area
FUAFunctional Urban Area
CSMConceptual Site Model
BCsBoundary Conditions
ICsInternal Conditions
CHConstant Head
SPStress Period
SIFSistema Informativo Falda

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.
Table A1. Pilot Points upper and lower variability bounds, relative to the homogeneous zone to which they belong.
ZoneUnitLower Bound (m/s)Upper Bound (m/s)
Zone 1Aquifer A 10 5 10 3 10 2
Zone 2Aquitard 10 8 10 6
Zone 3Aquifer B 10 5 10 3 10 2
Zone 4Intermediate lenses 10 8 10 6
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).
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).
Water 11 01745 g0a1

Appendix B. Further Results of Flow-Transport Model

Figure A2. Scatter plot observed vs. simulated head relative to the overall model.
Figure A2. Scatter plot observed vs. simulated head relative to the overall model.
Water 11 01745 g0a2
Figure A3. Computed heads for (a) SP0 (1989), (b) SP2 (2004), and (c) SP5 (2018).
Figure A3. Computed heads for (a) SP0 (1989), (b) SP2 (2004), and (c) SP5 (2018).
Water 11 01745 g0a3
Figure A4. Comparison between measured and computed concentrations (B-T curves) for Monitoring Wells (a) Monitor_07, (b) Monitor_99, and (c) Monitor_19. The computed value is given by the average value simulated in the layer actually screened by the well.
Figure A4. Comparison between measured and computed concentrations (B-T curves) for Monitoring Wells (a) Monitor_07, (b) Monitor_99, and (c) Monitor_19. The computed value is given by the average value simulated in the layer actually screened by the well.
Water 11 01745 g0a4

References

  1. Alberti, L.; Azzellino, A.; Colombo, L.; Lombi, S. Use of cluster analysis to identify tetrachloroethylene pollution hotspots for the transport numerical model implementation in urban functional area of Milan, Italy. Int. Multidiscip. Sci. GeoConf. SGEM Surv. Geol. Min. Ecol. Manag. 2016, 1, 723–729. [Google Scholar]
  2. Werner, B. European Waters: Current Status and Future Challenges: Synthesis; Publications of the European Union; European Environment Agency: Copenhagen, Denmark, 2012. [Google Scholar]
  3. Civita, M.; Dal Pr, A.; Francani, V. Proposta di classificazione e mappatura della qualità delle acque sotterranee. Inquinamento. 1993, pp. 8–17. Available online: http://hdl.handle.net/11311/525710 (accessed on 22 August 2019).
  4. Environmental Ministry of Italy. Legislative Decree No. 152 Approving the Code on the Environment; FAO: Rome, Italy, 2006.
  5. Alberti, L.; Francani, V. Studio idrogeologico sulle cause del sollevamento della falda nell’area milanese. GEAM Geoingegneria Ambientale e Mineraria 2001, 104, 257–264. [Google Scholar]
  6. ARPA Lombardia. PROGETTO PLUMES - Integrazione. Available online: http://www.regione.lombardia. it/wps/portal/istituzionale/HP/DettaglioRedazionale/istituzione/direzioni-generali/direzione-generaleambiente-e-clima/piano-per-inquinamento-diuso (accessed on 5 June 2019).
  7. Stefania, G.; Rotiroti, M.; Fumagalli, L.; Zanotti, C.; Bonomi, T. Numerical Modeling of Remediation Scenarios of a Groundwater Cr (VI) Plume in an Alpine Valley Aquifer. Geosciences 2018, 8, 209. [Google Scholar] [CrossRef]
  8. Alberti, L. Model calibration using the automatic parameter estimation procedure (PEST) of the North-eastern zone of the Milan Functional Urban Area (Italy). Acque Sotter. 2018, 7, 27–38. [Google Scholar] [CrossRef]
  9. Pedretti, D.; Masetti, M.; Beretta, G.P.; Vitiello, M. A revised conceptual model to reproduce the distribution of chlorinated solvents in the Rho Aquifer (Italy). Groundw. Monit. Remediat. 2013, 33, 69–77. [Google Scholar] [CrossRef]
  10. Harbaugh, A.W.; Banta, E.R.; Hill, M.C.; McDonald, M.G. MODFLOW-2000, The U. S. Geological Survey Modular Ground-Water Model-User Guide to Modularization Concepts and the Ground-Water Flow Process; Open-file Report; U. S. Geological Survey: Reston, VA, USA, 2000; p. 134.
  11. Zheng, C.; Wang, P.P. MT3DMS: A Modular Three-Dimensional Multispecies Transport Model for Simulation of Advection, Dispersion, and Chemical Reactions of Contaminants in Groundwater Systems; Documentation and User’s Guide; Technical Report; Alabama University: Tuscaloosa, AL, USA, 1999. [Google Scholar]
  12. Bonomi, T. Database development and 3D modeling of textural variations in heterogeneous, unconsolidated aquifer media: Application to the Milan plain. Comput. Geosci. 2009, 35, 134–145. [Google Scholar] [CrossRef]
  13. Beretta, G.P.; Civita, M.; Francani, V. Idrogeologia per il Disinquinamento Delle Acque Sotterranee: Tecniche per lo Studio e la Progettazione Degli Interventi di Prevenzione, Controllo, Bonifica e Recupero; Pitagora Editrice Bologna: Milano, Italy, 1992. [Google Scholar]
  14. Alberti, L.; Cantone, M.; Colombo, L.; Lombi, S.; Alessandra, P. Numerical modeling of regional groundwater flow in the AddaTicino Basin: Advances and new results. Rendiconti online della Società Geologica Italiana 2016, 41, 10–13. [Google Scholar] [CrossRef]
  15. Cavallin, A.; Francani, V.; Mazzarella, S. Studio idrogeologico della pianura compresa fra Adda e Ticino. Costruzioni 1983, 326–327, 39. [Google Scholar]
  16. Department of Earth and Environmental Sciences, University of Milan-Bicocca. Tangram, Database per Pozzi [Well Database]. 1993. Available online: http://www.tangram.samit.unimib.it/ (accessed on 1 July 2019).
  17. Leapfrog. Available online: https://www.leapfrog3d.com (accessed on 1 July 2019).
  18. Regione Lombardia; ENI Divisione Agip. Geologia Degli Acquiferi Padani Della Regione Lombardia; ELCA: Florence, Italy, 2002.
  19. Colombo, L.; Gattinoni, P.; Scesi, L. Influence of underground structures and infrastructures on the groundwater level in the urban area of Milan, Italy. Int. J. Sustain. Dev. Plan. 2017, 12, 176–184. [Google Scholar] [CrossRef]
  20. Doherty, J.; Brebber, L.; Whyte, P. PEST: Model-independent parameter estimation. Watermark Comput. Corinda Aust. 1994, 122, 336. [Google Scholar]
  21. Thornthwaite, C.W. An approach toward a rational classification of climate. Geogr. Rev. 1948, 38, 55–94. [Google Scholar] [CrossRef]
  22. Alberti, L.; Alimi, H.; Ertel, T.; Trefiletti, P.; Pietrini, I. Fingerprinting and groundwater model application to evaluate hydraulic barrier efficiency (Italy). Environ. Forensics 2015, 16, 217–230. [Google Scholar] [CrossRef]
  23. Provincia di Milano (Assessorato Ambiente). Prelievi Idrici da Falda nel Milanese 1989; Sistema Informativo Falda: Milano, Italy, 1989.
  24. Provincia di Milano (Assessorato Ambiente). Prelievi Idrici da Falda nel Milanese 1990; Sistema Informativo Falda: Milano, Italy, 1990.
  25. Raffaelli, L.; Raimondi, P.; Rord, G.; Provincia di Milano. Valutazioni Sull’innalzamento Della Falda Nella Città di Milano nei Primi Anni ’90; Provincia di MilanoAssessorato All’ Ambiente: Milan, Italy, 1996.
  26. Lombardia-Eupolis, R. Piano di Tutela delle Acque. 2016. Available online: https://www.regione.lombardia.it/wps/portal/istituzionale/HP/DettaglioRedazionale/servizi-e-informazioni/Enti-e-Operatori/territorio/governo-delle-acque/piano-tutela-acque-pta-2016/piano-tutela-acque-pta-2016 (accessed on 1 July 2019).
  27. Doherty, J. Calibration and Uncertainty Analysis for Complex Environmental Models; Watermark Numerical Computing: Brisbane, Australia, 2015. [Google Scholar]
  28. Anderson, M.P.; Woessner, W.W.; Hunt, R.J. Applied Groundwater Modeling: Simulation of Flow and Advective Transport; Academic Press: Cambridge, MA, USA, 2015. [Google Scholar]
  29. ASTM D5981-96(2008), Standard Guide for Calibrating a Groundwater Flow Model Application (Withdrawn 2017); ASTM International: West Conshohocken, PA, USA, 2008; Available online: www.astm.org (accessed on 1 July 2019).
  30. Wels, C.; Mackie, D.; Scibek, J. Guidelines for Groundwater Modelling to Asses Impacts of Proposed Natural Resource Development Activities; Technical report; Report n°194001; British Columbia Minstry of Environment—Water Protection and Sustainability Branch: Victoria, BC, Canada, 2012.
  31. Alberti, L.; Colombo, L.; Francani, V. The groundwater flow velocity distribution in the urban areas: A case study. Ital. J. Eng. Geol. Environ. 2014, 14, 17–26. [Google Scholar]
Figure 1. (a) Functional Urban Area (FUA) within the Lombardy Region and Milano Province, (b) FUA and pilot area (Municipalities involved in the AMIIGA project are colored by violet whereas the Milano City is line filled. Hydrogeological cross sections are represented by lines).
Figure 1. (a) Functional Urban Area (FUA) within the Lombardy Region and Milano Province, (b) FUA and pilot area (Municipalities involved in the AMIIGA project are colored by violet whereas the Milano City is line filled. Hydrogeological cross sections are represented by lines).
Water 11 01745 g001
Figure 2. Hydrogeological cross-section (a) N-S (blue line) and (b) E-W (red line) of the Pilot Area represented in Figure 1b. The light-blue represents the Aquifer A whereas the green represents the Aquifer B.
Figure 2. Hydrogeological cross-section (a) N-S (blue line) and (b) E-W (red line) of the Pilot Area represented in Figure 1b. The light-blue represents the Aquifer A whereas the green represents the Aquifer B.
Water 11 01745 g002aWater 11 01745 g002b
Figure 3. Boundary Conditions and Internal Conditions: Constant head (blue), rivers and canal (light-blue), and withdrawals differentiated by exploited Aquifer. The transport model is developed within the area of interest.
Figure 3. Boundary Conditions and Internal Conditions: Constant head (blue), rivers and canal (light-blue), and withdrawals differentiated by exploited Aquifer. The transport model is developed within the area of interest.
Water 11 01745 g003
Figure 4. Flow chart of WMG program.
Figure 4. Flow chart of WMG program.
Water 11 01745 g004
Figure 5. Scaled Root Mean Square Error (RMSE) of each stress period (SP).
Figure 5. Scaled Root Mean Square Error (RMSE) of each stress period (SP).
Water 11 01745 g005
Figure 6. Monitoring wells and Monitoring hydrographs location in the model domain.
Figure 6. Monitoring wells and Monitoring hydrographs location in the model domain.
Water 11 01745 g006
Figure 7. Comparison between observed and calculated hydrographs (presented in Figure 6 for some monitoring points in the domain: (a) FOG57, (b) FOG56, (e) Pz5 are located within the area of interest while (c) Well_3A, (d) Pz265, and (f) Well_3B are outside the area. A and B in brackets indicate the screened aquifer.
Figure 7. Comparison between observed and calculated hydrographs (presented in Figure 6 for some monitoring points in the domain: (a) FOG57, (b) FOG56, (e) Pz5 are located within the area of interest while (c) Well_3A, (d) Pz265, and (f) Well_3B are outside the area. A and B in brackets indicate the screened aquifer.
Water 11 01745 g007
Figure 8. Comparisons between transport simulations in the Aquifer B carried out in ARPA Plume 2015 and AMIIGA 2018 (average concentration of layers from 3 to 9). Figures (a,b) represent 1989 (35 years of simulation) while figures (c,d) represent 2014.
Figure 8. Comparisons between transport simulations in the Aquifer B carried out in ARPA Plume 2015 and AMIIGA 2018 (average concentration of layers from 3 to 9). Figures (a,b) represent 1989 (35 years of simulation) while figures (c,d) represent 2014.
Water 11 01745 g008aWater 11 01745 g008b
Figure 9. The rows represent the simulated plume for the same layer while the columns represent the same year (2021, 2027, and 2040, respectively). Figures (ac) are relative to layer 1 (Aquifer A, shallow), while figures (df) represent layer 7 (Aquifer B, deep).
Figure 9. The rows represent the simulated plume for the same layer while the columns represent the same year (2021, 2027, and 2040, respectively). Figures (ac) are relative to layer 1 (Aquifer A, shallow), while figures (df) represent layer 7 (Aquifer B, deep).
Water 11 01745 g009
Figure 10. Breakthrough curves of the contaminant registered in (a) Monitoring Well 07, (b) Monitoring Well 99, and (c) Monitoring Well 19.
Figure 10. Breakthrough curves of the contaminant registered in (a) Monitoring Well 07, (b) Monitoring Well 99, and (c) Monitoring Well 19.
Water 11 01745 g010
Figure 11. Proposed monitoring networks for (a) Aquifer A and (b) Aquifer B.
Figure 11. Proposed monitoring networks for (a) Aquifer A and (b) Aquifer B.
Water 11 01745 g011
Table 1. Stress Period definition.
Table 1. Stress Period definition.
PeriodStress PeriodDuration (year)Time Step (-)
1980–1989Steady-state/SP 010-
1990–1999SP 11010
2000–2004SP 2510
2005–2009SP 3510
2010–2014SP 4510
2015–2018SP 5410
Table 2. Constant Head sources.
Table 2. Constant Head sources.
Stress PeriodPiezometric MapSourceReliability
SP 0Regional piezometry map 1982Regione LombardiaLow
SP 1Regional piezometry map 1996Regione LombardiaMedium
SP 2Regional piezometry map 2003Regione LombardiaMedium
SP 3Regional piezometry map 2007Regione LombardiaMedium
SP 4Regional piezometry map 2014Area Vasta model [26]Medium
SP 5Piezometric campaign (March 2018)AMIIGA ProjectMedium
Table 3. Upper and lower variability bounds for the two well macro-groups defined.
Table 3. Upper and lower variability bounds for the two well macro-groups defined.
Macro-GroupGroup NumberLower BoundUpper Bound
Municipalities10–560.52
Milano-Pumping station57–720.991.01
Table 4. Time-frame Scenario.
Table 4. Time-frame Scenario.
Stress PeriodTimeTime-Frame
SP 620211-Short Term
SP 720272-Middle Term
SP 820403-Long Term
Table 5. Statistics relative to each simulated aquifer.
Table 5. Statistics relative to each simulated aquifer.
ValuesABTotal
Absolute Residual Mean1.972.362.24
Sum of Squared Residuals101928323852
Observation number150308458
RMS2.613.032.90
Scaled Absolute Residual Mean0.030.020.02
Scaled RMS Error0.040.030.03
Table 6. Total water budget (inflow-outflow in m 3 /s) for each stress period (SP).
Table 6. Total water budget (inflow-outflow in m 3 /s) for each stress period (SP).
ComponentSP1SP2SP3SP4SP5
InOutInOutInOutInOutInOut
Storage0.050.090.990.210.011.510.370.110.180.08
CH14.819.727.567.5114.0711.9711.3510.2111.5610.18
Wells0.0010.350.005.870.005.950.007.610.007.43
River leakage0.230.000.230.000.230.000.230.000.230.00
Recharge5.260.005.260.005.230.005.460.005.740.00
Total (m 3 /s)20.3520.1614.0413.5919.5419.4317.4117.9217.7117.69
Table 7. Differences between ARPA Plumes 2015 and AMIIGA 2018 models.
Table 7. Differences between ARPA Plumes 2015 and AMIIGA 2018 models.
FeatureARPA Plumes 2015AMIIGA 2018
Domain and grid100 × 100 km 2 , 20 × 20 m cells 21 × 18.8 km 2 , 100 × 100 m cells
Vertical discretization (layers)39
Modeled sources21
BCs (Constant Heads)ConstantTime-varying
ICs (Wells abstraction, Recharge)ConstantTime-varying

Share and Cite

MDPI and ACS Style

Colombo, L.; Alberti, L.; Mazzon, P.; Formentin, G. Transient Flow and Transport Modelling of an Historical CHC Source in North-West Milano. Water 2019, 11, 1745. https://doi.org/10.3390/w11091745

AMA Style

Colombo L, Alberti L, Mazzon P, Formentin G. Transient Flow and Transport Modelling of an Historical CHC Source in North-West Milano. Water. 2019; 11(9):1745. https://doi.org/10.3390/w11091745

Chicago/Turabian Style

Colombo, Loris, Luca Alberti, Pietro Mazzon, and Giovanni Formentin. 2019. "Transient Flow and Transport Modelling of an Historical CHC Source in North-West Milano" Water 11, no. 9: 1745. https://doi.org/10.3390/w11091745

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop