Stochastic Particle Tracking Application in Different Urban Areas in Central Europe: The Milano (IT) and Jaworzno (PL) Case Study to Secure the Drinking Water Resources

Abstract

Urban areas are typically characterized by the presence of industrial sites, which are often sources of groundwater contamination, posing a serious threat for the groundwater. In such cases, a crucial step is to find a link between the contaminant sources and freshwater supply wells at risk. As a part of the AMIIGA Project, two different stochastic approaches were applied to assess drinking water supply wells vulnerability in Functional Urban Areas in the presence of several chlorinated hydrocarbons sources in an alluvial aquifer in Milano and a pesticide mega site in a complex geological setting in Poland. In the first case study, the innovative Pilot Point Null-Space Monte Carlo forward particle tracking was used, applying a forward solution instead of the classical backtracking, while in the second case was chosen the classical Monte Carlo methodology. Both case studies represent useful application examples, allowing an effective prioritization of expensive remediation actions in order to protect freshwater wells.

1. Introduction

The problem of groundwater resource quality degradation in highly urbanized areas is one of the most important environmental issues at both the European [1,2] and national levels (Legislative Decree 152/2006 for Italy and Dz.U. 2017 poz. 1566 (Water Law), DzU z 2006, nr 123, poz. 858, DzU z 2006, nr 129 poz. 902 (Environmental Protection Law) for Poland) during the last three decades. Recently, national and regional regulations are considering the necessity to develop plans for the remediation and management of the most industrialized urban areas affected by groundwater contamination. Functional Urban Areas (FUA, [3]), i.e., large areas composed by an urban core and its commuting zone, compose a unique contamination core where the groundwater is affected by many different kinds of sources, and therefore are difficult to manage. In general, the causes of groundwater contamination can be classified in three different classes [4,5]: (a) Point Sources (PS), hot spots of contamination corresponding to areas releasing plumes of high concentration [6]; (b) Multiple Point Sources (MPS) which are unknown and for which the problem can be faced by using a multi-methodological approach as developed in [7,8,9]; (c) Non-Point Sources (NPS), where the contaminant load comes from the development of anthropogenic activities on large areas (for example pesticides from agricultural practices [10]). In urban contexts, PS contamination (i.e., brownfield, old chemical plant, refineries) are more frequent than NPS (i.e., leakages from agricultural areas [5]). PSs are usually located inside a contaminated site that represents active or abandoned industrial area. Sometimes, some of these PSs generate sewage or effluents that were discharged via surface, lagoons, tanks, pipes or in sewers, near or within well protection areas. In order to protect the groundwater resources, in particular water supply wells, it would be necessary to identify the main sources responsible to the contamination and to track the main contaminant movements. This could help Public Authorities to activate emergency safety plans and reclamation procedures avoiding high costs of water clean-up, which are eventually charged by Water Managers to citizens (e.g., active carbon filters and associated management).

As the FUA is composed of a mixture of urban, agricultural and industrial areas, two different kinds of PSs can be distinguished:

-

Large mega sites, which are located in a former industrial and often abandoned sub-region, historically known (sometimes in operation for more than half a century) and whose responsibility to groundwater contamination is well recognized. There are many cases in the literature [11,12,13] where groundwater contamination by chlorinated hydrocarbons (CHCs) or pesticides is a major cause of deteriorating drinking water quality. At a local scale, a mega site can be associated to several PSs but at a regional scale it should be considered as a unique large contaminated area.

-

Several point sources, which are corresponding to several industrial sites whose contribution has been discontinuous during time due to different industrial activity (i.e., change of final product, change of use of product). In this case, it is difficult to find the liable part as the contamination spreading from these sites can create a superimposition of different plumes [14,15].

