Evaluation of MNA in A Chlorinated Solvents-Contaminated Aquifer Using Reactive Transport Modeling Coupled with Isotopic Fractionation Analysis

: Groundwater contamination by chlorinated hydrocarbons is a worldwide problem that poses important challenges in remediation processes. In Italy, the Legislative Decree 152/06 deﬁnes the water quality limits to be obtained during the cleanup process. In situ bioremediation techniques are becoming increasingly important due to their affordability and, under the right conditions, because they can be more effective than conventional methodologies. In the initial feasibility study phase, the numerical modeling supports the reliability of each technique. Two different codes, BIOCHLOR and PHREEQC were discussed and compared assuming different ﬁeld conditions. Isotopic Fractionation-Reactive Transport Models were then developed in one synthetic and one simple ﬁeld case. From the results, the two codes were in agreement and also able to demonstrate the Monitored Natural Attenuation processes occurring at the dismissed site located in Italy. Finally, the PHREEQC model was used to forecast the remediation time frame by MNA, hypothesizing a complete source cleanup: a remediation time frame of about 10–11 years was achieved by means of natural attenuation processes.


Introduction
Chlorinated Hydrocarbons (CHCs) are globally distributed anthropogenic contaminants that affect the groundwaters at industrial sites. Tetrachloroethene (PCE), Trichloroethene (TCE) and Dichloroethene (DCE) were widely used in dry cleaning and degreasing operations, but had a large variety of other uses [1][2][3]. The Italian Legislative Decree 152/06, which adopts the European Directive 2006/118/CE, introduces the water quality objectives that must be fulfilled as a result of remediation processes, which correspond to the Contamination Concentration Threshold (CCT) for the different pollutants.
Nowadays, the in situ remediation methodologies for contaminated groundwaters are preferred because they avoid off-site handling of contaminated materials, resulting in less expensive operations [4]. Among these, bioremediation techniques are assuming more importance due to the relatively low costs and good removal efficiencies. Bioremediation involves the use of microorganisms to degrade the contaminants into non-harmful compounds; it can be divided into two main categories: Monitored Natural Attenuation (MNA), which is a passive methodology, and Enhanced Natural Attenuation (ENA), which involves active methodologies.
MNA is based on the reliance on natural processes to achieve site-specific remediation objectives within a time frame that is reasonable compared to that offered by other more active methods [5]. The key concepts of this method are to take advantage of preexisting processes occurring at the site and to benefit from the absence of direct actions on the contamination. The processes ascribable to natural attenuation are biodegradation, dispersion, dilution, volatilization, sorption, and biochemical stabilization. Among such processes, biodegradation plays an important role since it involves the reduction in contaminant mass, while the others involve dilution or sequestration processes. The knowledge and understanding of such processes, which are strictly related to the site-specific conditions, is essential to evaluate whether MNA can be a suitable methodology for the site remediation [6]; nevertheless, it requires a detailed site characterization aimed at the evaluation of the contaminant's characteristics and the environment's suitability for the natural attenuation processes [5]. The method minimizes the amount of material to be handled, mainly because all that are needed are the set-up of the observation wells and the simple monitoring of the concentration over time and the personnel exposure to the contaminations, resulting in a low-cost and safe method. For instance, considering the chlorinated solvents case, a common degradation pathway in reducing redox conditions is the reductive dechlorination, a bacteria-mediated process by Dehalococcoides strains [7,8]. In oxic conditions, a possible degradation pathway is provided by the cometabolic oxidation mediated by methanotrophic bacteria, even if it is reported to be significantly effective only in the plume fringe or along sections where the plume discharges to surface water bodies [4,9]. Compound Specific Isotope Analysis (CSIA) can improve MNA interpretation results [10][11][12][13][14][15], supporting the identification of the biodegradation processes and distinguishing them from other natural attenuation mechanisms, by means of the isotopic fractionation effect. However, many manuals and handbooks are present in the literature, encouraging the evaluation and the choice of MNA as a remediation strategy [5,16,17].
MNA is not accepted worldwide as a proper remediation technique, although its popularity is increasing, and more and more countries are including it in their national regulations on remediation processes. In the literature, [16,17] presented a comprehensive analysis of regulations in European countries and the USA, concerning MNA and the occurrence of this strategy as a remediation action. Countries such as the USA and the Flanders region in Belgium introduced MNA as a remediation strategy and often apply it, such as Germany does, but with the difference that is not yet clearly regulated. Finally, countries such as Finland do not recognize it as proper remediation action. However, MNA is generally coupled with other active remediation processes and usually is implemented for residual contaminations or low concentration plumes, mainly due to the quick response needed to protect sensible targets (e.g., a pump & treat procedure to protect water supply wells) that MNA is not able to fulfill.
Several case studies concerning the MNA approach applied alone or combined with other treatments are present: [18] applied the U.S. EPA's screening protocol for MNA evaluation in the Superfund site of Brooklawn (Baton Rouge, LA, USA), where a CHC (chlorinated ethane and ethene compounds) contamination existed, highlighting the importance of implementing a detailed conceptual model. In the study, a Reactive Transport Model (RTM) through the U.S. EPA's BIOCHLOR [19] was also developed and several simulations assuming different scenarios were performed, carrying out that MNA could be one of the feasible alternatives for the site remediation. In the literature, [4] presented the Twin Cities Army Ammunition Plant (TCAAP) case study (Arden Hills, MN, USA), where MNA was successfully applied to clean up the residual contamination by trichloroethene (TCE) and 1,1,1-tichloroethane (1,1,1-TCA) after the application of active remediation strategies (pump & treat, P&T). Additionally, [20] presented the unique Italian case study available in the literature, showing evidence of degradation by natural attenuation for a BTEXS plume in the NIT site and highlighting the need for evaluation protocols before the application of this methodology in formal plans. Some studies include a comprehensive analyses of several plumes managed through MNA [21,22], other studies on Carbon Tetrachloride and Chloroform plume [23], CHCs [24][25][26][27], and hydrocarbons [28][29][30] and abiotic MNA [31].
Numerical modeling enables evaluation and prediction of the results of MNA and ENA approaches: some numerical codes can be applied with this goal and the purpose of the code selection is to identify appropriate models to develop an Isotopic Fractionation-Reactive Transport Model (IF-RTMs) able to simulate MNA or ENA processes. In the literature, a variety of codes are able to reproduce reactive transport conditions is present [32][33][34][35]. Some implement analytical solutions, such as BIOCHLOR [19] and BIOSCREEN [36], while others are numerical codes, such as PHREEQC [37], PHT3D ( [38][39][40][41], which is a code coupling PHREEQC and MT3D [42] a solute transport code based on a MODFLOW [43]), RT3D [44], OpenGeoSys [45,46], TOUGHREACT [47], MIN3P [48][49][50] and many others. However, not all the codes can simulate the isotopic fractionation (IF), e.g., BIOCHLOR-ISO [51,52] and PHREEQC. BIOCHLOR-ISO is an extension of the EPA'S BIOCHLOR, an analytical model able to reproduce both the natural attenuation processes and the isotopic fractionation occurring during biodegradation through an extension of the analytical solution proposed by [53]; the codes also include the stable isotope ratios analysis for carbon and chlorine isotopes in chloroethenes [54,55]. BIOCHLOR-ISO exists in a version that considers only carbon isotopes (2015 version, [51]) and a more recent version where the dual-isotope approach (carbon-chlorine, version 2016, [52]) is implemented. PHREEQC is a 1D geochemical code that allows the definition of customized reactions defining any type of compound as a new chemical species. Such versatility makes PHREEQC a simple numerical model able to reproduce many possible situations occurring in the groundwater. In the literature, many examples of RTMs implemented by means of PHREEQC are discussed. In [56], an IF-RTM for a chlorinated ethenes plume is presented, while [57] discussed a model able to reproduce 3D-CSIA data. Beyond CHCs, PHREEQC, for example, was used to simulate behavior of arsenic in groundwater [58][59][60][61][62], the behavior of soluble reactive phosphorus [63], to model carbon and hydrogen isotopes during propanol and toluene volatilization [64], heavy metal dissolution [65] and to simulate electrokinetic transport and biogeochemical reactions in porous media [66].
Nowadays, it is becoming clear that RTMs and CSIA are important tools for the assessment of bioremediation processes that can help Public Authorities in accepting MNA as remediation strategy. Indeed, the combination of these tools can strengthen the confidence of designers and controllers in the existence of active biodegradation processes, allowing a better evaluation of the feasibility and design of remediation actions. For this reason, it is needed to perform field-case studies in which these tools are applied.
In the present study, the analytical and numerical codes that are suitable for modeling both the reactive transport and the isotopic fractionation (IF-RTM) reproducing the natural attenuation processes were initially selected. Then, an initial comparison between the numerical code (PHREEQC) and the analytical one (BIOCHLOR) was carried out, in order to further validate the first code.
The aim of this study is to evaluate the feasibility of MNA as remediation strategy in a real dismissed contaminated site in southern Italy by means of data interpretation and numerical modelling, in order to demonstrate the advantages of the combined application of RTMs and CSIA. First, a CSM was developed and then a numerical model was implemented both for descriptive and forecast purposes. After an initial calibration of the degradation parameters, the numerical model was used to forecast the remediation time frame as a result of the removal of the pollution source.

