Analytical and Numerical Methods for a Preliminary Assessment of the Remediation Time of Pump and Treat Systems

: Several remediation technologies are currently used to address groundwater pollution. “Pump and treat” (P&T) is probably one of the most widely applied, being a process where contaminated groundwater is extracted from the subsurface by pumping and then treated before it is discharged or reinjected into the aquifer. Despite being a very adaptable technology, groundwater remediation is often achieved in long and unsustainable times because of limitations due to the hydrogeological setting and contaminant properties. Therefore, the cost–beneﬁt analysis over time results in an ine ﬃ cient system and a preliminary evaluation of the clean-up time is crucial. The aim of the paper is to compare, in an integrated manner, the application of some models to estimate the time to compliance of a P&T system in relation to the speciﬁc hydrogeological condition. Analytical solutions are analyzed and applied to an industrial site and to a synthetic case. For both cases, batch ﬂushing and the advection-dispersion-retardation (ADR) model underestimate remediation times comparing the results to real or simulated monitoring data, whereas the Square Root model provided more reliable remediation times. Finally, for the synthetic case, the reliability of analytical approaches and the e ﬀ ects of matrix di ﬀ usion are tested on the basis of a numerical groundwater transport model speciﬁcally implemented, which conﬁrm the results of the analytical methods and the strong inﬂuence of the matrix di ﬀ usion on the results. Author Contributions: Conceptualization, L.A.; methodology, F.R. and M.A.; software, F.R. and M.A.; validation, L.A., F.R., and M.A.; investigation, F.R. and M.A.; resources, L.A. and F.R.; data curation, M.A.; writing—original preparation, M.A.; writing—review and editing, L.A. and F.R.; supervision, L.A.; project administration,