These two kinds of contamination sources have been taken into consideration within the Interreg Central Europe project (AMIIGA–CE032, 2016–2019). The aim of the project was to address the problem of groundwater contamination in European FUAs (Figure 1), in order to identify the areas most affected by CHCs contamination. For this reason, the main objective of the present paper is to develop stochastic approaches in heavily industrialized areas in order to recognise which water supply wells are directly involved by contamination and secure those wells before the possible arrival of contamination. The paper presents two application examples, one located in the Milano FUA (northern Italy, Site 1) and one in Jaworzno (southern Poland, Site 2) (Figure 2). Such areas have a high relevance in central Europe, due to their historical contaminations. Milano historically hosts many industries and chemical plants located upgradient of the city pumping stations, posing a relevant threat for the freshwater supply while Jaworzno hosts a Mega Site that produced pesticides and contaminated the groundwater. The Site 1 is characterized by many potential sources (brownfields and industrial sites) and potential receptors (pumping stations) in a relatively simple hydrogeological setting, while the Site 2 presents only one and well-defined source (pesticides mega site) and a unique supply well, in a complex and uncertain geological context. Therefore, two different methodologies were applied. In the first case, a Null-Space Monte Carlo (NSMC) coupled with forward particle tracking (PT) method was chosen. Differently from previous works that were focused on a backward particle tracking analysis to search the potential sources [4,16,17,18], here the forward tracking is used to find the sensible receptors that are probably impacted by the contamination (i.e., water supply wells). In the second case, a classical Monte Carlo (MC) coupled with PT was preferred, due to the most complex hydrogeological conditions and the lack of data availability.

This procedure can be useful to Public Authorities to avoid huge costs of management and sanitation that often are charged to private citizens due to the not-individuation of the contamination responsible.

2. Materials and Methods

2.1. Methodology Description: A Decisional Framework to Secure Water Supply Wells

The methodology developed during the AMIIGA Project is schematized in Figure 2. In order to identify the water supply wells more affected by contamination, two different stochastic methods were applied in two areas located in Central Europe. The first location is the northern sector of Milano (Italy, Lombardy Region), characterized by the presence of many industrial sites producing/using CHCs that affected the quality of groundwater; this first area is identified as Site 1. The second area is located in the southern sector of Poland (Jaworzno Region), where a pesticide mega site has been contaminating groundwater resources since many years; this site is called Site 2. A groundwater flow model developed in MODFLOW (MODFLOW2005, [19]) was implemented, calibrated and finally used to perform a deterministic particle tracking simulation starting from a known contamination source (MODPATH, [20]), in order to obtain a preliminary identification of water supply wells potentially affected by the contamination. Since groundwater models contain many intrinsic uncertainties (i.e., the direction of flow and so the path of contaminants [17,21,22]), a stochastic approach has been developed for both areas considering the uncertainties linked to the hydraulic conductivity (main parameter governing the flow), and so the main flow direction. The main steps are summed up in Figure 2 and explained below:

Groundwater models using the MODFLOW code [19] are extensively used for different applications [23,24,25,26,27,28]. Basically, the model is able to reproduce the water table and the main direction of the flow taking into account the conceptual model of the studied area (aquifer geometry, aquiclude extension and hydrogeological parameter reconstruction, i.e., hydraulic conductivity and porosity). Moreover, the water balance of the model can help to understand the main variables involved in the hydrogeological balance (i.e., recharge, boundary conditions, abstraction wells) within the study area. Then, by means of MODPATH [20], it is possible to track the contaminant under the hypothesis of pure advective flow. The code calculates the forward and backward particles path, identifying possible sensitive receptors downgradient of the source or the more probable source areas.

Stochastic Particle Tracking—Null-Space Monte Carlo (PT-NSMC): for Site 1 (IT), deterministic modelling (i.e., transport model) is a widespread and useful tool to investigate contaminant plumes and related sources; nevertheless, when it has to deal with areas affected by several plumes, the deterministic approach fails because the modeler is unable to study in a deterministic way the influence of many contaminated sources to the water supply wells. Stochastic modeling is the most powerful tool for source apportionment. One of the most used methods is combining PEST code [29] with the NSMC method [4,17,30,31,32]. It is possible to provide a mechanism for the exploration of the uncertainty associated to measurement errors, enabling the model to fill the observed data with a suitable stochastic descriptor of parameter variability within the FUA. The combination of forward PT and stochastic simulations is employed to narrow the possible area affected by the contamination detected in water supply wells resulting from the considered sources. In the present work, the stochastic forward technique is applied instead of the backtracking one. This choice allowed to find the water supply wells that are evenly more affected by the contamination coming from several PSs. For Site 2, due to the presence of strong uncertainties in the conceptual model, the multiple simulations approach has been applied in order to indicate the areas which should be secured from the contamination.