Materials and Methods
PHREEQC was selected as the IF-RTM numerical code, as it was deemed suitable to represent the 1D conditions of the site and BIOCHLOR-ISO was selected as the analytical code since it is a well-established tool to evaluate the MNA processes. The MNA is strongly based upon a full conceptual model, borehole drillings, monitoring and characterization of the polluted site; however, a preliminary analysis carried out by means of numerical models can be useful to understand the feasibility of the approach.

Synthetic Case
The first step of this study was to compare the results of PHREEQC and BIOCHLOR-ISO for a benchmark case, concerning a groundwater contamination by chlorethenes. In accordance with the example discussed in [51] based on the work of [67], PHREEQC codes were developed (the code is provided as Supplementary file S4) for two different cases, summarized in Table 1. Since PHREEQC can simulate only 1D transport, the benchmark example was modified to compare the results by setting the transversal parameters (transversal and vertical dispersion) to 0 and the source dimensions equal to unity. The simulation ran for 30 years to reach steady-state conditions. A summary of the model parameters used in each case is presented in the Supplementary Materials (Table S1). The two versions of BIOCHLOR-ISO give the same results for the IF-RTM simulations, while BIOCHLOR-ISO version 2015 is not able to model the solute sorption; therefore, for Case 2 only BIOCHLOR-ISO version 2016 was used. The parameters used in the BIOCHLOR-ISO model are provided in Supplementary file S1. A comprehensive User's Guide for BIOCHLOR-ISO version 2015 is reported as Supplementary file S2. The models developed in PHREEQC need more inputs. First, the chemical species (which are PCE, TCE, DCE and VC) must be defined through the KEYWORDS (specific word of the programming language that defines an internal algorithm in PHREEQC code) "SOLU-TION_MASTER_SPECIES" and "SOLUTION_SPECIES", assigning a name to each of them. For a detailed description of the code and KEYWORDS, see the Supplementary file S4. Second, the initial concentration in each cell at the boundary side (concentrations of the four compounds, called "infilling solution") is defined by the keyword "SOLUTION". Third, through "KINETICS", the stoichiometry and the degradation parameters needed in the calculations must be defined. Fourth, the degradation reactions must be specified: e.g., assuming a 1:1 stoichiometry for each time step, the model calculates, a certain amount of 1,2-cis-DCE produced by the degradation of its parent compound (TCE) and a certain amount of degraded 1,2-cis-DCE, which becomes its daughter compound (VC). This behavior can be modeled following a first order degradation sequential process, as described in [19].
The concentration of a generic compound i, at a given time t, for a known initial concentration C(t = 0), is given by the negative exponential: where λ i is the degradation rate constant of the related compound. The keyword "RATE" defines a reaction for each compound, through a Basic code. The code reads the moles of each compound present in the cell, then reads the degradation parameter from the "KINETICS" keyword and calculates the increase/decrease in moles (depending on the sign defined in the stoichiometry) in the specific time step integrating the equation in time.
Finally, by means of the "TRANSPORT" keyword, the usual groundwater flow and transport parameters are defined.
In Case 1, a more complex code was written because of the presence of different isotopic species. Indeed, the "heavy" and the "light" isotope of each compound must be defined as a separate chemical species to carry out the calculations. The initial concentration of the two isotopes can be calculated from the initial isotopic composition, assuming that the sum of the two carbon (C) isotopes gives the total amount of C: where 13 C and 12 C are the concentrations (in mol/l) of the heavy and light C isotope, respectively, R comp is the isotopic ratio, R STD is the standard isotopic ratio (Vienna Pee Dee Belemnite, shortened as VPDB, 0.0112372), and δ 13 C is the isotopic signature. Next, because the reaction rates change from light to heavy isotope (slower for the last one), the equations must be adapted. First, in KINETICS, the reaction rate is provided for the light isotope species and the enrichment factor (ε) is provided for the heavy one. Then, the equations are written for the two isotopes (modified from [56]).
where, for a generic compound i, 13 C and 12 C are the concentrations (in mol/l) of the heavy and light C isotope, respectively, λ i is the degradation rate constant, α i is the fractionation factor, and t is the time. Finally, the isotopic composition is calculated through Equation (4) once the abundance of the two isotopes is known. In Case 2, sorption mechanisms were considered. Assuming a linear isotherm, the sorption effect was modeled through the introduction of a retardation factor (R d ), calculated as: where ρ b is the bulk density (g/cm 3 ), θ is the porosity, f oc is the fraction of organic carbon and k oc is the partition coefficient (cm 3 /g). The product f oc × k oc is often indicated as K d (distribution coefficient). The BIOCHLOR-ISO 2016 solution only allows the application of an overall retardation factor (R d,model ), which was calculated as the average value of all components, since each compound has its specific retardation factor depending on the partition coefficient (k oc ). To compare the PHREEQC results with the BIOCHLOR-ISO 2016 ones, the overall retardation factor was applied to both the seepage velocity (v r ) and the degradation rate constants (λ i ). Therefore, the retarded velocity v ret and the new degradation rate constants (λ i,ret ) were computed as:

Field Case
The codes were then applied to a real case study, representing an industrial site in southern Italy. The plant started its activity in the early 1960s, remaining in production until 1979 when the activities ceased. During this period, the plant produced VCM/PVC (Monomer Vinyl Chloride and Poly Vinyl Chloride), methanol, and caustic soda. At the beginning of the 1990, all the industrial facilities were demolished, and a first site characterization of the soils and groundwater was carried out. Then, a groundwater monitoring program was activated, and a large number of concentration data of the pollutants was collected, evidencing the current presence of a secondary source composed by TCE, 1,2-cis-DCE and VC. However, based on the historical production information, the primary source was probably caused by a TCE storage facility. Among the different contaminants present in the area, the chlorinated ethenes were the most diffused and therefore were studied in detail. In particular, the conceptual model implemented in the numerical codes concerns the NE side of the site, where the contaminant's plume is well detected by a series of monitoring wells displaced along the natural groundwater flow direction (Figure 1b), running along a slurry wall.

Hydrogeological Setting
The site is located in a valley and has a river to the east, in a fairly flat area. The geology is characterized almost completely by quaternary and recent alluvial deposits in the valley floor, while argillites have emerged on the sides, creating a low permeability boundary. The characterization surveys identified three different stratigraphic layers, as already shown in [68] (Figure 1a):

