A Multidisciplinary Approach to Evaluate the Effectiveness of Natural Attenuation at a Contaminated Site

: This study evaluates the natural attenuation of chlorinated hydrocarbons as remediation action in a contaminated site downtown the city of Parma (Italy). To achieve this goal, a combination of new investigation methods (bio-molecular analysis, compound speciﬁc isotope analysis, phytoscreening) has been proposed. The approach (named circular multi step) allows to: fully understand the phenomena that occur at the study site, design new investigation activities, and manage best practices. Consequently, each step of the approach improves the conceptual and numerical models with new knowledge. The activities carried out at the study site allowed to detect a contamination of perchloroethylene in a large part of the city of Parma and, of main importance, underneath a kindergarten. The results of the study did not show signiﬁcant natural attenuation of chlorinated hydrocarbons and that the detected contamination could refer to the same unknown contaminant source. Furthermore, the innovative phytoscreening technique was applied to assess the presence of chlorinated hydrocarbons at the ground level. The plume spread was estimated through numerical modeling starting from potential contaminant sources. This study enhances the knowledge of groundwater ﬂow and contamination in Parma and allows authorities to design new investigation/reclamation activities through management actions. 2017, no markers were found in P6, vcrA was detected only in P5 and P19 with middle quantity. Low-middle quantity of DHC-RT were detected at P5, P14, P15, and P19. Dre markers were identiﬁed in P5, P14, P15, P16, and P19. These results allow to state that in December 2017 only slightly active organohalide respiration was detected but no VC reductases was found.


Introduction
At the end of the last century, industries located in urban areas gradually concluded their activities and production or were relocated outside the city boundaries. These activities, as a result of improper processing of dangerous substances or originated by production processes, may have generated a soil and groundwater pollution [1]. The mitigation of groundwater pollution in urban areas is a complex, long-term, and very expensive process. Even with considerable investments, the reclamation may not provide a relevant improvement in groundwater quality. A further difficulty in reclamation processes is the movement of contaminant across the aquifer. Sources of contamination located in downtown cities could affect the quality of groundwater of the downgradient suburbs (contamination goes beyond sites and administrative borders).
The approach of the Groundwater Framework Directive of the European Union (Directive 2006/118/EC) is to identify, to assess, and manage pollution in groundwater bodies. For this reason, EU, within the INTERREG programme [2], funded several projects during the last 15 years. This work is part of the AMIIGA CE32 (Integrated Approach to Management of Groundwater Quality in Functional Urban Area) project; for more details, see [1].
During the last decades, the industrial use of chlorinated hydrocarbons (CHC) has been various and multiform: industrial solvent (dry cleaning and metal stain remover),