In this paper, PT+MC was used for both sites but in different mode: in Site 1, the NSMC technique has been preferred, due to the fact that in the study area there is availability of a lot of aquifer tests, many piezometer data (i.e., head observation targets), many hydrogeological logs used to reconstruct the main aquifers and several contaminants concentration data. In Site 2, due to limited information concerning the conductivity of rocks downstream the contaminated site, over 60 model variants were performed sampling the values in an interval centered within the hydraulic conductivity of the native model (classical MC technique).

2.2. Site 1 (S1)—Milano FUA: A PT-MC Approach to Assess the Area Involved by Several Point Sources

2.2.1. Study Area, Hydrogeology and Conceptual Model

The study area (named S1) is located in the Lombardy Region (Figure 3), within the Po Plain, and includes the N-W sector of Milano and some surrounding municipalities (Pilot Area in the AMIIGA Project, PA). The total area is 395 km² and lies at the center of one of the most urbanized and industrialized areas in Europe. Such area is widely affected by heavy contamination caused either by some known plumes [26] and by diffuse pollution [7]. In the north-western part of Milano Province, there are many suspected contaminated sites (as automotive, refineries, chemical plant, chromium plating, galvanic) which had been or are still active during most of the Twentieth century. Some of these sites have been subjected to site characterization and reclamation since the 1990’s, hence a significant number of stratigraphic, hydrogeological and concentration data are available covering a period of about thirty years.