•
Sandy silts: consisting of recent Holocenic alluvial sediments, this layer is composed of sandy silts and silty fine sands, with an average thickness of 5 m. The hydraulic conductivity was derived from an in situ pumping test and estimated to be around 10 −7 m/s. This layer is unsaturated. • Sandy gravels: also consisting of recent Holocenic alluvial sediments, this layer is composed of polygenic gravels with sub-rounded clasts in sandy or silty sandy matrix. The hydraulic conductivity of the layer was derived from an in situ pumping test and estimated to be around 10 −4 m/s. This layer hosts the phreatic aquifer that is characterized by a saturated thickness between 3 m and 5.30 m; the total thickness of the layer in the area is about 5-10 m. • Silty clay: consisting of Pliocenic sediments, this layer reaches considerable thickness (greater than 50 m). This is considered the bottom of the aquifer (low permeability boundary), since is consisting of not permeable alluvial materials with a hydraulic conductivity estimated (from literature analyses) around 10 −8 -10 −11 m/s.
The groundwater flow is mainly from NO towards SE, with an average hydraulic gradient of 0.005 and an average hydraulic conductivity of about 5 × 10 −4 m/s.
Locally, the piezometric map shows that the water runs along the containment wall turning southwards in correspondence of the piezometer p10, flowing towards p11 and p12. Globally, in the area, groundwater table distance from the ground surface varies from 5 (upgradient of the wall) to 10.5 m.

Analysis of Historical Concentrations and Isotopic Composition
The database of the historical concentrations available for the site spans from 2001 to 2019, however, for the purposes of the present work only the analytic campaigns performed from 2015 were used (Figure 2), due to differences in the analytical methods and the presence of many missing values in the database of the 4 monitoring wells for the previous years. The concentration analysis showed that the wells along the plume path are p9, p10, p11, p12, PN1, PN2 and PN3, however the last three points were excluded from the analyses, since only one analytical campaign and no isotopic composition data were available. Furthermore, the position of well p10 can influence the detected concentration, because it is not perfectly positioned along the hypothesized plume center line. Since no contamination was detected in the wells upstream of p9 and the concentrations in this well were similar to the one measured in well PN1, it was assumed that p9 may likely represents the source concentrations.
The presence of PCE, TCE, 1,2-cis-DCE, VC, and ethene indicates that reductive dechlorination is active in the area. The concentrations of all components, except for the 1,2-cis-DCE, tend to decrease with distance from the source, whereas they remain almost constant in the different piezometers over time (except for seasonal variations or exceptional years). The almost constant values of 1,2-cis-DCE over space can be explained by a balanced production and consumption of the compound, following the reductive dechlorination degradation chain. While the concentrations of PCE are relatively low (even though in p9 it is still above the CCT equal to 1.1 µg/L), the other compounds show significant abundance. TCE and VC concentrations exceed the CCT defined by the Italian law (152/06, 1.5 µg/L and 0.5 µg/L respectively) in each piezometer, while 1,2-cis-DCE is always below the threshold value (152/06 sets a limit of 60 µg/L only for total 1,2-DCE).
The analyses of isotope values is very important and can support the interpretation of the concentration values. Study [69] showed strongly depleted carbon signatures and chlorine and hydrogen signatures similar to literature values of the 2 field cases; this was the proof that by means of a multi-element CSIA approach, the identification of sources and site-specific processes affecting chloroethene transformation in groundwater can be better interpreted. Concerning the measured carbon isotopic composition (shown in Figure 3), the available data is limited to one analytic survey in 2014 for TCE, 1,2-cis-DCE, and VC only (PCE shows very low concentrations). Moving away from the source, there is significant isotopic enrichment (i.e., less negative isotopic composition), indicating active biodegradation of these compounds. During biodegradations, the lighter isotope is consumed more rapidly because it is easier for bacteria to break bonds (metabolize) containing the light isotope [14]. This fractionation process is not significantly influenced by sorption or dilution [70][71][72]. There is enrichment phenomenon for TCE and 1,2-cis-DCE, whereas VC tends to maintain a constant isotopic composition downgradient of well p10, showing that the degradation stops (no more fractionation occurs). The isotopic composition for all compounds is particularly negative compared to the literature values; this specific characteristic is due to the production process that started from a local methane gas reservoir, which is characterized by very depleted levels [73,74]. Concerning the concentrations and isotopic composition analysis, some essential conclusions can be drawn. The concentration values can be assumed to be representative of a steady state condition because the values are constant over the time. Similar to concentration values, the isotopic composition can be assumed to be in steady-state conditions as initial hypothesis, missing new isotopic analyses. A clear presence of the compounds, which belong to the reductive dechlorination chain of chloroethenes, was shown indicating that this biodegradation reaction is occurring. The CSIA shows evidence of biodegradation for all the chlorinated solvents and the abundance of such compounds entails a potentially harmful situation and marks the site as "contaminated", given that the concentration of some compounds is above the Italian CCTs.