Introduction
In Western Europe only, more than 300,000 sites whose aquifer need be remediated have been identified [1]. In Italy, 12,482 sites under regional jurisdiction, together with 41 sites under national jurisdiction (Site of National Interest, SIN), are currently subjected to administrative proceedings [2].
Groundwater contamination in industrial sites is mainly due to leakages from plants or tanks and to the incorrect disposal of untreated wastes on soils. In general, the main pollutants found in groundwater are aromatic hydrocarbons (e.g., Benzene, Toluene, Ethylbenzene and Xylene, BTEX compounds) and chlorinated hydrocarbons (e.g., Tetrachloroethylene, Trichloroethylene) present also as Non-Aqueous Phase Liquids (NAPLs, [3][4][5]), as well as heavy metals.
A correct and thorough definition of the Conceptual Site Model (CSM) of the polluted area, obtained by organizing all of the available data and information about the contamination and the hydrogeological setting and the potential targets, is crucial to select the most appropriate remediation technique. Currently, many remediation techniques are available, classified according to where the contaminated matrix is treated (in situ, on site or off site) and how the contamination removal is achieved (physical, chemical, or biological) [6]. system design phase and in the performance evaluation of an existing system. When the transport phenomena are driven by advection and sorption, the batch flushing model [38] is noteworthy. In this model, the physical process of water renewal due to a pumping system is simulated by means of a subsequent series of flushing stages with clean water in sufficient quantity to fill the pore volume within a given portion of the aquifer. The batch flushing model was applied to the Lone Pine Superfund site in the US, where a simple aquifer system exists and the contamination was due to dissolved chlorinated solvents [39]. The batch flushing model becomes more effective when coupled with particle tracking analysis, although it involves numerical simulations [40]. Another interesting analytical approach is the advection-dispersion-retardation (ADR) model developed by [41], where clean water flows through the contaminated area and the resulting contaminant flushing within the control volume is simulated using an adjusted version of the Ogata-Banks equation. This model can be applied when the diffusion phenomenon is negligible in the transport processes. Nevertheless, in most contaminated sites, the remediation time depends on matrix diffusion processes. Matrix diffusion was widely studied by [20,[42][43][44], which by means of contaminant numerical transport models showed how this phenomenon is the cause of long clean-up times when low hydraulic conductivity zones are present. In cases where fine sediments are present and matrix diffusion is a relevant mechanism, the square root model (SRM) and the Dandy-Sale model (DSM) [21,[45][46][47], as implemented in the Matrix Diffusion Toolkit (a freeware Microsoft Office-based tool) [48], can be used to assess the performance of the P&T in the long-term.
In this paper, some available methodologies to assess the remediation time of a P&T system were applied to a real field case and to a synthetic case. Analytical and easy to implement solutions were chosen depending on the contamination transport phenomena occurring at the site; the batch flushing model, the ADR model, and the SRM were applied both to the field case and benchmark case. Some ad hoc numerical models were implemented only for the benchmark case (in MODFLOW/MT3DMS (version 2000, United States Geological Survey, Reston, VA, USA)) in order to study the back diffusion role on clean-up time, since low hydraulic conductivity zones are often present in real cases. The numerical simulation results allowed for the verification of the back diffusion effects and to validate the analytical SR model approach.

Field Case
The first case study refers to a former chemical industry, located in Italy, that from the early 1960s until 1990 manufactured PVC (polyvinyl chloride), chlorine, and sodium hydroxide, among others. After site decommissioning, the site characterization activities carried out highlighted a high groundwater contamination due to chlorinated hydrocarbons (mainly 1,1-dichloroethylene or 1,2-dichloroethylene or vinyl chloride). This fact led to the implementation of a physical source control measure in a portion of the contaminated area.
The site stratigraphy can be described as follows: (1) Sandy silts and fine silty sands with an average thickness of about 3-5 m and of brown-light brown color; (2) Sandy gravels in a sandy or silty-sandy matrix, of gray light brown color, with an average thickness of about 5-10 m; (3) Silty clays of gray-light blue color with an average thickness greater than 50 m, with thin sandy layers.
The silty clay layer acts as a base for groundwater, the latter being mainly hosted in the sandy gravel unit; the shallow silty sandy layer is an unsaturated (or vadose) zone, since no groundwater flow is present. The data recorded during several piezometric and analytical field campaigns allowed definition of the hydrogeological setting and groundwater quality. The flow is mainly directed towards North North West-South South East, where a river is located. allowed definition of the hydrogeological setting and groundwater quality. The flow is mainly directed towards North North West-South South East, where a river is located. Figure 1a shows the site stratigraphy while Figure 1b shows the piezometric map obtained by interpolating the hydraulic head values measured in the piezometers; the map shows only the piezometers relevant from a contamination point of view; the other points used for the piezometric interpolation are not represented. Despite the physical source control implementation near the monitoring wells p1, p2, p3, p5, p7, and p8 (see Figure 1), concentrations of 1,1-DCE, 1,2-DCE, and VC exceeding regulatory limits were continuously observed. Comparing the contaminant distributions over time in each piezometer, no decrease in the detected concentrations was observed, but a fluctuation around an average value. Based on the site characterization analyses, the presence of an active source (probably active for many years) is likely assumed in the proximity of p1; therefore, the aquifer contamination can be supposed to be in a steady-state condition. For the sake of simplicity, the analyses discussed in this paper refer to one contaminant only, 1,2-DCE, whose groundwater remediation limit is 60 μg/L, according to the Italian regulation. The analysis was specifically focused on 1,2-DCE because the concentrations were higher than the legal limit, unlike 1,1-DCE and VC, whose concentrations were just above the legal limit. The maximum 1,2-DCE concentration in the last 20 years was 20,000 μg/L, measured in piezometer p2 in 2006. Table 1 shows the average values of the parameters implemented in the analytical model calculations. Kd represents the distribution coefficient, which for hydrophobic, nonpolar organic contaminants is derived from: where Koc is the organic carbon partition coefficient for 1,2-DCE and foc is the organic carbon content of the aquifer. Both parameters can be obtained from laboratory analyses of core material and literature sources. The retardation factor is linked to the lithology and is calculated as: where ρ is the bulk density and θ the porosity. Despite the physical source control implementation near the monitoring wells p1, p2, p3, p5, p7, and p8 (see Figure 1), concentrations of 1,1-DCE, 1,2-DCE, and VC exceeding regulatory limits were continuously observed. Comparing the contaminant distributions over time in each piezometer, no decrease in the detected concentrations was observed, but a fluctuation around an average value. Based on the site characterization analyses, the presence of an active source (probably active for many years) is likely assumed in the proximity of p1; therefore, the aquifer contamination can be supposed to be in a steady-state condition. For the sake of simplicity, the analyses discussed in this paper refer to one contaminant only, 1,2-DCE, whose groundwater remediation limit is 60 µg/L, according to the Italian regulation. The analysis was specifically focused on 1,2-DCE because the concentrations were higher than the legal limit, unlike 1,1-DCE and VC, whose concentrations were just above the legal limit. The maximum 1,2-DCE concentration in the last 20 years was 20,000 µg/L, measured in piezometer p2 in 2006. Table 1 shows the average values of the parameters implemented in the analytical model calculations. K d represents the distribution coefficient, which for hydrophobic, nonpolar organic contaminants is derived from: where K oc is the organic carbon partition coefficient for 1,2-DCE and f oc is the organic carbon content of the aquifer. Both parameters can be obtained from laboratory analyses of core material and literature sources. The retardation factor is linked to the lithology and is calculated as: where ρ b is the bulk density and θ the porosity. Seepage velocities vary according to the different hydraulic condition implemented. In one case, the unperturbed aquifer conditions with a hydraulic gradient of 0.44% were simulated, whereas in the second case, the hydraulic gradient (equal to 0.45%), calculated considering the whole length of the contaminated area, was slightly influenced by a hypothetical hydraulic barrier made up of eight wells (total withdrawal of 20 m 3 /s) located on a 300 m line between piezometers p4 and p6, approximately 320 m downgradient of p2.

Synthethic Case
The second case study represents a synthetic case, specifically implemented to evaluate the clean-up time in the presence of low permeability (clay) layers. The simulated domain reproduces a rectangular aquifer portion 1 km long, 450 m wide, and 11 m thick. The aquifer material is mainly composed of sand, with the presence of a single square clay lense with an area of 100 × 100 m 2 and a thickness of 2.5 m, located 150 m downgradient of the model domain's west boundary ( Figure 2).  Seepage velocities vary according to the different hydraulic condition implemented. In one case, the unperturbed aquifer conditions with a hydraulic gradient of 0.44% were simulated, whereas in the second case, the hydraulic gradient (equal to 0.45%), calculated considering the whole length of the contaminated area, was slightly influenced by a hypothetical hydraulic barrier made up of eight wells (total withdrawal of 20 m 3 /s) located on a 300 m line between piezometers p4 and p6, approximately 320 m downgradient of p2.

Synthethic Case
The second case study represents a synthetic case, specifically implemented to evaluate the clean-up time in the presence of low permeability (clay) layers. The simulated domain reproduces a rectangular aquifer portion 1 km long, 450 m wide, and 11 m thick. The aquifer material is mainly composed of sand, with the presence of a single square clay lense with an area of 100 × 100 m 2 and a thickness of 2.5 m, located 150 m downgradient of the model domain's west boundary ( Figure 2). The hydraulic head boundary conditions were implemented as constant heads conditions, of 11 and 9.6 m a.s.l. for the western and eastern boundaries, respectively. Such conditions impose a groundwater flow with a west-east direction and a hydraulic gradient of 0.14%.
In the aquifer, a Trichloroethylene (TCE) spill of 500 kg dispersed over a square area of 36 m 2 was assumed to occur during a period of 5 years. TCE was chosen because it is a widespread contaminant and therefore could well represent many real situations; in any case, the choice of another compound does not affect the overall methodology. The hydraulic head boundary conditions were implemented as constant heads conditions, of 11 and 9.6 m a.s.l. for the western and eastern boundaries, respectively. Such conditions impose a groundwater flow with a west-east direction and a hydraulic gradient of 0.14%.
In the aquifer, a Trichloroethylene (TCE) spill of 500 kg dispersed over a square area of 36 m 2 was assumed to occur during a period of 5 years. TCE was chosen because it is a widespread contaminant and therefore could well represent many real situations; in any case, the choice of another compound does not affect the overall methodology. The simulated scenarios follow: (1) Scenario A: the contaminant spreads vertically into the hydraulically transmissive material (sand) until it reaches the top of the clay lense; then it proceeds horizontally according to groundwater flow. The whole TCE mass is released in the transmissive zone ( Figure 3a). (2) Scenario B: the contaminant spreads vertically into the sand up to the top of the clay lense, then it accumulates into the clay layer; 96% of the TCE mass remains in the transmissive layer while 4% migrates into the clay lense ( Figure 3b). (3) Scenario C: 96% of the TCE mass is inside the sand aquifer, while 4% accumulates in a small central portion of the clay lense with thickness of 0.5 m (Figure 3c). (4) Scenario D: no clay lense is present and the aquifer is therefore homogeneous; the aquifer thickness is 4.5 m and the whole TCE mass is released into the sand (Figure 3d).
Water 2020, 12, x FOR PEER REVIEW 6 of 25 The simulated scenarios follow: (1) Scenario A: the contaminant spreads vertically into the hydraulically transmissive material (sand) until it reaches the top of the clay lense; then it proceeds horizontally according to groundwater flow. The whole TCE mass is released in the transmissive zone ( Figure 3a). (2) Scenario B: the contaminant spreads vertically into the sand up to the top of the clay lense, then it accumulates into the clay layer; 96% of the TCE mass remains in the transmissive layer while 4% migrates into the clay lense ( Figure 3b). (3) Scenario C: 96% of the TCE mass is inside the sand aquifer, while 4% accumulates in a small central portion of the clay lense with thickness of 0.5 m (Figure 3c). (4) Scenario D: no clay lense is present and the aquifer is therefore homogeneous; the aquifer thickness is 4.5 m and the whole TCE mass is released into the sand (Figure 3d). The contamination source is assumed to be active for five years after which it is instantaneously removed. However, the concentrations in the aquifer downgradient of the source remain high, so two remediation solutions are hypothesized: (1) Solution 1: monitored natural attenuation (i.e., no remediation system is activated), and the contamination evolution is monitored by means of a piezometer monitoring network. (2) Solution 2: a P&T system with one pumping well located 140 m downgradient of the source and screened over the entire aquifer thickness is implemented together with monitoring piezometers. For scenarios A, B, and C, the pumping rate is 40 m 3 /day, while for scenario D is 7.2 m 3 /day.
The hydrogeological and physical-chemical properties of the aquifer material and the contaminant are shown in Table 2. The contamination source is assumed to be active for five years after which it is instantaneously removed. However, the concentrations in the aquifer downgradient of the source remain high, so two remediation solutions are hypothesized: (1) Solution 1: monitored natural attenuation (i.e., no remediation system is activated), and the contamination evolution is monitored by means of a piezometer monitoring network. (2) Solution 2: a P&T system with one pumping well located 140 m downgradient of the source and screened over the entire aquifer thickness is implemented together with monitoring piezometers. For scenarios A, B, and C, the pumping rate is 40 m 3 /day, while for scenario D is 7.2 m 3 /day.
The hydrogeological and physical-chemical properties of the aquifer material and the contaminant are shown in Table 2.

Analytical and Numerical Approaches
Among all the analytical methods and numerical models available in the literature to evaluate the remediation time of a P&T system, the easiest to implement and most efficient analytical procedures were selected to assess the feasibility of a designed P&T. The theory behind the chosen approaches (batch flushing model, ADR model, and the square root model) is hereafter presented. Finally, the synthetic case implementation in a numerical model is discussed.

Batch Flushing Model
The batch flushing model is implicitly based on the number of pore volumes, and accounts only for advection and sorption processes, neglecting dispersion. The pore volume (PV) is the volume of groundwater within a contamination plume, which is defined as: where b is the plume thickness (-), θ is the formation porosity (-), and A is the area of the plume (m 2 ). If the thickness is relatively uniform, then the following relationship applies: where B is the average thickness of the plume (m). The batch flushing model [38,39] is a useful approach for estimating the clean-up time in simple aquifer systems with chemicals for which interaction with the solid matrix can be represented by linear sorption. It simulates a situation in which clean water is circulated through a region that initially contains contaminated ground water; the approach assumes simple advective displacement of contaminated water, thus neglecting dispersive transport. It further assumes that the contaminant concentration in the influent water is always zero as the water enters the contaminated region, but it adjusts instantly to a concentration in equilibrium with the remaining sorbed contaminant mass after the water enters the contaminated region.
The replacement of water occurring due to the P&T system is simulated through a series of consecutive discrete flushing periods. During each flushing period, enough clean water is introduced, at a known rate, to fill the pore space in a given volume of aquifer. Soil and water concentrations at each flushing period i can be calculated by the following equations, assuming instantaneous linear equilibrium between solid and liquid phase concentrations: where C s (i) is the soil total volatile organics concentration in (mg/kg) after i flushes, and C w is the concentration of total volatile organics in the water in equilibrium with the soil (mg/L). Once C w (i−1) is evaluated as the ratio between C s (i) and K d , the value for C w (i) can be entered into Equation (5) to calculate the soil concentration after the next flush. This is repeated until the soil and groundwater reach the desired concentrations.
Batch flushing described in [38] is an explicit finite difference approximation of the solution to the governing differential equation. As this approximation is relatively difficult to use and is subject to numerical errors, [7,39] suggest using the exact solution to the governing differential equation, which is the following expression for the number of flushes (or pore volumes) required to reach the clean-up concentration C f , starting from an initial contaminant concentration in groundwater of C 0 , with R as the retardation factor: From Equation (6), the relationship between the decrease in concentration as a function of the number of pore volume flushes (N PV ) is obtained: The total groundwater volume (m 3 ) to be extracted (V e ) to reach a concentration in groundwater lower than the clean-up level is: and since the duration of each flushing period (t) can be obtained dividing the contaminated volume PV by the discharge rate of the P&T system: PV can be calculated as the product of the P&T system discharge rate (Q) and the flushing duration (t) and then substituting it into Equation (8). Hence, N PV can be calculated as the ratio between V e and Q·t, which can be entered into Equation (7) to obtain: The total time required to reach remediation by a P&T system made up of n wells is then given by: In case the P&T system is not yet active and/or the attenuation effect on the contamination due to the natural groundwater flow rate is to be assessed, the volume of water that naturally flows through a given section of the contaminated area can be calculated as: where K x is the hydraulic conductivity (m/s), A is the contaminated area given by the product of the width of the contamination plume and the saturated thickness of the contaminated aquifer, and "i" is the hydraulic gradient near the contaminated area. The batch flushing model tends to underestimate clean-up time [7] for several reasons: the assumption of chemical equilibrium between sorbed and dissolved contaminant mass; the assumption of linear sorption/desorption processes; the hypothesis that the concentration of the contaminant in the water used to flush the aquifer is less than or equal to the desired clean-up level, and that regardless of concentration, this level remains constant during the entire flushing process; the fact that it does not account for processes due to heterogeneities present in the domain; and the contamination source is assumed exhausted from the moment in which the estimation of the clean-up time starts.

Advection-Dispersion-Retardation Model
Another analytical model that accounts for advection, dispersion, and retardation processes but neglects matrix diffusion is the advection-dispersion-retardation (ADR) model, described in [41]. The authors present an equation to simulate the flow of clean groundwater through a contaminated area and the resulting flushing of contamination within the volume of control, using a modified version from Ogata-Banks, in which longitudinal dispersivity equals 10% of the flow path length: where C(t) is the contaminant concentration at time t (µg/L), C 0 is the highest concentration observed during the first two years of monitoring data (µg/L), and their ratio represents the normalized contaminant concentration. Knowing the normalized contaminant concentration and the retardation factor, the number of pore volumes to be extracted (N PV ) can be determined by solving Equation (13) iteratively. Multiplying this value for the pore volume (PV), the total amount of groundwater to be pumped is obtained. From the latter, based on the P&T system discharge, the time to reach site remediation is derived. The ADR model can be applied at those sites where the governing transport processes do not imply the presence of diffusion. For this model, the same indications as for the batch flushing apply, with the exception that the ADR model also accounts for dispersive phenomena. The formulation of the ADR model presented by [41] assumes that the contaminant is present in the dissolved and sorbed phase only, and that no other sources are present.

Square Root Model
The third analytical approach selected refers to the models known as the square root model (SRM) and Dandy-Sale model (DSM), which allow for the assessment of matrix diffusion effects and to understand whether remediation limits can be reached in a reasonable timeframe.
As a matter of fact, in complex hydrogeological sites, where the aquifer is heterogeneous and fine materials (clay or silt) are present, or NAPL contamination is revealed, the remediation time becomes longer because of matrix diffusion processes. In this case, the remediation timeframe can be evaluated by means of these models, both implemented in the Matrix Diffusion Toolkit. The Toolkit was developed on the basis of previous works (among many [49]) by GSI Environmental Inc. in cooperation with the Colorado State University for ESTCP (Environmental Security Technology Certification Program) [48]. The Toolkit is freely available software with a Microsoft Office Excel GUI and can perform analyses with both models by evaluating the effects of matrix diffusion in the source area, in the plume, and downgradient of a permeable reactive barrier.
The Toolkit uses a simplified conceptual model of a two-layer aquifer system (a transmissive layer and a low permeability nontransmissive layer) where there are two different time periods: the "loading period", where there is a constant concentration of contaminants in the transmissive zone that drives contaminants into the low permeability zone, and the "release period", where the transmissive zone is assumed to have no concentration, and an upper-range estimate of release from the low permeability zone is generated. The low permeability zone is at least 1 m thick and there is no degradation in the low-k zone. The SRM assumes that the nontransmissive zone is "loaded" through a horizontal area placed exactly above it, as shown in Figure 4a, while the DSM assumes that the source area is a vertical plane (parameter L) and assesses the effects of the diffusion only downgradient of this plane, as shown in Figure 4b. Concentration results in the transmissive zone from both the SRM and DSM are based on estimates of mass discharge leaving the low permeability zone. Concentrations are calculated by assuming a 10-foot screened interval. Since in this paper the SRM is used, the following considerations refer to this model only. The Toolkit provides an instantaneous mass discharge from the entire area (A) during the release period. This mass discharge from the entire low permeability zone is assumed to be transported instantaneously to the downgradient edge of the modeled area (there is no advection or travel time component of the SRM).
At any time t, the concentration of contaminants in the transmissive zone can be estimated using the equation: where C(t) is the plume concentration in the transmissive zone at time t (μg/L), Md the mass discharge from the low permeability layer into the transmissive layer (μg/s), vD the Darcy velocity of the transmissive compartment (m/s), H the screened interval of the hypothetical well (m), and W the width of the modeled area (m). Under the previous hypotheses, at any time t, the mass discharge into the low permeability layer underlying the pool can be estimated using the equation: where t is the time after the source was introduced (year), t' the time at which source was removed from the high-permeability compartment (year), θ the porosity of low permeability zone (-), Cs the mean plume concentration above the low permeability compartment (μg/L), A the area of low permeability compartment beneath the plume (m 2 ), R the retardation factor for the low permeability compartment (-), and De the effective aqueous phase diffusion coefficient in the low permeability compartment (m 2 /s) evaluated as the product of θ (where p is the apparent tortuosity factor exponent) and D0 (molecular diffusion coefficient in free water, m 2 /s).
As seen, some of the input data are similar to those used for existing solute transport models, (e.g., Darcy groundwater velocity, size of the modeled area, information on when the source started) while other input data may appear particular (such as the tortuosity of the low permeability materials where matrix diffusion occurs); these are provided in Tables 1 and 2. Among the inputs, the mean concentration value at the interface between transmissive and nontransmissive zones is important, as well as the extension of the contaminated area. The length and width of the modeled area for the SRM can be entered directly or selected using a historical contour map, preferably one with the highest historical concentrations observed during the monitoring record (drawing a downgradient transect line perpendicular to groundwater flow and an upgradient transect line perpendicular to groundwater flow). Since in this paper the SRM is used, the following considerations refer to this model only.
The Toolkit provides an instantaneous mass discharge from the entire area (A) during the release period. This mass discharge from the entire low permeability zone is assumed to be transported instantaneously to the downgradient edge of the modeled area (there is no advection or travel time component of the SRM).
At any time t, the concentration of contaminants in the transmissive zone can be estimated using the equation: where C(t) is the plume concentration in the transmissive zone at time t (µg/L), M d the mass discharge from the low permeability layer into the transmissive layer (µg/s), v D the Darcy velocity of the transmissive compartment (m/s), H the screened interval of the hypothetical well (m), and W the width of the modeled area (m). Under the previous hypotheses, at any time t, the mass discharge into the low permeability layer underlying the pool can be estimated using the equation: where t is the time after the source was introduced (year), t' the time at which source was removed from the high-permeability compartment (year), θ the porosity of low permeability zone (-), C s the mean plume concentration above the low permeability compartment (µg/L), A the area of low permeability compartment beneath the plume (m 2 ), R the retardation factor for the low permeability compartment (-), and D e the effective aqueous phase diffusion coefficient in the low permeability compartment (m 2 /s) evaluated as the product of θ p (where p is the apparent tortuosity factor exponent) and D 0 (molecular diffusion coefficient in free water, m 2 /s).
As seen, some of the input data are similar to those used for existing solute transport models, (e.g., Darcy groundwater velocity, size of the modeled area, information on when the source started) while other input data may appear particular (such as the tortuosity of the low permeability materials where matrix diffusion occurs); these are provided in Tables 1 and 2. Among the inputs, the mean concentration value at the interface between transmissive and nontransmissive zones is important, as well as the extension of the contaminated area. The length and width of the modeled area for the SRM can be entered directly or selected using a historical contour map, preferably one with the highest historical concentrations observed during the monitoring record (drawing a downgradient transect line perpendicular to groundwater flow and an upgradient transect line perpendicular to groundwater flow).
For the field and synthetic cases, the simulations were run performing a "plume analysis", and the length and width of the polluted area were chosen as shown in Figure 5. The blue box is usually drawn around the highest concentration contour downgradient of the source area, while the green box is usually drawn around the second highest concentration contour downgradient of the source area. The width of this green box is fixed on the basis of the transverse distance of two points referring to the same isocontour line. Because the model estimates the results on the right edge of the external rectangle, the length is defined on the basis of the position of the target point. For the field and synthetic cases, the simulations were run performing a "plume analysis", and the length and width of the polluted area were chosen as shown in Figure 5. The blue box is usually drawn around the highest concentration contour downgradient of the source area, while the green box is usually drawn around the second highest concentration contour downgradient of the source area. The width of this green box is fixed on the basis of the transverse distance of two points referring to the same isocontour line. Because the model estimates the results on the right edge of the external rectangle, the length is defined on the basis of the position of the target point. Both the loading concentration value in the modeled source area and the time when the contamination source started and stopped are crucial parameters to obtain good results. Since the contamination event often occurs much earlier than the installation of the monitoring wells, the assessment of the loading concentration value may be very difficult. In [48], two approaches are suggested to evaluate the value: the first involves studying the history of the site and information about the processes (e.g., presence of Dense Not Aqueous Phase Liquids or other contaminants released); the second involves setting the maximum concentration measured over time as the load concentration value in the study area.
By means of the Toolkit, it is therefore possible to estimate both the concentration of contaminants at a monitoring point in the transmissive zone downgradient of the modeled area and the mean contaminant concentrations leaving the source area, as well as the concentration of contaminants at each point in the low permeability zone at different times. Moreover, the SRM allows for the evaluation of the mass discharge in the low permeability unit. Hence, in one monitoring point in the transmissive zone, it is possible to estimate the time to reach concentration values equal or lower than the remediation targets.
In the face of a relative ease of implementation and use, the SRM has some limitations related to the hydrogeological context or the type of contamination. For instance, the SRM evaluate the effects of only one contaminant at a time. Up until now the effect of matrix diffusion in a low transmissive zone has been studied for chlorinated solvents and Methyl Tertiary Butyl Ether (MTBE), although matrix diffusion processes can be valid for each dissolved contaminant, but the overall impact could be different [48]. Moreover, additional interpretation and expertise is required to apply the model in fractured rocks [48]. Finally, the system is assumed to be of the "on/off" type, with a "loading" period that is extended for a defined time and then instantly switches to a "release" period, when the concentrations in the transmissive unit, not resulting from back diffusion, are considered null.

Numerical Model
The flow and solute transport numerical model in the aquifer were implemented by means of Both the loading concentration value in the modeled source area and the time when the contamination source started and stopped are crucial parameters to obtain good results. Since the contamination event often occurs much earlier than the installation of the monitoring wells, the assessment of the loading concentration value may be very difficult. In [48], two approaches are suggested to evaluate the value: the first involves studying the history of the site and information about the processes (e.g., presence of Dense Not Aqueous Phase Liquids or other contaminants released); the second involves setting the maximum concentration measured over time as the load concentration value in the study area.
By means of the Toolkit, it is therefore possible to estimate both the concentration of contaminants at a monitoring point in the transmissive zone downgradient of the modeled area and the mean contaminant concentrations leaving the source area, as well as the concentration of contaminants at each point in the low permeability zone at different times. Moreover, the SRM allows for the evaluation of the mass discharge in the low permeability unit. Hence, in one monitoring point in the transmissive zone, it is possible to estimate the time to reach concentration values equal or lower than the remediation targets.
In the face of a relative ease of implementation and use, the SRM has some limitations related to the hydrogeological context or the type of contamination. For instance, the SRM evaluate the effects of only one contaminant at a time. Up until now the effect of matrix diffusion in a low transmissive zone has been studied for chlorinated solvents and Methyl Tertiary Butyl Ether (MTBE), although matrix diffusion processes can be valid for each dissolved contaminant, but the overall impact could be different [48]. Moreover, additional interpretation and expertise is required to apply the model in fractured rocks [48]. Finally, the system is assumed to be of the "on/off" type, with a "loading" period that is extended for a defined time and then instantly switches to a "release" period, when the concentrations in the transmissive unit, not resulting from back diffusion, are considered null.

Numerical Model
The flow and solute transport numerical model in the aquifer were implemented by means of the finite difference numerical code MODFLOW/MT3DMS, in order to replicate the four scenarios described in Section 2.2. In the hydrogeology field, MODFLOW/MT3DMS are the most widely used programs to investigate environmental problems related to forecasting quantitative and qualitative impacts on groundwater resources. MODFLOW-2000 (version 2000, United States Geological Survey, Reston, VA, USA) [40], developed in 1988 by the USGS (United States Geological Survey), is the most popular code for groundwater flow simulation in a three-dimensional porous medium. The partial differential equation that describes the three-dimensional transient groundwater flow in heterogeneous and anisotropic media, along the main directions (x, y, z), is written as: where K xx , K yy , and K zz are the hydraulic conductivities along the main directions (m/s), t the time (s), h the piezometric head (m), q s the volumetric flow rate per unit volume of aquifer representing sources and sinks ((m 3 /s)/m 3 ), and Ss the specific storage (m −1 ). MT3DMS (Model Transport in 3 Dimensional Multispecies) model is a modular three-dimensional multispecies transport model for simulation of advection, dispersion, and chemical reactions of contaminants in groundwater systems [50], and it represents the evolution of MT3D. MT3DMS can be used to simulate concentration changes of miscible contaminants in groundwater considering advection, dispersion, diffusion, and some basic chemical reactions, with various types of boundary conditions and external sources or sinks. It can also accommodate very general spatial discretization schemes and transport boundary conditions and is designed for use with any block-centered finite difference flow model.
The partial differential equation describing the fate and transport of contaminants of species k in three-dimensional, transient groundwater flow systems, disregarding chemical reactions, can be written as follows [50]: where C k is the dissolved mass concentration of species k (kg/m 3 ), D ij the diffusion-dispersion matrix (m 2 /s), C s is the concentration of the sources or sinks (kg/m 3 ). Several numerical models (Table 3) were developed by implementing properties and hydraulic conditions of the four scenarios defined in Section 2.2. For each numerical model, the modeling domain was equal to the study area defined in Section 2.2, the grid was discretized through 175 rows and 230 columns ( Figure 6), with square cell sizes ranging from 10 m in external area to 1 m in the central area where the clay lense was present. Twenty layers were implemented as vertical discretization in models named A, B, and C; the first 1.5 m thick and the subsequent ones 0.5 m thick (total thickness equal to 11 m); layers 8 to 20 were deleted (total thickness equal to 4.5 m) in model D. For each model, the first layer is partially saturated, unlike subsequent saturated layers. Hydrogeological and physical properties assigned to the domain are described in Table 2; sand properties were assigned to the transmissive zone cells (layers from 1 to 20), whereas clay lense properties were assigned to the nontransmissive zone cells (refined area in layers [8][9][10][11][12]. Initially, models A, B, C, and D were run in a transient regime 5 years long to reproduce the six "loading" phases of contamination in aquifer. Subsequently, numerical simulations of models A1, A2, B1, B2, C1, C2, D1, and D2 were run in a transient regime 200 years long to reproduce two different remediation scenarios; A1, B1, C1, and D1 simulated a remediation by natural attenuation (i.e., natural groundwater flow), whereas models A2, B2, C2, and D2 simulated a remediation by means of a P&T system.
For groundwater flow boundary conditions, two constant head boundary conditions (Dirichelet condition) were assigned at the west and east boundaries (values defined in Section 2.2).
Internal conditions required some distinction according to the model developed: -Model A: the internal specified mass flow condition was used defining a constant water flow with a specified concentration (Neumann condition) to reproduce the source of contamination; it was inserted in a specific area of the domain (36 m 2 ) from layer 1 to 7. The concentration value (317.1 mg/L) was introduced into the aquifer as a TCE mass equal to 500 kg for 5 years, with a low flow rate (0.01 L/s, as not to influence the hydraulic head distribution). A1 and A2 models: the internal specified mass flow condition was no longer implemented. In the A1 model, no additional internal conditions were added, whereas in the A2 model, the Neumann condition (through "analytic element well") was added, 140 m downgradient of the contamination source area, to reproduce the pumping operation of one well (flow rate equal to 40 m 3 /day). -Model B: the same internal condition as model A was implemented at the same cells of model A, but from layer 1 to 12. A constant concentration value (304.4 mg/L to reproduce 480 kg TCE mass) in the transmissive zone and initial concentration value (3703.7 mg/L to reproduce 20 kg TCE) in the nontransmissive zone were assumed. B1 and B2 models: the internal specified mass flow condition was no longer implemented (i.e., primary source instantly removed). In B1, no additional internal conditions were added, whereas in B2 the same Neumann condition as in model A2 was implemented to simulate the pumping well. -Model C: the same internal condition as model B to reproduce the contamination source (480 kg TCE mass); a constant concentration value (304.4 mg/L) was assigned at the same cells as in model A, from layer 1 to 7; another constant concentration value (18,518.5 mg/L) was implemented only in layer 10 on the same area (reproducing 20 kg TCE). C1 and C2 models: the internal specified mass flow condition was no longer implemented; in C1 no additional internal conditions were added, whereas in C2 the same Neumann condition as in A2 was implemented to simulate the pumping well. Hydrogeological and physical properties assigned to the domain are described in Table 2; sand properties were assigned to the transmissive zone cells (layers from 1 to 20), whereas clay lense properties were assigned to the nontransmissive zone cells (refined area in layers [8][9][10][11][12]. Initially, models A, B, C, and D were run in a transient regime 5 years long to reproduce the six "loading" phases of contamination in aquifer. Subsequently, numerical simulations of models A1, A2, B1, B2, C1, C2, D1, and D2 were run in a transient regime 200 years long to reproduce two different remediation scenarios; A1, B1, C1, and D1 simulated a remediation by natural attenuation (i.e., natural groundwater flow), whereas models A2, B2, C2, and D2 simulated a remediation by means of a P&T system.
For groundwater flow boundary conditions, two constant head boundary conditions (Dirichelet condition) were assigned at the west and east boundaries (values defined in Section 2.2).
Internal conditions required some distinction according to the model developed: -Model A: the internal specified mass flow condition was used defining a constant water flow with a specified concentration (Neumann condition) to reproduce the source of contamination; it was inserted in a specific area of the domain (36 m 2 ) from layer 1 to 7. The concentration value (317.1 mg/L) was introduced into the aquifer as a TCE mass equal to 500 kg for 5 years, with a low flow rate (0.01 L/s, as not to influence the hydraulic head distribution). A1 and A2 models: the internal specified mass flow condition was no longer implemented. In the A1 model, no additional internal conditions were added, whereas in the A2 model, the Neumann condition (through "analytic element well") was added, 140 m downgradient of the contamination source area, to reproduce the pumping operation of one well (flow rate equal to 40 m 3 /day). -Model B: the same internal condition as model A was implemented at the same cells of model A, but from layer 1 to 12. A constant concentration value (304.4 mg/L to reproduce 480 kg TCE mass) in the transmissive zone and initial concentration value (3703.7 mg/L to reproduce 20 kg TCE) in the nontransmissive zone were assumed. B1 and B2 models: the internal specified mass flow condition was no longer implemented (i.e., primary source instantly removed). In B1, no additional internal conditions were added, whereas in B2 the same Neumann condition as in model A2 was implemented to simulate the pumping well. A, from layer 1 to 7; another constant concentration value (18,518.5 mg/L) was implemented only in layer 10 on the same area (reproducing 20 kg TCE). C1 and C2 models: the internal specified mass flow condition was no longer implemented; in C1 no additional internal conditions were added, whereas in C2 the same Neumann condition as in A2 was implemented to simulate the pumping well. -Model D: same Neumann condition as model A applied to simulate the contamination source. D1 and D2 models: the internal specified mass flow condition was deleted; in D1 no additional internal conditions were added, whereas in D2 the same the Neumann condition (through "analytic element well") as in A2 was implemented to simulate the pumping well, but in this model the flow rate equaled 7.2 m 3 /day.

Field Case Results
Since no P&T system has been applied yet to the field case, the preliminary application of the analytical approaches discussed allowed for evaluation of the remediation time needed to clean-up the site and to recognize the possible advantage of the implementation of a P&T compared to the case where no measure is taken.

Batch Flushing and ADR Models
Disregarding the potential for matrix diffusion posed by the clay layer at the bottom of the aquifer (Figure 1), the time to clean up the site was evaluated by means of the batch flushing and ADR models. The width of the polluted area had a great influence on the results, therefore a detailed analysis of the data measured during several hydrochemical monitoring campaigns was carried out. On the basis of the available data and of the distribution of the monitoring points, the contaminated area was assessed to be 389 m wide and 885 m long. Using Equation (4), the pore volume value, which is representative of the volume of contaminated water, was calculated.
The 1,2-DCE initial concentration was the maximum recorded value in the piezometer p2 (about 20,000 µg/L), during the last 20 years. Knowing the retardation factor (Table 1) and the clean-up goal set by the Italian regulation (60 µg/L), the number of flushings was calculated by means of Equation (6). Subsequently, the total volume of water to be extracted was calculated by means of Equation (8) (Table 4). Knowing the size of the polluted area, the hydraulic conductivity of the aquifer (Table 1) and the hydraulic gradients, both the groundwater flow rates, and the remediation times were evaluated by means of Equations (11) and (12) for the cases with (Solution 2) and without (Solution 1) the P&T system (Table 4).
For the batch flushing model, about 7 years were needed to reduce 1,2-DCE from the initial concentration of 20,000 µg/L to the desired level of 60 µg/L. No significant difference was observed between natural flushing and the application of a P&T system, due to the negligible increase in the hydraulic gradient caused by pumping in the modeled area (about 900 m 2 ).
As for the ADR model, solving Equation (13) allowed for the identification of the value of the number of flushings (Table 5). Using the calculated N PV , the water volume needed to be flushed out (V e ) was evaluated and, finally, from the ratio between V e and groundwater flow rates, the remediation times were obtained.  The remediation time is reduced to 3.6 years using the ADR model. As seen before, both analytical approaches can be easily and quickly implemented and allowed for evaluating remediation times, even prior to the installation of a P&T system, requiring a minimum number of parameters.
Comparing the results of both models to real monitoring data show that they underestimate remediation times, because for batch flushing in 2013 and for the ADR model in 2010, the clean-up goal should have been reached, while p2 still had concentrations higher than the regulatory limit.
In any case, even assuming matrix diffusion is not a relevant phenomenon at the site, the application of the batch flushing and ADR models in this case showed that the installation of a P&T system would not cause any appreciable decrease in the remediation times. This occurs because the hydraulic gradient was assessed considering the whole length of the contaminated area (about 900 m) was not high and the piezometric drawdown, due to the high groundwater flow velocity, were not so relevant.

Square Root Model
A detailed examination of the stratigraphy in the field case shows that a clay layer is present at the bottom of the aquifer; therefore, it is important to evaluate how matrix diffusion could affect the remediation time. Hence, the SRM, as implemented in the Matrix Diffusion Toolkit, was used.
Following the data input sequence of the Toolkit (into the different Sections), in Sections 1 and 2 the reference system and the type of analysis ("plume analysis") were selected; in Sections 3 and 4 the values of hydraulic conductivity, total porosity, Darcy velocity, molecular diffusion coefficient, and retardation factor (Table 1) were assigned. Lacking detailed information on the nontransmissive zone, the same values of the synthetic case were assigned ( Table 2). In Section 5 the details of the polluted area were implemented: the highest historical concentration (equal to 20,000 µg/L) was assigned to an area of 150 × 150 m 2 and a secondary (optional) representative concentration (equal to 4634 µg/L) was assigned to a larger area (570 × 200 m 2 ). Four simulations were carried out considering piezometer p5 as the reference monitoring well, since very few and incoherent information was available on the source duration (loading period), being that the site was no longer in operation after 1990, but revealing the analysis of contaminant concentrations are still and active source. For this reason, in Section 6, two different loading periods were implemented; the first one (simulations SRM_1 and SRM_3 in Figure 7) involved a contamination source active for 30 years (from 1960 to 1990, i.e., the time period the chemical industry was active), the second one (simulations SRM_2 and SRM_4 in Figure 7) a contamination source was active for 53 years (from 1960 to 2013, i.e., assuming the source was not completely removed/contained and was still partially active even after manufacturing ceased). The simulations were run without (SRM_1 and SRM_2) and with (SRM_3 and SRM_4) a P&T system, which, as with the other models, slightly increased the hydraulic gradient (and therefore the groundwater velocity). The concentration distribution over time in the piezometer p5, 350 m from the source located in the proximity of p2, is shown in Figure 7, where simulated data (orange, blue, green and pink colors) are plotted together with the available detected concentrations in p5 (red points). Numerical simulations run with a loading period of 53 years in the absence (SRM_2) or in the presence (SRM_4) of a P&T system fit the measured concentrations very well, proving the scenario with an active source until 2013 as more realistic. For this scenario, about 28 years starting from now (2020) would be needed to reach aquifer restoration.
Compared to the batch flushing and ADR models, which neglect back diffusion, the remediation time increased by an order of magnitude. This case highlights the influence of matrix diffusion on remediation time and emphasizes how helpful a sound conceptual site model is in selecting the suitable approach. As the good fit between simulated and observed data in Figure 7 shows, the presence of clay at the aquifer bottom, in the case under analysis, plays an important role in trapping and then releasing contaminants due to back diffusion. For this site, this estimate is more sensible and truthful than the previous analytical approaches, as measured data show.

Benchmark Case Results
In the following section, the batch flushing, ADR, and SRM solutions are compared with a synthetic case study created by means of a numerical model (MODFLOW/MT3DMS). The hydrogeological structure was intentionally maintained simple in order to have the chance to easily compare the different solutions with the numerical one.

Numerical Modeling
The output of the numerical solute transport simulations of Model A, as well as of the other models, provided a concentration value in each domain cell. For instance, Figure 8 shows the polluted area after 5 years, at the end of the loading period. The polluted area, considering the isoconcentration line equal to 1000 mg/L, had a width of 60 m and a length of 170 m. Numerical simulations run with a loading period of 53 years in the absence (SRM_2) or in the presence (SRM_4) of a P&T system fit the measured concentrations very well, proving the scenario with an active source until 2013 as more realistic. For this scenario, about 28 years starting from now (2020) would be needed to reach aquifer restoration.
Compared to the batch flushing and ADR models, which neglect back diffusion, the remediation time increased by an order of magnitude. This case highlights the influence of matrix diffusion on remediation time and emphasizes how helpful a sound conceptual site model is in selecting the suitable approach. As the good fit between simulated and observed data in Figure 7 shows, the presence of clay at the aquifer bottom, in the case under analysis, plays an important role in trapping and then releasing contaminants due to back diffusion. For this site, this estimate is more sensible and truthful than the previous analytical approaches, as measured data show.

Benchmark Case Results
In the following section, the batch flushing, ADR, and SRM solutions are compared with a synthetic case study created by means of a numerical model (MODFLOW/MT3DMS). The hydrogeological structure was intentionally maintained simple in order to have the chance to easily compare the different solutions with the numerical one.

Numerical Modeling
The output of the numerical solute transport simulations of Model A, as well as of the other models, provided a concentration value in each domain cell. For instance, Figure 8 shows the polluted area after 5 years, at the end of the loading period. The polluted area, considering the isoconcentration line equal to 1000 mg/L, had a width of 60 m and a length of 170 m. In the following graphs, the results of the first (5 years) and second transient (200 years) numerical models were combined to show the complete concentration distribution at the end of the 205 years. More specifically, the attention is focused on the observation wells called OBS1 and OBS3, located 10 and 70 m downgradient from the contamination source, respectively.
In Figure 9a,b the concentration values of the first 12 layers of models A1 and A2 are shown for OBS1 and OBS3. In OBS1 well, the concentrations increased up to a maximum value, but, after source removal, they rapidly decreased by 50% over the following 3 years. In OBS3 well, farther than OBS1, the maximum concentration value was not as high and the time to reach it increased due to the propagation of the plume in the transmissive area.
A focus on the concentration distributions in the nontransmissive layers and vertical concentration distribution along an aquifer section is shown in Figure 10a,b. In the following graphs, the results of the first (5 years) and second transient (200 years) numerical models were combined to show the complete concentration distribution at the end of the 205 years. More specifically, the attention is focused on the observation wells called OBS1 and OBS3, located 10 and 70 m downgradient from the contamination source, respectively.
In Figure 9a,b the concentration values of the first 12 layers of models A1 and A2 are shown for OBS1 and OBS3. In the following graphs, the results of the first (5 years) and second transient (200 years) numerical models were combined to show the complete concentration distribution at the end of the 205 years. More specifically, the attention is focused on the observation wells called OBS1 and OBS3, located 10 and 70 m downgradient from the contamination source, respectively.
In Figure 9a,b the concentration values of the first 12 layers of models A1 and A2 are shown for OBS1 and OBS3. In OBS1 well, the concentrations increased up to a maximum value, but, after source removal, they rapidly decreased by 50% over the following 3 years. In OBS3 well, farther than OBS1, the maximum concentration value was not as high and the time to reach it increased due to the propagation of the plume in the transmissive area.
A focus on the concentration distributions in the nontransmissive layers and vertical concentration distribution along an aquifer section is shown in Figure 10a,b. In OBS1 well, the concentrations increased up to a maximum value, but, after source removal, they rapidly decreased by 50% over the following 3 years. In OBS3 well, farther than OBS1, the maximum concentration value was not as high and the time to reach it increased due to the propagation of the plume in the transmissive area.
A focus on the concentration distributions in the nontransmissive layers and vertical concentration distribution along an aquifer section is shown in Figure 10a In Figure 10a, layers 9 and 12 (even with the lowest concentration values) initially tend to store part of the TCE mass, but then they gradually release it because of back diffusion, when the concentrations in the transmissive layers tend to decrease. Layers 10 and 11, farther from the transmissive layers than layers 9 and 12, showed a similar behavior, but with lower concentration values, proving a continuous storing of the TCE mass. This was verified as shown in Figure 10b, observing the vertical distribution of the concentrations (2.5-8 m farther from the aquifer bottom) in OBS1 at different times. During the first years of simulation, the concentration values increased both in the sand and the clay because of the TCE propagation. Subsequently, the concentrations in the aquifer decreased rapidly, whereas the concentration peak moved into the low permeability layer, in accordance with [20], and took up the central portion of the lense with a value of 2500 μg/L, at the end of the 205 years.
The time when the TCE concentration values decreased below a target concentration value provides the remediation time of each model. In Italy, with respect to groundwater, two limits apply: the contamination threshold concentration (CTC) for groundwater remediation set by the legislative decree (D.Lgs. 152/2006) and the drinking water standard set by D.Lgs. 31/2001. For TCE, the CTC is 1.5 μg/L, whereas the drinking level is 10 μg/L (as a sum of TCE and PCE). The time to reach the remediation target is different at each point of the modeling domain. Considering 10 μg/L as the main target, according to model A1, in OBS3, more than 178 years were needed to achieve concentration values below the drinking water concentration; moreover, at the end of 205 years, the mean concentration value in the transmissive zone was equal to 6 μg/L in OBS1 and 8 μg/L in OBS3. These values are lower than the drinking water concentration, but still higher than the CTC (Table  6). In the A2 model, the concentration decrease phase was evident because of the effect of the pumping well; in OBS1, groundwater concentrations were reduced by 50% (about 1 year), compared to the maximum concentration value of the A1 model. In OBS3, the concentrations reached values lower than the drinking water level in 68 years, less than model A1, due to the effects of the pumping well (Table 6). However, the concentration values at the end of the 205 years were still higher than the CTC and this is a proof of the ineffectiveness of the P&T system due to the occurrence of the back diffusion and of the tailing effects caused by the presence of the lense. In Figure 10a, layers 9 and 12 (even with the lowest concentration values) initially tend to store part of the TCE mass, but then they gradually release it because of back diffusion, when the concentrations in the transmissive layers tend to decrease. Layers 10 and 11, farther from the transmissive layers than layers 9 and 12, showed a similar behavior, but with lower concentration values, proving a continuous storing of the TCE mass. This was verified as shown in Figure 10b, observing the vertical distribution of the concentrations (2.5-8 m farther from the aquifer bottom) in OBS1 at different times. During the first years of simulation, the concentration values increased both in the sand and the clay because of the TCE propagation. Subsequently, the concentrations in the aquifer decreased rapidly, whereas the concentration peak moved into the low permeability layer, in accordance with [20], and took up the central portion of the lense with a value of 2500 µg/L, at the end of the 205 years.
The time when the TCE concentration values decreased below a target concentration value provides the remediation time of each model. In Italy, with respect to groundwater, two limits apply: the contamination threshold concentration (CTC) for groundwater remediation set by the legislative decree (D.Lgs. 152/2006) and the drinking water standard set by D.Lgs. 31/2001. For TCE, the CTC is 1.5 µg/L, whereas the drinking level is 10 µg/L (as a sum of TCE and PCE). The time to reach the remediation target is different at each point of the modeling domain. Considering 10 µg/L as the main target, according to model A1, in OBS3, more than 178 years were needed to achieve concentration values below the drinking water concentration; moreover, at the end of 205 years, the mean concentration value in the transmissive zone was equal to 6 µg/L in OBS1 and 8 µg/L in OBS3. These values are lower than the drinking water concentration, but still higher than the CTC (Table 6). In the A2 model, the concentration decrease phase was evident because of the effect of the pumping well; in OBS1, groundwater concentrations were reduced by 50% (about 1 year), compared to the maximum concentration value of the A1 model. In OBS3, the concentrations reached values lower than the drinking water level in 68 years, less than model A1, due to the effects of the pumping well (Table 6). However, the concentration values at the end of the 205 years were still higher than the CTC and this is a proof of the ineffectiveness of the P&T system due to the occurrence of the back diffusion and of the tailing effects caused by the presence of the lense.
The results of the B1 model showed concentration values significantly higher in the proximity of the clay lense; in the nontransmissive layer, on the other hand, a rapid decrease in concentrations over time was obtained and the maximum concentration value was located in the middle of the clay lense. Even in OBS3 the concentrations were very high and therefore above the drinking water concentration limit; at the end of the 205 years, an average concentration value in the aquifer of 110 µg/L was achieved. Similar behavior occurred in model B2 simulations; the influence of the pumping well caused both a faster arrival of the peak concentration and a subsequent decrease in the aquifer concentrations. At the end of the 205 years, however, the average concentration value is higher (44 µg/L) than the drinking water concentration ( Table 6). The simulations of the C1 and C2 models showed similar results to the B1 and B2 models, with a specific slow propagation of the contaminant towards the central layers of the lense (9 and 11), and higher concentration values than in the B1 and B2 models. Even in the C1 and C2 models, at the end of the 205 years the average concentrations in OBS3 were 215 and 84 µg/L (Table 6), respectively, therefore much higher than the drinking water concentration.
Different results were achieved in the D1 and D2 simulations (homogenous aquifer), where the time to reach values lower than the drinking water concentration was shorter (40 years), since the clay lense was absent. Moreover, the average concentrations in OBS1 and OBS3, at the end of 205 years, were significantly lower than the CTC values.
In summary, the TCE mass stored in the nontransmissive layer plays an essential role in extending the remediation times. Moreover, as proved by the results of models C1 and C2, the position of the contaminant inside the lense also had a great influence on the long-term concentration values. Overall, the pumping well supports the flushing of groundwater in the polluted area, but results are insufficient to reduce the concentration values below the regulatory limits in a reasonable time frame, where fine materials are present. According to the results of models A, B, and C, the failure to achieve remediation targets in a short time is caused by the effects of back diffusion and probably also by the simple configuration of the P&T system simulated.

Batch Flushing and ADR Models
The remediation time was also evaluated by means of the batch flushing model to be compared to the results of numerical models D1 and D2. The polluted area, 65 m plume length and 50 m plume width, was defined on the basis of the numerical simulations results, this being a synthetic case without any hydrochemical monitoring campaign. The number of pore volumes was calculated by means of Equation (4).
The initial concentration was obtained assuming 100% TCE mass released in the polluted area; the ratio between TCE mass and polluted aquifer volume provided the value for the initial concentration. Knowing the retardation factor (Table 2) and the TCE drinking water concentration (10 µg/L), the number of flushings were calculated by means of Equation (6). Subsequently, the total volume of water needed to be extracted was calculated by means of Equation (8) ( Table 7). Knowing the extension of the polluted area, the hydraulic conductivity of the aquifer (Table 2), and the hydraulic gradients, both the groundwater flow rates and the remediation times were evaluated by means of Equations (11) and (12) for the cases with (Model D2) and without (Model D1) the P&T system ( Table 7).
The results of the batch flushing model, which accounts for advection and sorption only, were in accordance with D1 and D2 model results; depending on the chosen strategy (P&T or no active measure), the remediation times were always lower than 40 years.
Subsequently, the ADR model was implemented solving Equation (13) and identifying the value of the number of flushings (Table 8). Using the calculated N PV , the water volume needed to be flushed out (V e ) was evaluated and finally, from the ratio between V e and groundwater flow rates, the remediation times were obtained. As in the field case, the remediation times of the ADR model were lower than those of the batch flushing model, reducing from 38.6 to 15.4 years for the natural groundwater flow case and from 30.9 to 12.3 years for the P&T case. The low remediation times provided by the ADR model for this specific case could lead to consider the P&T a sustainable option; however, the difference between P&T and no active measure is not so remarkable as to implement such a system, unless a hydraulic containment is needed.

Square Root Model and Numerical Model
By means of the analytical approach SRM (implemented in the Matrix Diffusion Toolkit), which takes into account matrix diffusion, an evaluation of the remediation times considering the influence of fine materials was carried out and then compared with the results of models A1 and A2. Four different simulations were run in order to reproduce the remediation of the polluted aquifer by means of natural groundwater flow and P&T system in OBS1 (a distance of 10 m from the contamination source) and in OBS3 (a distance of 70 m).
Following the data input sequence of the Toolkit, in Sections 1 and 2 the reference system and the type of analysis ("plume analysis") were selected; in Sections 3 and 4 the values of hydraulic conductivity, total porosity, Darcy velocity, molecular diffusion coefficient, and retardation factor ( Table 2) were assigned. In Section 5, from the results of models A1 and A2, both concentration values (set equal to the mean value of the concentrations observed in OBS1 and in OBS3 over time) and the plume length and width (corresponding to those observed in OBS1 and OBS3) were input. When simulating the behavior of OBS1, a concentration value of 105,000 µg/L, a plume source size of 6 × 6 m 2 and a plume size of 10 × 10 m 2 were used. For the simulation of OBS3, a concentration value of 25,000 µg/L, a plume source size of 20 × 20 m 2 , and a plume size of 40 × 70 m 2 were used. For each simulation, the loading period was set equal to the contamination event duration, which was 5 years.
The distribution of the concentrations obtained from the analytical (SRM) and numerical (A1) models for the natural groundwater flow case are shown in Figure 11a,b, whereas those of the analytical (SRM) and numerical (A2) models for the P&T case are shown in Figure 11c,d. Water 2020, 12, x FOR PEER REVIEW 21 of 25 Figure 11. P&T vs. natural attenuation. Concentration distribution over time for SRM and numerical models in OBS1, cases (a,c); and in OBS3, cases (b,d).
As observed in Figure 11, in OBS1 146 years and 87 years were needed to reach concentrations below the drinking water level (10 μg/L) by means of natural groundwater flow or P&T remediation, respectively. In OBS3, remediation times were longer: 170 and 103 years to achieve concentration values lower than the target limit, as observed for A1 and A2 models in Table 6. For each simulated case, more than 200 years were needed to have concentrations lower than the regulatory limit (1.5 μg/L), as observed in the results of the numerical modeling.
Moreover, the SRM allowed for the evaluation of the contaminant mass value reached in the nontransmissive layer at the end of the loading period: 3.5 kg and 15 kg were present in OBS1 and OBS3, respectively, in accordance with the total TCE mass injected in the clay lense of the modeling domain.
The concentration distribution over time was driven by TCE mass release from fine material layers because of back diffusion, which was the main reason for the extended remediation times in both scenarios, with respect to the cases where no low permeability layers were present. In addition, the SRM approach suggested that the P&T installation was not advisable because of the very long estimated remediation times and consequent huge costs.
Compared to the batch flushing and ADR models, which neglect back diffusion, the remediation time increased by an order of magnitude, proving the great influence of the fine materials on the results. Therefore, the remediation time evaluated by means of SRM is more realistic when nontransmissive layers are present, as is typical for many aquifers.
The good agreement between the numerical and analytical results highlights the effectiveness of the latter in predicting remediation times, provided that the right model is used in accordance with the hydrogeological setting.

Conclusions
The P&T system is a very widespread remediation technology because of its simplicity and versatility. However, these systems have some limitations, which frequently do not allow for the achievement of remediation targets in reasonable times. Therefore, a preliminary methodology allows for an estimate of the remediation time of P&T systems and to understand if the planned strategy is feasible or an alternative methodology is needed, or if it is better to couple the P&T with a source control technology to enhance remediation. Indeed, often the hydrogeological setting (e.g., Figure 11. P&T vs. natural attenuation. Concentration distribution over time for SRM and numerical models in OBS1, cases (a,c); and in OBS3, cases (b,d).
As observed in Figure 11, in OBS1 146 years and 87 years were needed to reach concentrations below the drinking water level (10 µg/L) by means of natural groundwater flow or P&T remediation, respectively. In OBS3, remediation times were longer: 170 and 103 years to achieve concentration values lower than the target limit, as observed for A1 and A2 models in Table 6. For each simulated case, more than 200 years were needed to have concentrations lower than the regulatory limit (1.5 µg/L), as observed in the results of the numerical modeling.
Moreover, the SRM allowed for the evaluation of the contaminant mass value reached in the nontransmissive layer at the end of the loading period: 3.5 kg and 15 kg were present in OBS1 and OBS3, respectively, in accordance with the total TCE mass injected in the clay lense of the modeling domain.
The concentration distribution over time was driven by TCE mass release from fine material layers because of back diffusion, which was the main reason for the extended remediation times in both scenarios, with respect to the cases where no low permeability layers were present. In addition, the SRM approach suggested that the P&T installation was not advisable because of the very long estimated remediation times and consequent huge costs.
Compared to the batch flushing and ADR models, which neglect back diffusion, the remediation time increased by an order of magnitude, proving the great influence of the fine materials on the results. Therefore, the remediation time evaluated by means of SRM is more realistic when nontransmissive layers are present, as is typical for many aquifers.
The good agreement between the numerical and analytical results highlights the effectiveness of the latter in predicting remediation times, provided that the right model is used in accordance with the hydrogeological setting.

Conclusions
The P&T system is a very widespread remediation technology because of its simplicity and versatility. However, these systems have some limitations, which frequently do not allow for the achievement of remediation targets in reasonable times. Therefore, a preliminary methodology allows for an estimate of the remediation time of P&T systems and to understand if the planned strategy is feasible or an alternative methodology is needed, or if it is better to couple the P&T with a source control technology to enhance remediation. Indeed, often the hydrogeological setting (e.g., presence of areas with low permeability) or the type of contaminants (e.g., NAPL) lead to a prolongation of the remediation times of a P&T system.
In the present work, some easy to apply analytical approaches were used to evaluate the remediation time of a P&T system for a synthetic case study; the batch flushing (considering advection and sorption terms), ADR (considering advection sorption and dispersion terms), and SRM (also considering matrix diffusion) models were compared with a numerical solution. These analytical models were applied to a synthetic case to show how they are able to provide truthful results if applied in accordance with the hydrogeological model, without the need to implement a specific numerical model.
This study represents a first result showing that in the case of homogeneous and hydraulically transmissive domains, the use of the batch flushing and ADR models provide realistic evaluation of the remediation times due to the P&T system. In contrast, for cases where low permeability zones are present and matrix diffusion plays a key role, more complex models, such as SRM, need to be adopted.
The analytical solutions were also applied at a field scale for 1,2-DCE contamination under natural attenuation and P&T conditions. Batch flushing and ADR model underestimated remediation times (less than 7 and 4 years, respectively) when compared to the results of actual monitoring data, whereas SRM provided remediation times (less than 30 years) that were more truthful. Moreover, for this case where a contamination event is ongoing, the installation of a P&T system would not cause any appreciable decrease in the remediation times. That means P&T can be a good solution for plume containment, but is expected to be unable to improve the clean-up time if not combined with other remediation techniques.
In the synthetic case, several numerical solute transport models were implemented in MODFLOW/MT3DMS in order to simulate the effects of natural groundwater flow and a P&T system, with and without the presence of hydrogeological heterogeneities (i.e., low permeability zones). The results of the numerical model, in terms of time to reach a target concentration value, were compared to those of the analytical solutions, thus assessing the reliability of the analytical methods. In particular, the influence of back diffusion, which is among the main reasons P&T systems often come out as ineffective and unsustainable over time, was carefully analyzed to understand whether the SRM (easier and less expensive and time consuming than a numerical model) could provide reasonable outcomes.
When the aquifer is homogeneous, the batch flushing model results were comparable to the numerical models: 38 vs. 39.5 years (natural groundwater flow) and 31.6 vs. 35 years (P&T remediation). When hydrogeological heterogeneities play an important role, this method fails to provide a reasonable estimate of the remediation time, while the SRM results are much more realistic compared to the numerical models: 170 vs. 178 years (natural groundwater flow) and 103 vs. 110 years (P&T remediation).
This study was necessary to assess the advantages or disadvantages of the different approaches. It should be seen as the first step to understand whether a specific analytical approach could be applied to a "higher level" (laboratory scale or field case). In the future, the SRM should be thoroughly compared with a real case where a P&T system is active, in order to better assess the real applicability of this method.
Concluding, when P&T is active, concentrations in groundwater decrease faster than with natural groundwater flow; despite this aspect, the presence of hydraulically nontransmissive layers did not allow for the attainment of concentration targets in a reasonable time. This result suggests that concentration values in the transmissive zone are mainly influenced by mass contaminant release from the low permeability layers because of back diffusion. This phenomenon cannot be simulated by means of the batch flushing or ADR models; therefore, if based on a correct definition of the initial parameters, SRM can be considered a valid and easy to implement an analytical approach to evaluate remediation times whenever a contaminant has sufficient time to penetrate into low permeability layers.