The main aquifers in the Lombardy Plain are hosted by a sequence of plio-pleistocenic sediments that filled the Neogene Po plain foredeep reaching a maximum thickness of approximately 500 m [33,34]. At the base of the sequence (Pliocenic sediments), pelagic sediments are mainly composed by clay and silt. At the top, gravel and sand are predominant and fed the basin under the control of middle to late Pleistocene glacial cycles. The origin of the latter is alluvial and glacio-fluvial and the hosted aquifers reveal high transmissivity values. Four main aquifers (named from A to D) are recognized in the Lombardy Plain but only 3 have been well investigated through drillings in the Milano area [35,36,37,38]. The case study discussed in this paper focuses only on the aquifers A and B, which are separated by a clay layer in the south part of the domain (the aquitard is represented by a brown line in (Figure 3), that becomes discontinuous and then disappears northward. Based on the previous studies [18,26], the hydrogeological conceptual model has been updated with the most recent stratigraphic data retrieved from local authorities (Environmental Protection Agency (ARPA) and Lombardy Region).

The groundwater flow model used in this paper comes from previous works [18,26], however some essential details are reported in the Supplemental Material (Table S3). Although the NSMC method was already presented, the main aim herein is to identify the main receptors (i.e., public wells) starting from several sources formerly identified by the use of cluster analyses [7,39].

The model was previously implemented by reconstructing the layer characteristics (thicknesses and hydraulic properties) and the low-permeability lens extensions using many stratigraphic logs available for the area [18,26]. The first model layer has an average thickness of 30 m and represents the shallow aquifer (Aquifer A), the most impacted by contamination nearby the suspected sources. The hydraulic conductivity ranges between 1 × 10⁻⁵ and 4 × 10⁻³ m/s. The aquitard (layer 2 and layer 3) is discontinuous in the studied area and it is located at a depth of about 40 m with a thickness of 5 m with a hydraulic conductivity ranging from 1 × 10⁻⁹ m/s and 1 × 10⁻⁷ m/s. Towards north the aquitard disappears (Figure 4), allowing the hydraulic connection between the shallow and semi-confined aquifer. Underneath, a second aquifer is present (layer from 3 to 9 is the Aquifer B with a hydraulic conductivity varying between 1 × 10⁻⁵ and 1.4 × 10⁻³ m/s.), with a thickness of about 60 m, that is considered semi-confined as includes discontinuous clay lenses that only locally confine the system.

The Boundary Conditions are selected to be Constant Head (CH) and are considered suitable as the area impacted by the contamination is far away from such boundaries, being located in the Milano City area. The groundwater flow model was calibrated with PEST (using the Pilot Point (PP) technique [29]) considering 63 head targets relative to a piezometric survey executed in 2018. The absolute residual mean is about 0.74 m. In the present study, the hydraulic conductivity (K) was considered as the unique uncertain parameter because it is assumed to be the parameter most affecting the main flow direction. Therefore, in the following steps, the K is considered as a stochastic parameter (i.e., varying parameter in the NSMC procedure).

2.2.2. Application of the NSMC Method

The first step to apply the NSMC technique [40] is to calibrate the flow model (previously conducted in [26]. Secondly, initial K-values have to be attributed in each PP (i.e., for each layer) by kriging interpolation using the calibrated K-values at each PP position. An appropriate regularization is employed to obtain the minimum error variance solution to the non-unique inverse problem and to avoid the emergence of unrealistic parameter values [41]. Since the parameter that is assumed to most influence the flow direction of the path is the hydraulic conductivity, it is the main source of predictive uncertainty. The resulting uncertainty also affects the calibrated model that is used in the first step of the NSMC procedure and, for this reason, a prior covariance matrix centered on the calibrated parameter field is used to generate random parameter fields. To evaluate advective particle variability induced by the K-field distribution, random realizations of the model parameter values are generated based on a stochastic characterization of their variability and spatial correlation.

Procedure Overview

Once the model was built and calibrated, the NSMC analysis has been applied following the procedure widely presented and discussed in the literature [17,18,31]. A Flow diagram is presented as Supplementary Material (Figure S1).

The NSMC method is an efficient procedure able to produce a multitude of calibration-constrained stochastic parameter fields (see the flow chart in Figure S1, adopted in the AMIIGA Project and applied in [18]). The steps are already presented in previous works and here, only a short list is presented:

Step 1: Generation of a multitude (150) stochastic fields (i.e., K-field). PPs are attributed to the model based on a 2 km × 2 km uniform grid. Additional PPs are positioned near the aquifer test locations available in the contaminated sites. All 337 PPs are assigned to a single regularization group. A very wide range of K-values from 10⁻⁵ to 10⁻² m/s was considered for the more permeable zones while an interval from 10⁻⁹ to 10⁻⁷ m/s was chosen for clay lenses, in order to avoid any constraint in the generation of random values. Under the assumption of normality, the bounds for each parameter are assumed to lay within the 95% confidence interval. Hydraulic conductivity is considered to be isotropic in the horizontal plane (Kx = Ky) while the Kz is assumed equal to 0.1 the horizontal permeability. The spatial structure of the K-values is described by an exponential variogram.

Step 2: Each stochastic field is projected into the calibration null-space; therefore, the solution space component of this field is removed and substituted with the calibrated parameter field. The calibration dataset is composed of hydraulic head observations (63 head values) representing the measurements of the 2018 dataset.

Step 3: Because of the non-linearity of the problem solution, the null-space projected stochastic field requires re-calibration with an iterative procedure.

Step 4: Using MODPATH, it is possible to perform a forward tracking (considering only advective transport) of the particles positioned nearby the monitoring wells/piezometers suspected to be a potential contamination source (Figure 4). The starting locations are placed in industrial sites where high CHCs concentrations have been observed to be between 10 and 20 times the threshold value within the historical time series (e.g., the limit for Tetrachloroethylene in the Italian law is 1.1 µg/L). A MODPATH simulation is undertaken for each of the 150 randomly generated models.

At the end of the cycle, the total number of particles that passed through each model cell is summed for all the simulations and divided by the number of realizations (150). The result is a frequency map for all layers that highlights the area most affected by the passage of particles, i.e., by contaminated groundwater.

Stochastic PT Results

The stochastic PT frequency map is shown in Figure 4. The frequency of passage has been computed normalizing the number of particles for each cell with respect to the maximum passages in all the model cells. The map shows areas potentially affected by the contamination coming from the several suspected sources located upgradient of the contaminated wells. It must be highlighted that also in those areas where the frequency is low, the possibility that the source area might impact the groundwater quality is still not negligible [18]. For example, the yellow areas have the probability to be crossed by particle contamination with a frequency varying from 3% to 26% whereas the percentage in the green areas is very low (1–3%), but still it cannot be excluded that the groundwater is affected by the contamination. The passage frequency computed for each model cell is higher close to the sources and progressively decreases downgradient due to the different paths that the contaminant can follow due to the different sampled K-field. The resulting map generally shows a complex situation: in fact, especially near the Milano borders, there is a high probability of the pathlines coming from different sources to superimpose (i.e., plume coalescence). It is also clear how all the western most pumping stations are likely to be affected by such contamination.

The results analysis can be divided in 3 different zones:

(1)

Northern part of Milano: there are 2 suspected sources. Source n°S10 heavily impacts the groundwater within 3 km from it (yellow area), then the impact decreases (with lower, but non-zero probability) but can be still tracked far as 7 km downgradient of the area. On the contrary, source n°S9 has a lower impact and the probability to affect the groundwater downgradient decreases quickly with distance (less than 1% within 3 km).

(2)

Municipalities near Milano: the presence of a large industrialized area immediately in the outer border of Milano involved many suspected sources (from n°S3 to n°S8). Sources n°S5 and n°S8 seem to not affect the water supply wells in Milano, as the yellow area (frequency of passage is about 10%, representing the main direction of particles) flows outside the city. The PT coming from sources n°S4 and S7 are superimposed: the probability that contamination coming from them can affect Novara, San Siro and Tonezza pumping stations is not negligible. The main direction of the plume from source n°S3 is towards Cimabue and Chiusabella pumping stations, with a frequency close to 30%. The last source (n°S6) causes a contamination of the very close Vialba pumping station (the yellow area has an extension no more than 2 km from the source), even if the dismissed Espinasse pumping station, closed in the 80’s for the highly contaminated waters, is also interested by the particle arrival.

(3)

Milano municipality: some brownfields are located inside the city borders. The source n°S11 highly affects Cimabue pumping station, being just upgradient of this well field. The high frequency of passage (around 30%, orange area) is due to the superimposition with source n°S3, located immediately upgradient outside the city borders. Source n°S1 and S2 may affect both the Cimabue and the dismissed Espinasse Stations with a relatively high frequency (10–20%).

2.3. Site 2 (S2)—Jaworzno FUA: A Deterministic PT Approach to Assess the Area Involved by a “Mega Site” Source

2.3.1. Study Area, Hydrogeology and Conceptual Model

Jaworzno is an industrial city in southern Poland (Figure 1), which is an area of major concentration of industry in the country and one of the most industrialized areas of Europe. The major contamination source (i.e., mega site) in Jaworzno is the chemical industry plant Organika-AZOT together with the nearby areas where its waste was disposed, mainly from the 1960’s till 1980’s (named S2, Figure 5). This pesticide waste collection area is widely known and well identified, and it is even listed as a serious hot-spot posing a threat to the distant Baltic Sea, few hundreds of kilometers to the north. However, the threat that this PS may pose to groundwater intakes that operate within the city of Jaworzno and its FUA have never been evaluated in detail. From a geographical point of view, there are three different areas included into the model domain:

(1)

Quaternary valley of Wąwolnica river, starting from the Chemical Plant Organika-AZOT and further downstream to the Przemsza river valley

(2)

Triassic hills located south of the Wąwolnica valley

(3)

Quaternary valley of Przemsza river, located above and below its confluence with Wąwolnica.

The main rock formations in the modelled area are:

•

Triassic limestones and dolomites of middle Triassic of triple porosity fissured and karst character [42] (in Figure 6b is the green color).

•

Neogen—glacifluvial sands and gravels forming Quaternary aquifers interbedded with discontinuous impermeable clay layers–variable thickness, in post-glacial buried valleys up to a few tens of meters (in Figure 6b is represented by various zones of K-field marked by all other colors except green).

More complexity is added by presence of few hundred meters thick carboniferous strata underneath the model domain, which are under active exploitation of mineral deposits resulting in intensive drainage of the groundwater aquifers [43].

The model presents a regular grid with 550 rows and 450 columns composed by 10 m × 10 m cells, for a total of 358,023 active cells. The vertical discretization consists of three layers. The model covers the Quaternary and partly the Triassic aquifer and the boundary between them runs along the fault zone at the southern bank of Wąwolnica river. Therefore, both Triassic (from layer 1 to layer 3 only in southern-eastern part, represented by green zone on Figure 6) and Quaternary formations are represented in each model layer. Concerning the Quaternary formation, it hosts two aquifers represented by the first and the third model layers (thickness: 10–15 m and 3–25 m respectively), while the layer 2 represents an aquitard (3–9 m thick) that separates the two aquifers (the yellow zone of Figure 6). The aquitard disappears towards the center of the Przemsza river valley in the West.

The eastern and western horizontal boundaries of model were modeled using Type 3 (Cauchy) conditions (GHB package in MODFLOW), since the head are linked to the regional groundwater flow. On the contrary, the northern and southern boundaries were modelled with Type 2 conditions (No-Flow). Most of groundwater recharge in the area of the model takes place mainly through rainwater infiltration. The drainage occurs mainly through the rivers Przemsza and Wąwolnica as well as with a system of drainage ditches in the vicinity of the landfills located in the former sand pit “Rudna Góra”. Another element of drainage is the vertical leakage of groundwater from Quaternary aquifer into the Carboniferous strata. The Carboniferous is drained to large extent by an active hard coal mine located underneath the Quaternary aquifer, which is also drained by the mine in particular areas of preferential connection between Quaternary and Carboniferous formations. This particular area is represented by the DRAIN package, shown in Figure 7. Another important element of drainage of the model is also a groundwater intake “Jarosław Dąbrowski” located in a closed mine shaft. This intake works in the Carboniferous aquifer (WELL package in Figure 7), but for similar reasons as the mine drainage described above, it affects the water balance of the Quaternary strata in the model domain; therefore, it was also included in the third layer of the model. The model calibration has been done in steady-state, using 16 head observation targets. The main parameters subjected to a trial-and-error procedure were hydraulic conductivity, recharge and hydraulic conductivity of riverbed deposits. The details are reported in the Supplemental Material (Table S1). The initial values of model parameters were taken from previous studies [44].

Anderson et al. [45] suggest the reference of scaled RMSE below 10% as the indication of a well-calibrated model, value respected by the model with a final scaled RMSE equal to 3% (Figure 4). Concerning the mass m² balance of the whole model, the main inflow comes from GHB (46%) and recharge from precipitation (32%), while main outflow is due to drainage by the rivers Przemsza and Wąwolnica (54%). Some further details are provided in Table S2 in Supplementary Material.

2.3.2. MC PT Procedure and Results

The PT simulation was performed using MODPATH [20]. The particles were placed at the plane running across the area contaminated with pesticides, located in the valley of the Wąwolnica river near the of Organika-AZOT chemical plant and its nearby waste collection areas (Figure 5). Then, the particles have been tracked forward in order to delineate the zone that is potentially affected by the contamination. The particle tracking simulations have been performed for:

•

Native model representing the calibrated flow model (Figure 8), where the k-hydraulic conductivity field fit the optimized head level (black paths in Figure 9b).

•

48 combinations of hydraulic conductivity values for the area of Przemsza buried river valley (pink and light blue zones in Figure 6), also producing converging model variants with acceptable model error (orange paths in Figure 9b).

•

12 combinations of hydraulic conductivity values for the area of Przemsza buried river valley (pink and light blue zones in Figure 6) for which model is not converging, but the nRMSE value is still at acceptable level—still below 10% (brown paths in Figure 9b).

•

The Figure 9a presents the normalized root mean square error (nRMSE) for all collected 61 model variants vs the K-field of Quaternary deposits in Przemsza valley downstream the contaminated site. The native model (i.e., the best calibrated model) conductivity value is represented by the black dot in Figure 9a.

Figure 9b shows all the possible pathlines obtained through the particle tracking analysis fields starting from the contaminated site. Different pathlines are given by different of K-fields and then, consider a different flow regime.

The modelling results indicate that in the best calibrated model, the groundwater contaminated with pesticides is not flowing directly into the nearby water well (black lines). However, considering the other 48 possible combinations of the parameters obtained by the MC PT method that are also producing converging model variants with acceptable model error, the contaminated groundwater is likely to flow very close to the water intake (orange lines). Moreover, including also another 12 parameter combinations (still obtained with the MC PT procedure) for which model is not converging, but with a model error still acceptable, for some extreme variants the contamination may even reach the pumping well of groundwater intake (brown lines). Therefore, it is vital to perform regular wellhead protection monitoring at the groundwater intake in order to be prepared in case of pesticide contamination approaches the pumping well used for public water supply.

3. Discussion and Conclusions

The current paper presents two different applications of the particle tracking methods applied to different investigated groundwater systems. In both cases multiple model variants have been analysed in order to explore the intrinsic uncertainty of the groundwater flow controlling parameters. In this case, the hydraulic conductivity is considered the property that most influences the flow direction, therefore it was considered the uncertain parameter. The methods result in a variety of possible solutions instead of only one deterministic best fit model, allowing to understand all the possible contaminant paths. Such methods can be an efficient alternative to solute transport models, because of the consideration of parameter uncertainty and because they require less parameters than a transport model; indeed, only a calibrated flow field is required while transport parameters (often hard to estimate) are not necessary. The two different approaches adapt to different contexts of geology, data availability, kind of source and number of potential sensible receptors.

In case of S1 (Milano, Italy) were considered many potential sources (active and dismissed industrial sites, chemical plants and brownfields) from which the contamination can flow through alluvial aquifers, with many sensible receptors located downgradient (i.e., public pumping stations). In this case, the stochastic PT simulation have been performed in order to calculate the particle passage frequency for each cell. This allowed the identification of the areas potentially affected by the contamination coming from the several suspected sources and the most probable responsible sources. In this approach, the more frequently the particles pass from a given model cell, the more likely it is that the area represented by this cell is affected by contamination. It is known that the water supply well stations in Milano, due to the high drainage induced by the high-water extraction demanding in the semi-confined aquifer, are more affected by the CHCs contamination coming from different sources. Moreover, the map represented in Figure 4 could be used by Public Authorities to assess the most suitable and useful intervention to early detect the contamination (i.e., monitoring networks) and to remediate the contamination if already arrived at the pumping stations (e.g., using pumping wells upstream the pumping stations as a hydraulic barrier). This can prevent the huge costs that Water Managers would have to cover in case the contamination enters on the water supply network. In this context, the outcomes of this method are certainly comparable to a transport model (including advection-dispersion) as the extension of non-zero probability areas is quite similar to the plume simulations as presented in [9].

In case of S2 (Jaworzno, Poland), only a single and well-known source is considered (mega site) in a complex hydrogeologic context, with a sensible receptor represented by a pumping well. The modelling results of best calibrated model indicate that the pesticide plume is not flowing directly into the nearby water well. Due to the complex geology and conceptual model of the area, the Monte Carlo simulation is a first step to assess the areas interested by the contaminant pathways, considering the uncertainties of the geological parameters. In this paper, it is the first time that a Gaussian sampling of K-field was assessed in Jaworzno FUA. Considering a large combination of K-field for which the model error is still at acceptable level, the threat for drinking water intake should be noted, as some of the particle pathlines reach the well (Figure 9b). Therefore, a wellhead protection monitoring program is suggested.

In future, relying on the current results, the level of protection should be more studied and determined also based on other uncertainties that influence the flow (i.e., pumping well, recharge and time-variant boundary conditions). Thus, possible improvements should focus on the implementation of tests that could significantly reduce the uncertainty of model predictions (e.g., tracer tests [46]), and consider uncertainties analyses for the other parameter [17]. However, the stochastic results can be still compared with solute transport models that also consider dispersion, diffusion and sorption, to achieve a better understanding of the contaminant movements.

In conclusion, this work reports two hydrogeologically different examples useful for water managers and local authorities to understand the advantages of MC methods in order to manage the water quality in urban areas.

Two different methodologies (NSMC + PT and MC + PT) were applied in two hydrogeological FUAs in Europe to map the uncertainties linked to flow path of the shallower aquifers. The obtained results can be very useful in the design of an effective monitoring system to prevent the pollution risk in both studied areas, even though two different approaches were used. Starting from pollution sources (several like in Milano (S1), or mega site like in Jaworzno (S2)), it is possible to identify the areas more suitable for continuous monitoring of groundwater quality to secure the well protected areas. Therefore, the results of this research can be used by the local government (i.e., environmental agency or public authorities) and decision makers to address groundwater quality issues. Moreover, the obtained maps are useful to avoid human health risk caused by drinking water polluted groundwater and determine if more precautions or measure should be taken to remediate the contamination present in aquifers.