Conceptual Site Model Setup
Based both on the hydrogeological setting and on the concentrations and isotopic composition analyses, a Conceptual Site Model (CSM) was developed to represent a simple but comprehensive overview of the situation, in order to show the potential contaminant's pathways, the interfering elements and the potentially sensible receptors. The sandy gravel layer, where the groundwater flows, was assumed homogeneous since no evidence of sub-layers with different permeabilities was found during the field surveys. The source potential location was identified to be a few meters upgradient of well p9, a hypothesis confirmed by a recent additional characterization campaign that shows no contamination in the upgradient wells and soils (Figure 1, clean wells). Some wells were recently drilled close to the historical ones during this characterization procedure (PN1, PN2 and PN3) showing that the plume spreads also eastwards with generally lower concentrations, defining the plume fringe. However, these concentration data were not taken into account during the analysis, since only one measure over the time was available for such points. A leakage of chlorinated solvents from storage tanks or improper removal in the proximity of the production facilities was assumed. Due to the chemical-physical characteristics of these compounds (DNAPL), they either deeply penetrated in the aquifer causing a contaminant pool or were adsorbed on a low-permeability layer, which continuously and slowly released the pollutants in the aquifer (e.g., following a back-diffusion mechanism, such as in [68]). The CHCs move downgradient of the source in accordance with the natural groundwater flow direction along the containment wall, towards the site boundaries, showing important decreases in the concentration values. Therefore, since no active measures are employed at the site, strong natural attenuation effects are evidently active in the area enclosing the above shown piezometers. In particular, the reductive dechlorination almost completely degrades the higher chlorinated compounds into final non-harmful compounds (ethene).

Analytical and Numerical Model Setup
Implementing the analytical and numerical models, some hypotheses are necessary to simplify the complexity of IF-RTM and to compare the two models: • 1D simulations, due to the simple hydrogeological setting • homogeneous aquifer properties • the degradation follows a 1st order kinetic with a 1:1 molar stoichiometry • degradation rate constants and enrichment factors are set constant • the only degradation process is the reductive dechlorination, with a simplified reaction chain, i.e., only TCE, 1,2-cis-DCE, and VC are produced and consumed • the degradation takes action only in the dissolved phase • the source is located in proximity to well p9 and contains TCE, 1,2-cis-DCE, and VC.
PCE was not considered because of the low concentrations. The homogeneous aquifer properties are summarized in Table 2, while the initial estimated values of the degradation rate constants (following [75]) and enrichment factors (set as the median value calculated from the statistical analysis in Table S4, Supplementary file S1) are presented in the (Table S3, Supplementary file S1), together with the initial concentration and isotopic composition in the source (well p9). In particular, the hydraulic conductivity was derived from 11 pumping tests available for the site, whereas the effective porosity and longitudinal dispersivity was achieved using literature analyses. The initial concentration value in p9 for the different compounds was selected to be equal to the average values measured since 2015 (which is the year from which complete and consistent analytical measurements are available), while the initial isotopic composition was set equal to the only available isotopic campaign, assuming the steady-state condition hypothesis. Furthermore, the effect of the sorption was considered. Applying Equation (7), the retardation factors were evaluated for each compound and an average value R d,model was calculated equal to 3.74 (see Table S4 in Supplementary file S1).
The spatial and temporal discretizations in the PHREEQC model are shown in Table 3, while the complete code is available in the Supplementary file S4. In the forecast model, the number of cells was varied and increased to 150, extending the domain up to 1800 m. Furthermore, the shift number was adapted in order to reach the complete cleanup time (about 11 years, when all the concentrations are below the relative CCT) and to get solutions at specific times.
In this case, the sequential dechlorination considers only TCE, 1,2-cis-DCE, and VC, therefore the PHREEQC code was adapted. The initial concentration of the light and heavy isotopes was calculated for each compound following Equations (2)-(4), starting from the isotopic composition measured in p9 (Figure 3).