Integrated Multi Step Circular Approach
The multi-step circular approach shown in Figure 1 was applied in the present work to evaluate the feasibility of the MNA for CHCs as reclamation action. Figure 1 shows four main steps with different actions (a., b., c., . . . ) to be performed and green lines with arrows indicating the path to follow. The first step is the site characterization that consists in collecting the already available data regarding geology, hydrology, contamination, and potential sources. In order to achieve this goal, existing data (not in digital format) owned by Local and Regional Authorities has been examined and organized into a geographical information system (GIS) database. The data consisted in reports of: reclamation actions closed or ongoing, borehole logs collected for geological purposes, results of sampling campaigns carried out to monitor the quality of groundwater. The GIS allows an easy comprehension of the information and a fast visualization on thematic maps. After the analysis of existing data, new investigation actions, such as boreholes, monitoring wells (MWs), and hydraulic tests are designed. The results of the first step are: the implementation of a new monitoring network and the development of a preliminary conceptual model of the site.
Then, moving to the second step (following the green arrow on left of Figure 1) thanks to the new monitoring network, new data are collected to fully understand the degree of contamination and to evaluate the bacteria that are acting at the study site. To reach this goal, several innovative analyses could be carried out such as biological molecular tools (BMTs), compound specific isotope analysis (CSIA), and phytoscreening. The new data improve the site characterization and consequently the conceptual model. For this reason, it may be necessary to design additional investigation actions and go back to step 1 (follow the green arrow on right of Figure 1).
In order to evaluate the biodegradation, biological molecular tools (BMTs) can be applied [15][16][17][18]. BMTs are used to identify and characterize the bacteria that play an important role in the biodegradation of the contaminants in groundwater. In particular, they evaluate the capability of indigenous microbial consortia to bio-degrade [19][20][21][22][23] CHCs at the study sites [13]. The researchers of the Technical University of Liberec (that was an AMIIGA project partner) developed the guidelines of the BMTs [24] and performed the analyses for all partners. The community bacteria that exist at the monitoring wells have been investigated; in particular, the study focused on organohalide respiring bacteria. In addition, nitrifying and denitrifying bacteria, sulphate-reducing bacteria, iron oxidizing and reducing bacteria, and BTEX (benzene, toluene, ethylbenzene and xylene) degraders were investigated in order to have a complete view of the investigated site. Then, moving to the second step (following the green arrow on left of Figure 1) thanks to the new monitoring network, new data are collected to fully understand the degree of contamination and to evaluate the bacteria that are acting at the study site. To reach this goal, several innovative analyses could be carried out such as biological molecular tools (BMTs), compound specific isotope analysis (CSIA), and phytoscreening. The new data improve the site characterization and consequently the conceptual model. For this reason, it may be necessary to design additional investigation actions and go back to step 1 (follow the green arrow on right of Figure 1).
In order to evaluate the biodegradation, biological molecular tools (BMTs) can be applied [15][16][17][18]. BMTs are used to identify and characterize the bacteria that play an important role in the biodegradation of the contaminants in groundwater. In particular, they evaluate the capability of indigenous microbial consortia to bio-degrade [19][20][21][22][23] CHCs at the study sites [13]. The researchers of the Technical University of Liberec (that was an AMIIGA project partner) developed the guidelines of the BMTs [24] and performed the analyses for all partners. The community bacteria that exist at the monitoring wells have been investigated; in particular, the study focused on organohalide respiring bacteria. In addition, nitrifying and denitrifying bacteria, sulphate-reducing bacteria, iron oxidizing and reducing bacteria, and BTEX (benzene, toluene, ethylbenzene and xylene) degraders were investigated in order to have a complete view of the investigated site.
CSIA has two-fold: (1) it can provide evidence of contaminant's degradation in groundwater, i.e., if biodegradation is occurring, by measuring change in the ratio of stable isotopes in the contaminant [25][26][27][28][29][30][31]; (2) it quantifies the isotopic composition of a specified contaminant and can help in distinguishing sources of organic compounds. In fact, the isotopic signatures of chemicals show differences between manufacturers [27][28][29]32,33] [34][35][36] showed successful applications of CSIA. The guidelines of CSIA [24] and the analyses for all partners were carried out by the researchers of Politecnico di Milano (that was an AMIIGA project partner).
Phytoscreening [37][38][39] technique is based on the ability of the root system of numerous tree species to absorb the chlorinated volatile organic compounds (VOC), and especially of perchloroethylene (PCE), trichlorethylene (TCE), and aromatic hydrocarbons (BTEX) present in the subsoil. Trees can be considered a living pump; the VOCs that are CSIA has two-fold: (1) it can provide evidence of contaminant's degradation in groundwater, i.e., if biodegradation is occurring, by measuring change in the ratio of stable isotopes in the contaminant [25][26][27][28][29][30][31]; (2) it quantifies the isotopic composition of a specified contaminant and can help in distinguishing sources of organic compounds. In fact, the isotopic signatures of chemicals show differences between manufacturers [27][28][29]32,33].  [34][35][36] showed successful applications of CSIA. The guidelines of CSIA [24] and the analyses for all partners were carried out by the researchers of Politecnico di Milano (that was an AMIIGA project partner).
Phytoscreening [37][38][39] technique is based on the ability of the root system of numerous tree species to absorb the chlorinated volatile organic compounds (VOC), and especially of perchloroethylene (PCE), trichlorethylene (TCE), and aromatic hydrocarbons (BTEX) present in the subsoil. Trees can be considered a living pump; the VOCs that are present in the aquifer or in interstitial gases can be absorbed by the root system of trees, transferred to the sap and transported to the trunk, branches, and leaves. The information they provide may allow the carrying out of assessments both for a possible risk assessment and for the identification of the actual values at the points of exposure of the different receptors. The selection of the trees to be sampled must be carried out among species that do not produce plant resins containing volatile terpenes; monoterpenes, such as pinene from conifers, produce problems with instrumentation. The trees that have shown a greater suitability for phytoscreening purposes are those belonging to broad-leaved trees. Populus (poplar) and Salix (willow) but also Quercus pubescens (downy oak), Platanus acerifolia (plane tree), Tilia platyphyllos (linden), Juglans regia (walnut) are suitable specimens. The analytical methods applied to evaluate the collected samples are the same used for the analysis of soils for the identification of VOCs. The first analysis phase involves the extraction of the sample in methanol; while in the second phase, a characteristic volume of vapor is collected from the sampler and analyzed in the gas chromatograph. The sampling of the organic matrix of the tree, was carried out using an incremental auger (Pressler's gimlet, φ = 0.5 cm, length 10-12 cm); the trees were drilled at 1 m from ground level, selecting the sector of the trunk according to the preferential flow direction of the contaminated groundwater. The samples were analyzed in the laboratory by a purge and trap process [40] followed by a gas chromatography-mass spectrometry [41] with detection limit of 0.0005 ppm. After collecting the tree micro-cores, it was possible to verify in situ the presence of contamination in the trees through a semi-quantitative method, colorimetric vials were used to collect the gases at the holes (in vivo sampling) due to the removal of the micro-cores. This sampling technique was developed by ARTA Abruzzo [42] to evaluate the concentrations of contaminants in their gaseous phase in situ.
Once the first two steps reach the convergence, i.e., the expert on the basis of the available budget evaluates that data are enough to develop a reliable conceptual model, a numerical model of the studied aquifer can be developed (third step of Figure 1). The groundwater flow model allows to estimate the mean groundwater flow direction and rate. Moreover, starting the particle tracking from the contaminated monitoring wells backward in time (backtracking), it is possible to delimit the most potential source area [43][44][45]. The results of the numerical model can imply a revision and/or an improvement of the conceptual one and consequently the need of other investigation actions (go back to step 1 or 2, green arrow on left).
In the fourth step (Figure 1), the existing and new data are considered to evaluate the natural attenuation as potential remediation action; moreover, the collected data allow to identify the need of new investigation activities, to design reclamation activities and or management actions that may imply to go back to previous steps.

Study Area
This work takes origin in early 2000 after a reclamation procedure of an old gas station located downtown the city of Parma. Parma is located in northern Italy, in the Emilia-Romagna Region between the Apennine Mountains and the Po Valley ( Figure 2). The area is surrounded by Enza River to the east, Taro River to the west, Po River to the north, and Apennine Mountains to the south. Parma covers a flat area of about 260 km 2 with an average altitude of approximately 57 m above sea level. Parma and Baganza Rivers flow across the study area ( Figure 2). The study area ( Figure 3) is located in the territory of the Municipality of Parma: it is densely urbanized, with a concentration of residential buildings, trading, and services. A contamination of CHCs has been found in the shallow aquifer close to MWs P5 and P24 (see Figure 3).

Hydrogeology
The Quaternary sedimentary deposits ( Figure 4) present ribbon-like and lens-like structures that lengthen following the stream that laid them down and that are mostly made of heterogeneous clastic sediments, from gravelly to clayey (with sudden spatial variations). Locally, the hydro-stratigraphic system is composed of aquifer reservoirs, juxtapose and overlapping, and partially or totally isolated by barriers of permeability made of sedimentary, clayey bodies see Figure 4 for a sketch view. The studied aquifer is the shallow confined one in which groundwater flows from south to north (see Figure 2). The shallow groundwater is recharged upgradient of the urbanized area ( Figure 4) where: i. the shallow aquifer crops out; ii. the rivers are hydraulically connected to the groundwater; iii. the nearby relieves (Apennines) laterally feed the alluvial aquifer [46][47][48][49]. Iacumin et al.  The first significant groundwater level, in the study area, is generally located around 10 m depth (direction N-NE). The investigated aquifer is used for agriculture upgradient the city of Parma. It is important to notice that the aquifer that is used as the water supply for drinking purposes is below 36 m deep the ground level and separated by the shallow one through several meters of clay.
During the AMIIGA project, 7 new monitoring wells (P14-P20) were drilled, up to 25 m from the ground level and screened at different depths in the same shallow aquifer (upgradient and downgradient-SW-NE direction-the study area). During the realization of the monitoring wells, borehole logs have been collected in order to improve the local stratigraphy. The local ( Figure 5) stratigraphic representation, based on existing and new boreholes, shows the following sequence: 2-3 m of backfill, 10-15 m of clay, silty clay, silty sand, and 10-20 m of sandy gravel. Some well tests have been performed at the new monitoring wells to evaluate the hydraulic properties of the aquifer (about 10 −5 m/s) and of the lower confinement (range: 10 −7 -10 −8 m/s). In addition, pumping tests performed on P13 [52] and in an upgradient area [53] confirmed the same order of magnitude of the hydraulic conductivity. The estimated values are summarized in Figure 5. The new MWs have been added to the monitoring network of the shallow aquifer.

Hydrogeology
The Quaternary sedimentary deposits ( Figure 4) present ribbon-like and lens-like structures that lengthen following the stream that laid them down and that are mostly made of heterogeneous clastic sediments, from gravelly to clayey (with sudden spatial shallow groundwater is recharged upgradient of the urbanized area ( Figure 4) where: i. the shallow aquifer crops out; ii. the rivers are hydraulically connected to the groundwater; iii. the nearby relieves (Apennines) laterally feed the alluvial aquifer [46][47][48][49]. Iacumin et al. (2009) and Zanini et al. (2019) [50,51] have confirmed this hypothesis through analyses on stable isotope δ 18 O and δ 2 H collecting samples at rivers, MWs, and precipitation among the study area. The observed δ 18 O and δ 2 H values indicate that rivers and groundwater from the Apennines feed the shallow aquifers within the Parma plain. The first significant groundwater level, in the study area, is generally located around 10 m depth (direction N-NE). The investigated aquifer is used for agriculture upgradient the city of Parma. It is important to notice that the aquifer that is used as the water supply for drinking purposes is below 36 m deep the ground level and separated by the shallow one through several meters of clay.
During the AMIIGA project, 7 new monitoring wells (P14-P20) were drilled, up to 25 m from the ground level and screened at different depths in the same shallow aquifer (upgradient and downgradient-SW-NE direction-the study area). During the realization of the monitoring wells, borehole logs have been collected in order to improve the local stratigraphy. The local ( Figure 5) stratigraphic representation, based on existing and new boreholes, shows the following sequence: 2-3 m of backfill, 10-15 m of clay, silty clay, silty sand, and 10-20 m of sandy gravel. Some well tests have been performed at the new monitoring wells to evaluate the hydraulic properties of the aquifer (about 10 −5 m/s) and of the lower confinement (range: 10 −7 -10 −8 m/s). In addition, pumping tests performed on P13 [52] and in an upgradient area [53] confirmed the same order of magnitude of the hydraulic conductivity. The estimated values are summarized in Figure 5. The new MWs have been added to the monitoring network of the shallow aquifer. In order to evaluate the seasonal variability of the groundwater level, 14 monitoring campaigns were carried out from July 2017 to October 2019 (crosses in Figure 6). In addition, automatic level probes with datalogger were installed at P13, P19, and P24 (see Figure 6). It is clear that the three monitoring points have the same trend, they increase during the spring, reach a peak in early summer, and then they decrease to a minimum value in autumn. It is interesting to remark that the three time series are synchronous; they reach the maximum and minimum peaks at the same time. This analysis confirms the confined nature of the aquifer. In fact, considering that P13 is about 2 km upgradient P19 and about 6 km P24, an unconfined aquifer would not respond so quickly to natural fluctuations. It is important to notice that the minimum reached at the monitoring well P24 in 2017, due to a couple of years of drought, is about 4 m below the minimum reached in 2018. In order to evaluate the seasonal variability of the groundwater level, 14 monitoring campaigns were carried out from July 2017 to October 2019 (crosses in Figure 6). In addition, automatic level probes with datalogger were installed at P13, P19, and P24 (see Figure 6). It is clear that the three monitoring points have the same trend, they increase during the spring, reach a peak in early summer, and then they decrease to a minimum value in autumn. It is interesting to remark that the three time series are synchronous; they reach the maximum and minimum peaks at the same time. This analysis confirms the confined nature of the aquifer. In fact, considering that P13 is about 2 km upgradient P19 and about 6 km P24, an unconfined aquifer would not respond so quickly to natural fluctuations. It is important to notice that the minimum reached at the monitoring well P24 in 2017, due to a couple of years of drought, is about 4 m below the minimum reached in 2018.
In Figure 6, the hydraulic head level at the river station of Baganza River is also depicted, see Figure 3 for the location. The river variability, in the investigated period, is about 1 m. Comparing the hydraulic head monitored at the river and the ones monitored at groundwater monitoring points, it is clear that the river has no effects on the groundwater system at the study site scale. Moreover, P19 is just 200 m far from the Baganza River and no appreciable influences, during the monitoring period, was observed on the groundwater level. the spring, reach a peak in early summer, and then they decrease to a minimum value in autumn. It is interesting to remark that the three time series are synchronous; they reach the maximum and minimum peaks at the same time. This analysis confirms the confined nature of the aquifer. In fact, considering that P13 is about 2 km upgradient P19 and about 6 km P24, an unconfined aquifer would not respond so quickly to natural fluctuations. It is important to notice that the minimum reached at the monitoring well P24 in 2017, due to a couple of years of drought, is about 4 m below the minimum reached in 2018. Figure 6. Groundwater level observed at three monitoring wells through automatic probes and hydraulic level observed at the river station of Baganza River (see Figure 3). The crosses indicate the manual measurements. Figure 6. Groundwater level observed at three monitoring wells through automatic probes and hydraulic level observed at the river station of Baganza River (see Figure 3). The crosses indicate the manual measurements.

Contamination
The first information about groundwater pollution date back to 2002 and was detected close to the monitoring wells P5 and P24. The chemical analyses showed PCE concentration values in groundwater higher than the Italian law limits (1.1 µg/L), even in upgradient monitoring wells such as P20. During the investigated period, the PCE concentration values at monitoring wells increased up to 23 µg/L.
During 2013, a detailed historical analysis (that covers about 25 years) of the commercial activities, which potentially used PCE close to the old gas station has been carried out. Moreover, existing concentrations data at monitoring points were collected in a database.
Nitrite was always close or below the detection and nitrate was always below 50 mg/L (Italian law limits) in all monitoring wells in all sampling campaigns, for further details, see [51]. Ethylene, ethanediol, and trichloroacetate were always below the detection limits. Trichloromethane presents very low values on almost all the MWs. It was detected above the law limits just at P1 and P14 during the winter sampling campaigns but with values below 0.7 µg/L. TCE was always below law limits. 1,1-DCE was not detected, or it was below the law limits in all monitoring campaigns except the one of December 2017. In this campaign 1,1-DCE was above the law limits in almost all the MWs. 1,2-DCE was not detected, or it was below the law limits in all monitoring campaigns. Vinyl chloride was detected only at P16 during the second sampling campaign and at P5 during the fourth sapling campaign, but with very low values (<0.024 µg/L). PCE was the main pollutant observed at the monitoring network; for this reason, in Figures 7 and 8 are respectively shown the PCE observed concentrations in map view and versus time.
It is important to notice that no contamination has been found in MWs P19 and P13, located upgradient P18, during the sampling campaigns. P1 and P16 showed values always below the law limit. P3 was the main polluted MW and unfortunately, it is located in the yard of a kindergarten. Consequently, the regional environmental authority (ARPAE) carried out several additional investigations. In particular, an innovative phytoscreening analysis on trees close to the most polluted area, aimed at indirectly confirming the nontransport by volatilization from the underneath contaminated aquifer to the surface, has been applied.   It is important to notice that no contamination has been found in MWs P19 and P13, located upgradient P18, during the sampling campaigns. P1 and P16 showed values always below the law limit. P3 was the main polluted MW and unfortunately, it is located in the yard of a kindergarten. Consequently, the regional environmental authority (ARPAE) carried out several additional investigations. In particular, an innovative phytoscreening analysis on trees close to the most polluted area, aimed at indirectly confirming the non-transport by volatilization from the underneath contaminated aquifer to the surface, has been applied. Figure 8 shows the concentrations observed at monitoring wells versus time. It is clear that concentrations show a seasonal variation with high values in late winter and low values in late summer. The maximum observed concentration (23 μg/L) was detected at P3 during the last sampling campaign.

Evaluation of Biodegradation
With the aim at characterizing the biodegradation, bio-molecular analyses were carried out during two sampling campaigns (December 2017 and May 2018) on 6 MWs: P5, P6, P14, P15, P16, and P19. The results of the analyses showed that denitrifying bacteria, sulfate-reducing bacteria, and BTEX degraders were reliably detected [51]. The denitrification and BTEX degrading processes are ongoing.  40, which mean: very high abundance, high abundance, middle quantity, low quantity, and not detected bio molecular markers, for additional information see [54,55]. Regarding the sampling campaign of December 2017, no markers were found in P6, vcrA was detected only in P5 and P19 with middle quantity. Low-middle quantity of

Evaluation of Biodegradation
With the aim at characterizing the biodegradation, bio-molecular analyses were carried out during two sampling campaigns (December 2017 and May 2018) on 6 MWs: P5, P6, P14, P15, P16, and P19. The results of the analyses showed that denitrifying bacteria, sulfatereducing bacteria, and BTEX degraders were reliably detected [51]. The denitrification and BTEX degrading processes are ongoing.  15-25, 25-36, 36-40, 40, which mean: very high abundance, high abundance, middle quantity, low quantity, and not detected bio molecular markers, for additional information see [54,55]. Regarding the sampling campaign of December 2017, no markers were found in P6, vcrA was detected only in P5 and P19 with middle quantity. Low-middle quantity of DHC-RT were detected at P5, P14, P15, and P19. Dre markers were identified in P5, P14, P15, P16, and P19. These results allow to state that in December 2017 only slightly active organohalide respiration was detected but no VC reductases was found. DHC-RT were detected at P5, P14, P15, and P19. Dre markers were identified in P5, P14, P15, P16, and P19. These results allow to state that in December 2017 only slightly active organohalide respiration was detected but no VC reductases was found. Regarding the sampling campaign of May 2018, all organohalide-respiring bacteria were decreased or disappeared at all. No CHCs degradation markers have been found except Dre in low quantity at P16. These results indicate an absence of ongoing dehalogenation process.
In order to improve the knowledge on the contamination, CSIA analyses on PCE δ 13 C have been carried out on nine monitoring wells: P3, P5, P6, P14, P15, P16, P17, P18, and P19 during the sampling campaign of December 2017 and May 2018. In December 2017, P16, P17, and P19 and in May 2018, P14, P16, P17, P18, and P19 had a too low PCE concentration value to evaluate δ 13 C. Figure 10 shows the δ 13 C values vs. the observed concentration of PCE for both sampling campaigns. It is possible to notice that, except P18, all the MWs present δ 13 C values between −27 and −25. Looking at Figure 10, it seems that the monitoring wells P3, P5, P6, P14, and P15 have similar δ 13 C values, while P18 shows a different δ 13 C signature. Unfortunately, the low PCE concentration values did not provide any significant result. Further investigations are required.  Regarding the sampling campaign of May 2018, all organohalide-respiring bacteria were decreased or disappeared at all. No CHCs degradation markers have been found except Dre in low quantity at P16. These results indicate an absence of ongoing dehalogenation process.
In order to improve the knowledge on the contamination, CSIA analyses on PCE δ 13 C have been carried out on nine monitoring wells: P3, P5, P6, P14, P15, P16, P17, P18, and P19 during the sampling campaign of December 2017 and May 2018. In December 2017, P16, P17, and P19 and in May 2018, P14, P16, P17, P18, and P19 had a too low PCE concentration value to evaluate δ 13 C. Figure 10 shows the δ 13 C values vs. the observed concentration of PCE for both sampling campaigns. It is possible to notice that, except P18, all the MWs present δ 13 C values between −27 and −25. Looking at Figure 10, it seems that the monitoring wells P3, P5, P6, P14, and P15 have similar δ 13 C values, while P18 shows a different δ 13 C signature. Unfortunately, the low PCE concentration values did not provide any significant result. Further investigations are required. DHC-RT were detected at P5, P14, P15, and P19. Dre markers were identified in P5, P14, P15, P16, and P19. These results allow to state that in December 2017 only slightly active organohalide respiration was detected but no VC reductases was found. Regarding the sampling campaign of May 2018, all organohalide-respiring bacteria were decreased or disappeared at all. No CHCs degradation markers have been found except Dre in low quantity at P16. These results indicate an absence of ongoing dehalogenation process.
In order to improve the knowledge on the contamination, CSIA analyses on PCE δ 13 C have been carried out on nine monitoring wells: P3, P5, P6, P14, P15, P16, P17, P18, and P19 during the sampling campaign of December 2017 and May 2018. In December 2017, P16, P17, and P19 and in May 2018, P14, P16, P17, P18, and P19 had a too low PCE concentration value to evaluate δ 13 C. Figure 10 shows the δ 13 C values vs. the observed concentration of PCE for both sampling campaigns. It is possible to notice that, except P18, all the MWs present δ 13 C values between −27 and −25. Looking at Figure 10, it seems that the monitoring wells P3, P5, P6, P14, and P15 have similar δ 13 C values, while P18 shows a different δ 13 C signature. Unfortunately, the low PCE concentration values did not provide any significant result. Further investigations are required.

Groundwater Flow Model
Starting from the conceptual hydrogeological model described in Section 2.3, at first, a regional groundwater model [51] was developed, see Figure 1, with the aim at evaluating flow directions, fluxes, and boundary conditions for a local groundwater flow and transport model. Then, a local model, focused on the study area, was built, see Figure 1. The local model has two layers (see Figure 5): the first consists in the top confinement (clay-silt layer); the second is the investigated aquifer (sand-gravel). The bottom of the second layer consists in the continuous clay layer confinement that has been considered as no flow boundary. The geometry was built according to new borehole logs ( Figure 4) and existing cross sections [46]. The modeled area is 6.1 km in the south-north direction and 1.5 km in the east-west direction and covers a surface of about 9 km 2 . The local flow model was discretized through a finite difference grid that consist of 610 rows, 150 columns, and 2 layers. Each cell, in the plan view, is square with dimension 10 m × 10 m. The flow simulations were performed by means of MODFLOW 2005 [56]. The regional model [51] allowed to define the boundary conditions of the local model as: assigned hydraulic head at south (upgradient) and north (downgradient) and no flow boundary, at east and west.
Considering the available data, an unsteady state simulation has been performed from 16 December 2017 to 18 November 2018, with stress periods with a length of 7 days each. Upgradient and downgradient boundary conditions are assigned considering hydraulic head measurements carried out at monitoring wells P13 and P1.
Considering that, the investigated aquifer is covered by few meters of clay, the meteoric recharge was not considered. Significant withdrawals, due to irrigation, are observed upgradient the city of Parma; the adopted boundary conditions kept into account the influence of the irrigation on the groundwater levels. Other withdrawals (such as private wells) were assumed negligible with respect to the mean groundwater flow rate.
Boundary conditions (BCs) are of main importance in groundwater modeling. The selection of boundary conditions is the most critical and difficult aspect of the characterization of a groundwater system for conceptual analysis or numerical simulation [57]. Improper boundary conditions cause a non acceptable match between the model response to the hydraulic stresses and the behavior of the natural system. The correct choice allows to properly estimating, in each computation node, the hydraulic heads, velocities, and consequently the fluxes. Moreover, the knowledge of the velocities is crucial in transport modeling. In fact, both advection and dispersion processes depend on flow velocity.
Once the model was set up, the calibration of the hydraulic conductivity process started from the knowledge of few hydraulic tests available on the study area [51,52] and from well test performed on monitoring wells. The parameter estimation process was carried out through the software PEST [58] and evaluated by means of mean absolute error (MAE) and root mean square error (RMSE) [59,60]: where N è is the number of observations,Ĥ i is the computed hydraulic head level at the monitoring point i, and H i is the observed hydraulic head level. Considering 102 observations, MAE and RMSE result respectively 29 and 43 cm. RMSE is less than 3% of the maximum variability in hydraulic heads data; a value below 10% indicate that the model calibration is acceptable [59]. In this case, both the regional [51] and local model are in unsteady state and are able to well reproduce the observed hydraulic heads at monitoring wells. This suggests the proper choice of the BCs.
The estimated hydraulic conductivities (1.8-7.1 × 10 −5 m/s) are consistent with the alluvial nature of the aquifer and are in the same order of magnitude of the pumping and well tests performed [52].

Groundwater Transport Model
A particle tracking backward in time, using MODPATH [61], starting from monitoring wells, has been performed. The backtracking allows to identify the particle path and to delimit potential source areas. Then, the potential sources have been identified superimposing the list of small business companies that are/or were located in the study area and managed CHCs (such as: car garages, laundries, mechanical constructions) [62] to backtracking results. Finally, a simplified transport model was implemented. The transport model was built with MT3DMS [63] to evaluate the plume dimension of a potential source of PCE. As example, a conservative non-reactive contaminant was assumed with a constant injection (0.1 kg/d) at the potential source 222 (Figure 2) for ten years; then the plume evolution was observed for other 20 years. The dispersivity coefficients (longitudinal α L = 10 m, and transversal α T = 1 m) were assumed from literature values [64]. In order to evaluate the effect of advection and dispersion, it is possible to compute the Peclet Number [65] as: where |v| is the magnitude of the seepage velocity vector, L is a characteristic length (commonly taken as the grid cell width), and D is the dispersion coefficient. In advectiondominated problems, the Peclet number has a large value. Considering the longitudinal direction, in this case Pe = 1, that means that the transport is influenced by both advection and dispersion [6]. Figure 11 shows the plume evolution for 30 years with a time step of 5 years. It is clear that there is an important dilution of the concentration and a slow movement of the plume downgradient in addition to a longitudinal and transversal dispersion.

Groundwater Transport Model
A particle tracking backward in time, using MODPATH [61], starting from monitoring wells, has been performed. The backtracking allows to identify the particle path and to delimit potential source areas. Then, the potential sources have been identified superimposing the list of small business companies that are/or were located in the study area and managed CHCs (such as: car garages, laundries, mechanical constructions) [62] to backtracking results. Finally, a simplified transport model was implemented. The transport model was built with MT3DMS [63] to evaluate the plume dimension of a potential source of PCE. As example, a conservative non-reactive contaminant was assumed with a constant injection (0.1 kg/d) at the potential source 222 (Figure 2) for ten years; then the plume evolution was observed for other 20 years. The dispersivity coefficients (longitudinal αL = 10 m, and transversal αT = 1 m) were assumed from literature values [64]. In order to evaluate the effect of advection and dispersion, it is possible to compute the Peclet Number [65] as: where | | is the magnitude of the seepage velocity vector, L is a characteristic length (commonly taken as the grid cell width), and D is the dispersion coefficient. In advectiondominated problems, the Peclet number has a large value. Considering the longitudinal direction, in this case Pe = 1, that means that the transport is influenced by both advection and dispersion [6]. Figure 11 shows the plume evolution for 30 years with a time step of 5 years. It is clear that there is an important dilution of the concentration and a slow movement of the plume downgradient in addition to a longitudinal and transversal dispersion.  The results of the transport models allow to formulate the hypothesis that PCE concentrations observed at P3, P4, P5, P6, and P21 are due to an old contaminant plume (considering the plume development time) originating from a single unknown source. Moreover, the results of the CSIA support this hypothesis; in fact, no significant differences were found in the δ 13 C signatures at the monitoring wells during the two sampling campaigns.

Phytoscreening Results
At the kindergarten, six plants of species suitable for carrying out "in vivo" sampling were identified, two upgradient and four downgradient the site [40][41][42]. The tree samples can be compared with the site chemical data, deriving from the traditional activities of characterization of the soil, subsoil, groundwater, and soil-gas matrices. No contamination was found in the six samples, which means that there is no transport by volatilization of the contamination from the aquifer up to the ground level of the kindergarten.
Thanks to all collected data, ARPAE evaluated the health risk to which people are subject. Luckily, no risk for people has been assessed at the kindergarten.

Evaluation of MNA
BMT analyses showed that a light active organohalide respiration was detected in the sampling round of December 2017; while in May 2018, all organohalide-respiring bacteria were decreased or disappeared at all. In particular, no significant CHC degradation markers have been found. CSIA analyses, due to the low observed concentrations, did not provide additional results. In fact, the δ 13 C signature is very similar for all collected samples during the two sampling campaigns. These results indicate an absence of ongoing dehalogenation process and the need of a continuous monitoring of the concentrations.
Considering that the natural attenuation was not observed, the strength of the contamination (below 10 µg/L in all MWs except P3, P5, and P6 that is between 10 to 23 µg/L), the usage of the studied aquifer and the highly urbanized area, the periodical monitoring of the concentrations at MWs is now the best option for the management of the aquifer.

Conclusions
The assessment of groundwater contamination at a medium-wide scale supported the local authorities (Municipalities, Health and Environment Authorities, and Regional Authorities) and stakeholders to evaluate the current situation and to design/manage remediation activities considering that groundwater contamination moves between administrative borders. Thanks to the AMIIGA project (funded in the INTERREG framework), the multi-step circular approach has been proposed and applied to investigate the shallow aquifer in the city of Parma.
The multi-step approach allows to incorporate new data from one step to the next and return to the previous one if the new knowledge suggests a different conceptualization of the studied aquifer.
The investigations carried out during this work allowed to detect a light contamination of PCE in the shallow aquifer in a large part of the city of Parma. This new knowledge resulted in being very important, because in the center of the studied area, there is a kindergarten. All the collected data were considered to improve the existing hydrogeological model and then to convert it into a numerical one. At first, the flow model allowed to evaluate the groundwater flow direction, seasonality and quantify the fluxes of the investigated aquifer. Then, the transport model allowed to evaluate the spread of the plume starting from a hypothetical contaminant source. Bio molecular analysis showed an absence of natural attenuation of CHCs, while compound specific isotope analysis showed that all the collected samples might refer to the same unknown contaminant source.
The identified contamination of PCE compels local authorities to perform periodical monitoring campaigns, with the aim at reaching several goals: detection of sources and the responsible of the contamination (at this time unknown); evaluation of the effect on health of the increase of contamination at local scale (as carried out for kindergarten); assessment of the plume spread in order to plan the urban development of the town.
Results obtained in this work demonstrated the usefulness of the multistep approach in the site investigation. The outcomes are intended to be an initial contribution for a more exhaustive investigation about distribution of contamination patterns and source location in the shallow groundwater in Parma. The knowledge of the extent of the problem is an essential key-factor for authorities and allows to list new investigation/reclamation activities through management actions [66] and to collect the proper budget.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.