Model Calibration
The aim of the calibration process is to minimize an objective function that represents a least-squares estimate of the fit between field observations and computed values [76].
In the present study, in order to achieve a good agreement between the computed and the observed values, the model degradation rate constants and enrichment factors were simultaneously calibrated through the optimization code PEST [77], relying on the observed concentrations and isotopic composition data. The assisted calibration is a technique that allows the calibration to be sped up by widely spanning the parameters variability through the minimization of the error variance; it is often applied in optimization procedures of complex flow models, transport models and RTMs [59,[78][79][80][81][82][83][84]. The upper and lower calibration bounds for the degradation rate constants were set as the minimum and the maximum values for each compound indicated by [85], respectively. Boundary values for enrichment factors were set by applying the same logic: a statistical analysis concerning enrichment factors data coming from the BIOCHLOR-ISO database [51] and the website www.isodetect.de (accessed on 18 November 2020) (see Figure S2 and Table S5 in Supplementary file S1 and the database presented in Supplementary file S3) was carried out.
The parameters are subsequently updated by the code in order to minimize an objective function PHI (Equation (10)), defined as the sum of the squared residuals between computed and observed values: where y i and y obs are the concentration or isotope ratio computed and observed values, respectively. The observations were divided in two groups, one for concentration and one related to the isotope ratios, in order to quantify the error associated to each group. The resulting overall objective function value (PHI) was 248.28, with a contribution of the concentrations PHI_conc = 195.83 and of the isotope ratios of PHI_iso = 52.45. A second calibration was performed weighting the observation relative to well p10, considering that the position of this well is not strictly along the plume center line, and thus the detected values may be biased. The initial parameter values were set equal to the first calibration results. Selecting a weighting factor of 0.25 for these observations, the calibration led to an overall PHI value of 135.58, with a contribution of the concentrations PHI_conc = 85.92 and of the isotope ratios of PHI_iso = 49.66.

Predictive Model
The calibrated model was then used to forecast the remediation time frame according to the hypothesis of complete source removal by active methods. The model was divided by three Stress Periods (SPs, Table 4), representing the pre-source removal situation (204 days); the source removal work period (204 days), where the concentrations were assumed half of the initial ones; and the post-source removal period, where the source concentration was set to be null. Varying the last period length (SP3), the remediation time frame was selected to be as long as it was necessary to reach compound concentrations below the CCT set by the Italian law. The PHREEQC code was adapted both assuming a linear trend of the concentrations among the wells and introducing the corresponding concentration as initial concentration in each cell of the model (a total of 29 initial concentrations, one for each model cell). This allowed the initial conditions of SP1 for the presence of a well-developed and stable plume (June 2019) to be correctly reproduced. Moreover, the model domain was extended up to 2000 m from the source (the aquifer extends more than 2 km) in order to check the time frame for which concentrations were below the CCT in the aquifer.

Results
The PHREEQC code was initially compared to BIOCHLOR-ISO for two synthetic cases to validate the code with an analytical solution. In the first case, concentrations and isotopic composition were calculated along a contaminant plume downstream of the source consisting of chlorinated ethenes (PCE, TCE, 1,2-cis-DCE, and VC). Then, in the second case, the two models were compared by considering sorption effects through Breakthrough (BT) curves 200 m downstream of the source, maintaining the same source mix as in case 1. The two models were then used to implement an interpretative model of a real contaminated site, modeling concentrations, isotopic compositions, and sorption effects. Finally, the validated PHREEQC model was used for a predictive analysis in a Natural Attenuation scenario to forecast the remediation time frame.

Case 1
Concerning the synthetic case 1, Figure 4a compares the concentration values in space along the plume center line while Figure 4b shows the isotopic composition calculated by the two models. The simulated concentrations over the time and space achieved through PHREEQC agree almost perfectly to those achieved by BIOCHLOR-ISO, as well as the isotopic composition calculated by the two models. The computed RMSE between the two results is 0.0062 mM for the concentrations and 0.17‰ for the isotopic composition. These results show good agreement, therefore PHREEQC can be considered a suitable and reliable code for RTM simulations.

Case 2
Concerning the synthetic case n.2, Figure 5 shows the results of the comparison between the BT curves computed by PHREEQC and BIOCHLOR-ISO version 2016 at 200 m downgradient of the source. Figure 5a shows the BT curves without implementing the sorption effects, while Figure 5b considers sorption mechanisms through the implementation of a retardation factor R d = 2.14. Also, in this case, the model results show good agreement with an RMSE equal to 0.01 mM in both cases.

Interpretative Model
The first important result applying RTM to the real case was the identification of the most reliable degradation rate constant values by means of a calibration process. Both the trial-and-error and the assisted calibration (PEST) approaches can be used, but to reduce the time expenditure, the second approach was chosen. The calibration process is a crucial and necessary step to properly set the model for any following predictive simulations. Figure 6 shows the concentrations resulting from the model after the calibration process of the degradation rate constants in the Best Fit (BF) parameters case and the weighted (W) case. The results of BIOCHLOR-ISO and PHREEQC agreed well also in the real case (the comparison in the Best Fit case is shown in Supplementary file S1, Figure S1). The associate calibrated values and corresponding RMSE before and after the calibrations are shown in Table 5.  The results of the calibration of the enrichment factors are shown in Figure 7, while the calibrated values and the corresponding RMSE are shown in Table 6.  The calibration process satisfactorily adapted the parameters to the field-measured concentrations: the comparison between the simulated curves and the measured values agrees quite well and an improvement related to the pre-calibration simulations is clear. However, the VC isotopic values computed by PHREEQC show a first enrichment of the isotopic composition but, after p10, the curve decreases showing more depleted values. This behavior can be explained as an effect of VC production and consumption balance. Around p9, the VC shows higher concentrations that rapidly decrease towards p12, while the 1,2-cis-DCE shows almost constant values (around 40 µg/L) throughout the whole domain. Therefore, a huge amount of VC is degraded in the first part of the domain, resulting in an enrichment of the isotopic composition. However, downgradient of p10, the production of depleted VC by the 1,2-cis-DCE degradation results in an overall decrease of 13 C concentration and, therefore, in more depleted values, effect that is enhanced by the higher enrichment factor of 1,2-cis-DCE.
The PHREEQC results were also compared with the analytical solution of BIOCHLOR-ISO achieving good agreements; therefore, this is an additional validation of the numerical code for a real case.
The calibrated parameters were then compared with literature values to verify whether the values were acceptable in comparison with the existing field measurements. The degradation rate constants were compared with the values reported by [85], which performed a statistical analysis on 185 field or laboratory measured degradation rate constant data. All the calibrated parameters are included in the box-plot within the 25 • and 75 • percentile and are quite close to the median values (except for the TCE value, which is closer to the 75th percentile). Therefore, the calibration process can be considered consistent with the literature values.
In order to compare the calibrated enrichment factors with literature values, a statistical analysis was performed relying on 179 enrichment factors available in the literature (See Table S5 in Supplementary file S1), mainly from microcosm studies, both in aerobic and anaerobic conditions. The comparison is shown in Figure 8.  (Table S5 in Supplementary file S1) performed on the enrichment factor database presented in Supplementary file S3 (elaboration of data from [51] and www.isodetect.de, accessed on 18 November 2020), compared with calibrated enrichment factors in the Best Fit case (green dots).
The calibrated TCE enrichment factor is acceptable and comparable with anaerobic degradation. The 1,2-cis-DCE calibrated enrichment factor falls on the lower whisker of the anaerobic boxplot. This value seems to be extreme but is nonetheless acceptable. Moreover, the modeled isotopic composition of 1,2-cis-DCE shows the biggest error among the different compounds, as previously discussed, indicating complex degradation. Finally, the VC enrichment factor value falls between the boxes relative to the aerobic and anaerobic conditions, indicating that the degradation conditions can be mixed.

Predictive Model
Once the RTM parameters were calibrated, according to the detected concentrations and isotope signatures, PHREEQC was applied to predict the natural attenuation effects on chloroethenes following a hypothetical source removal. In Figure 9, the results for three different time frames show the concentrations change after the source remediation, comparing the achieved solution using the parameters calibrated for the Best Fit case (BF) and the weighted case (W).
The graphs show an important reduction in the concentrations along the plume line, following the contamination source removal: 500 m downstream of the source, after just 5 years TCE and VC concentrations are less than 3 µg/L and just above the CCT. A complete site remediation is predicted to be achieved after 3876 days (10.6 years) in the weighted calibration parameters case (W), whereas after 4080 days (11.2 years) in the Best Fit (BF) parameters case. Conversely, considering a non-degradative model, the concentrations decrease less rapidly with a remediation time greater than 30 years (see Figure S2 in Supplementary file S1). Suffix BF refers to the Best Fit parameters case, while suffix W refers to the weighted parameters case. The limit values are 1.5 µg/L for TCE and 0.5 µg/L for VC. The limit for 1,2-cis-DCE is not shown because, according to the Italian regulation, the limit is given only for total 1,2-DCE (cis + trans, 60 µg/L); however, 1,2-cis-DCE always shows concentrations below the CCT.

Discussion
The reported case study shows how the design of the remediation actions needs to assess the presence of Natural Attenuation processes and include them in the prediction modeling. The aim of the paper is to show that even simple 1D models can help to assess the field conditions and can be used only as a preliminary feasibility tool for Natural Attenuation remediation procedures. Additionally, 1D models can be useful for simple hydrogeological settings, such as the one discussed here. Nevertheless, it is important to follow a series of steps to achieve reliable predictions. After the site hydrogeological characterization, the analysis of monitored concentration over space and time, as well as the identification of the contaminants' attenuation and degradation products, represents the first essential step. Moreover, CSIA can provide an added value in confirming the presence of degradation processes, excluding that the concentration reduction is given by mechanical dispersion, dilution, sorption or volatilization. Monitoring the isotopes in water samples is not a common or widespread technique but is often considered an added value. Indeed, once this information is available, it can be used to support the remediation design starting from the modeling phase.
In the real case presented herein, the site characterization indicates the presence of chlorinated ethenes, mainly TCE, 1,2-cis-DCE, and VC, whose source is assumed to be located upgradient of well p9. The supporting CSIA analyses on C isotopes and the concentration trends clearly indicate the occurrence of a steady-state plume that is undergoing biodegradation due to reductive dechlorination processes. The hydrogeological setting is simple because the aquifer is included into a uniform sandy gravels layer, which rests on a silty clay layer at the bottom. The silty clay layer could have produced the back-diffusion phenomenon over time and from concentration and isotopes analyses, the situation can be approximated as in steady-state condition. Therefore, the hydrogeological and contamination background discussed makes the site suitable for a combined remediation approach between an active source removal method and MNA to clean up the residual concentrations.
Then, a reliable numerical code is necessary for the modelling phase. The present study confirms that PHREEQC is a multidisciplinary code that can be adapted to reproduce the Natural Attenuation processes in aquifers. During the design phase, this code is indeed a useful tool to support both the remediation techniques based on biodegradation and the best possible alternative. The modeling phase requires a proper calibration of the transport parameters. In the present case, the adopted simultaneous calibration (carried out using PEST) of the degradation rate constant and isotope enrichment factor represents an added value because it considers the relative effect that one parameter group has on the other. Using the assisted calibration instead of a trial-and-error approach allows easy determination of an optimized solution and is a method seldomly used in IF-RTMs calibration.
The calibrated degradation rate constant values were compared with the statistical analysis performed by [85], with the selected values within the 25th and 75th percentiles of the population and comparable to anaerobic degradation. Moreover, the calibrated enrichment factors were compared with the statistics calculated on the database presented as Supplementary file S3. The TCE values were coherent with anaerobic degradation, while the ones for VC and 1,2-cis-DCE show values that may involve mixed conditions. In particular, the 1,2-cis-DCE enrichment factor was close to the lower extreme boundary of the aerobic conditions and outside the lower whisker of the anaerobic values. This issue cannot be neglected since two factors could affect the results. First, only one isotopic survey is available for the site (in 2014) and second, the enrichment factors database includes many values derived from both laboratory and field studies, which are strongly dependent on specific site conditions. Therefore, no values can be excluded, even if they fall outside the whisker. To properly design an executive remediation project, a new isotopic survey and new microcosm experiments should be executed in order to further constrain the model parameter ranges.
A potential limitation of this model is the simple 1D simulation, which could neglect the transversal dispersivity, which in some cases determines a reduction in the mass flowing along the plume center line. This influences the values of the targets used during the calibration process and, consequently, PEST adapts the degradation rate constant values, adjusting them to represent the transversal dilution effect. This issue affects the absolute values of the degradation rate constants but not the solution reliability, since the concentration reduction due to transversal dispersion is included. However, the 1D approach allows easy representation of the reactive transport occurring at the site, owing to the homogeneous hydrogeological setting, and contemporarily interprets and represents the simple Natural Attenuation processes occurring at the site quite well. For a more complex hydrogeological case, the model can be further improved through representation of 2D conditions by introducing trace tests to better estimate the dispersivity, by planning more isotopic campaigns and extending the CSIA to chlorine atoms using a dual-isotope approach, or by more thoroughly investigating the biodegradation of 1,2-cis-DCE and VC.
If MNA is selected as the remediation strategy, in general and for more complex cases, it should be accompanied by a structured monitoring network and planning. For example, in order to better characterize the site, new wells can be drilled downgradient of p12 in order to check that no contaminants exit from the site boundary and a backup system can be provided (e.g., Pump-and-treat P&T), in case MNA does not work as expected. However, the model shows that in simple geological contexts MNA can be a preferable technique to active strategies, such as P&T, as it is both an effective and a low-cost option.

Conclusions
This paper evaluates the applicability of MNA as possible a remediation strategy for a site contaminated by chlorinated solvents by means of a 1D numerical model in a field case where hydrogeological conditions are simple. First, the analytical model BIOCHLOR-ISO and the numerical code PHREEQC were compared assuming two different conditions, i.e., considering an IF-RTM and the sorption effects. The comparison showed good agreement in each case, for both the concentrations and the isotopic composition; therefore, PHREEQC is a reliable code that can be used for IF-RTMs and is more appropriate than BIOCHLOR-ISO because it is able to simulate transient conditions. Moreover, this PHREEQC property allows the coupling with common 2D or 3D hydrogeological codes in order to carry out more detailed results.
Subsequently, a CSM of a chlorinated solvents-contaminated site located in southern Italy was developed and converted into a 1D steady-state IF-RTM. The degradation rate constants were initially estimated by means of the method presented by [86] and the enrichment factors were selected carrying out a statistical analysis on 179 values collected from literature research. The flow and transport parameter values were derived from field data. Then, the numerical models were calibrated by using the optimization code PEST, calibrating the degradation rate constants and the enrichment factors simultaneously.
The calibrated model was finally used to forecast the remediation time frame. Assuming the hypothesis of complete source removal, the remediation time frame was determined to be about 11 years, making the site suitable for MNA application. The computed remediation time is reasonably short and a MNA procedure can be easily applied in this case because of the implicit hypothesis and of the simplicity of the site (hydrogeological setting, few monitoring wells, etc.). The solution was compared with a non-degradative case, showing that the remediation time frame without degradation would have been much longer than 30 years. This is further evidence of how RTM can provide an estimate of remediation times with or without degradation processes, in order to help engineers or agencies to clean up polluted sites.
In conclusion, even if MNA is not yet diffused as a common remediation strategy in Italy, this application study shows the potential and value, which can be relevant for similar industrial sites. IF allows researchers to confirm that degradation (bond-braking reactions), not only dispersion, dilution and sorption are active; therefore, this strongly supports MNA. This study provides a stepwise methodology that can be applied to preliminary evaluate the feasibility of the MNA approach application. For this case, the model is also able to provide an estimation of the remediation time frame that can be helpful to site owners and Public Authorities in the process of selecting a suitable remediation strategy that is both effective in reasonable time frames and as low-cost as possible.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10.339 0/w13212945/s1, Supplementary file S1 contains supplementary figures and tables: Figure S1: Comparison between the PHREEQC solution (full line) and the BIOCHLOR-ISO solution (dashed line) in the Best Fit parameters case; Figure S2: Forecast model result for 30 years in the non-degradative case; Table S1: Model parameters for the synthetic case, Table S2: Concentration and isotopic composition of the source, Table S3: Initial values of the degradation rate constants and enrichment factors, source concentrations and isotopic compositions of the field case, Table S4: Statistical analysis of the enrichment factor database. In Supplementary file S2 is presented the BIOCHLOR-ISO User Manual. In Supplementary file S3 is presented the enrichment factors database. Supplementary file S4 contains the PHREEQC codes for the synthetic case, the descriptive model, and the forecast model.

Data Availability Statement:
The data used in the manuscript are partially available, due to privacy issues of contaminated sites. The data about the isotopic enrichment factors are available at the webpage www.isodetect.de, (accessed on 18 November 2020) and in [51], however in Supplementary file S3 are present all the references to the single articles reporting the values. The data about the synthetic cases were taken from [51].