A Comparative Assessment of Analytical Fate and Transport Models of Organic Contaminants in Unsaturated Soils

: Analytical models for the simulation of contaminants’ fate and transport in the unsaturated zone are used in many engineering applications concerning groundwater resource management and risk assessment. As a consequence, several scientiﬁc studies dealing with the development and application of analytical solutions have been carried out. Six models have been selected and compared based on common characteristics to identify pros and cons as well as to highlight any di ﬀ erence in the ﬁnal output. The analyzed models have been clustered into three groups according to the assumptions on contaminant source and physico-chemical mechanisms occurring during the transport. Comparative simulations were carried out with ﬁve target contaminants (Benzene, Benzo(a)pyrene, Vinyl Chloride, Trichloroethylene and Aldrin) with di ﬀ erent decay’s coe ﬃ cient, three types of soil (sand, loam and clay) and three di ﬀ erent thicknesses of the contaminant source. The calculated concentration at a given depth in the soil for the same contamination scenario varied greatly among the models. A signiﬁcant variability of the concentrations was shown due to the variation of contaminant and soil characteristics. As a general ﬁnding, the more advanced is the model, the lower the predicted concentrations; thus, models that are too simpliﬁed could lead to outcomes of some orders of magnitude greater than the advanced one. M.C., G.V. and C.G.; G.V.; curation, M.G.S.; and M.G.S.; and Author Contributions: Conceptualization, M.G.S., G.L., M.C., G.V., C.G., and L.C.; methodology, M.G.S., G.L. and M.C.; formal analysis, M.G.S., G.L. and M.C.; investigation, M.G.S., G.L., M.C., G.V. and C.G.; resources, G.V.; data curation, M.G.S.; writing—original draft preparation, M.G.S., G.L., M.C. and G.V.; visualization, M.G.S.; supervision, G.L., M.C., G.V., C.G., and L.C. All authors have read and agreed to the published version of the manuscript.


Introduction
The release of organic contaminants into the environment represents a growing threat for groundwater and human health, as evidenced by the concentrations measured in several groundwater bodies [1][2][3][4][5]. Groundwater contamination is a significant environmental issue because it can damage potable water supply, endangering human health [6,7], and can affect aquatic ecosystems [8,9]. Organic contaminants moving in groundwater come from both point sources (solid waste tips, landfills, leaking storage tanks and leaking sewers) and diffuse sources (agricultural activities and farmyard drainage) [10]. Most of the fate and transport processes take place into the unsaturated zone [11], commonly known as vadose zone, posing a great concern in terms of environmental preservation and protection [12]. Table 1. Analytical solutions of the selected models.

Group Model
Boundary and Initial Conditions Solution Reference ASTM, 2000 [42] II Enfield et al., 1982 [46] VI C w (z, 0) = 0 C w (0, t) = C w0 exp(−µ 3 t) 2001 [47] They share the following common features: (i) the assessment of the concentration of a single contaminant dissolved in aqueous phase that moves from a source located in the unsaturated zone to the water table; (ii) the negligible effect related to the movement of substances in the non-aqueous phase; (iii) the small number of geometrical, hydrogeological and chemical input data; (iv) the use of a set of simplifying assumptions that allows to obtain analytical solutions. These characteristics make them straightforward tools in many applications such as the analysis of groundwater vulnerability to pesticides contamination as well as the risk assessment of contaminated sites [48][49][50][51]. They differ mainly in their assumptions on the source and in the physico-chemical mechanisms as highlighted in Table 2.
The reference spatial domain assumed by the models is depicted in Figure 1. It consists of an unsaturated layer, lower bounded by an aquifer. The contaminant source, located in the upper region, is a parallelepiped shape with edges L 1 (thickness), X and Y. Upper and lower positions of the source along the z-axis are given by L f and L 2 , respectively, with reference to the water table and by L d with reference to the ground surface. The unsaturated zone is characterized by fixed homogeneous and isotropic hydrogeological properties θ, θ w , θ a , ρ b , K sat , f oc , I eff , according to the hypotheses of the models. The physical-chemical characteristics of the contaminants are indicated by H, D a 0 , D w 0 , S, K oc . The selected models, except Model I, consider the degradation processes of the contaminant in water phase using a first-order kinetic formulation, where λ is the degradation constant.
Sustainability 2020, 12, x FOR PEER REVIEW 2 of 26  The total concentration of the contaminant is expressed by Equation (1): whereas the concentration of the contaminant in the different phases (water; solid and air) is expressed by Equation (2): The mass of the contaminant sorbed by the soil fraction upon the mass of the affected soil C s and the mass of contaminant in vapor phase upon the volume of pore air C a are given by Equations (3) and (4), respectively: In all models, the partitioning of the chemicals among these three different phases is expressed by instantaneous linear equilibrium.
When low solute concentration occurs, a linear adsorption isotherm can be assumed [52]. Thus, the ratio between C s and C w is constant and equal to K d : where The liquid-vapor partition is evaluated through Henry's Law: On the basis of Equation (7), the dissolved concentration in pore water C w can be evaluated from C tot : where K sw is expressed by Equation (9) representing the dissolution of soil contaminants into infiltrating rainwater, according to the hypothesis of instantaneous and linear equilibrium. The maximum value of C w is the solubility (S) of the contaminant. Therefore, taking into account air, water and solid phases, the maximum C tot is Some of the selected models (Models III, IV, VI; Table 1) take into account the presence of a residual NAPL (Non-Aqueous Phase Liquid) in the source. This allows to consider an initial C tot higher than C tot;sat . For these models, when the measured C tot is higher than C tot;sat it is possible to consider C tot as a sum of two components: C tot = C tot;sat +C NAPL (11) where C NAPL is the concentration of the contaminant in non-aqueous phase. The initial dissolved concentration in pore water at the source C w0 can be evaluated as follows: C w0 = K SW ·C tot in case of C tot < C tot;sat (12) C w0 = S otherwise (13) In all the analyses here developed, C tot is assumed less than C tot;sat in order to use the same initial C tot for all the simulations.
Based on the common characteristics, the models were organized into three groups. The first group is composed by models with a continuous source (Models I and II), leading to a constant concentration after a certain period of time. The second group considers a decaying source. In this case, the transport equation is not solved (Models III), or it is solved for a steady state condition (Model IV). Finally, the third group considers a decaying source and a solution of the transport equation in transient conditions (Models V and VI). Specific features of the models will be carefully described in the following sections.

Model I
Model I was developed by ASTM (American Society for Testing and Materials) and described in Appendix 3 of the Standard Guide for Risk-Based Corrective Action [42]. It provides a very basic approach to estimate the dissolved concentration of contaminant in the water table. The model is based on the determination of the contaminant LF, which represents the steady ratio between the contaminant concentration in the soil source and the resultant concentration in groundwater. This approach may be applied both to organic and inorganic contaminants and it is adopted by various software for risk assessment, e.g., Risk-Based Corrective Action (RBCA) Tool Kit for Chemical Releases [53]. Its widespread use is related to its simplicity and a limited number of input parameters. The source is assumed uniform and constant, the only considered transport mechanism is the steady leaching from the source to the groundwater table.
The LF is composed by three coefficients: K SW (see Equation (9)) SAM and LDF. SAM, introduced by [53], represents the attenuation of contaminant concentration caused by its sorption related to the leaching through the clean soil. This coefficient is derived from the principle of mass conservation and it depends only on the geometry of the scheme. In fact, it is defined as the source thickness L 1 and the distance between the top of the source and the water table L 2 , ratio. Thus, it is possible to obtain C w at the water table as reported in Table 1-Model I. LDF evaluates the uniform dilution of the contaminant in groundwater, but this mechanism has not been taken into account.

Model II
Model II differs from the others as it considers combined gas and water phase diffusion in three dimensions [43]. For this reason, the model is particularly suitable to be applied when gas phase transport may be prevalent as for VOCs. It assumes mono-dimensional advection in water phase and three-dimensional diffusion in gas and water phases. In addition, it considers uniform and constant concentration of the contaminant in the source, neglecting the attenuation processes. Furthermore, soil sorption and first-order degradation take place in the unsaturated zone.
The starting point of the transport equations for water and gas phases is where R is assumed as follows: Sustainability 2020, 12, 2949 7 of 24 Troldborg and co-workers [44] obtained the following Equation (16) where where The dispersion coefficients of the diagonal tensor D', appearing in the matrix representation (18), are given by where D w e and D a e are estimated according to [54] whereas v is calculated by means of the following equation: The set of initial and boundary conditions is shown in Table 1-Model II. The solution of Equation (16) was developed by Sagar [55] and Wexler [56] and has been reported in Table 1-Model II as well.
Troldborg [43] pointed out that the evaluation of the contaminant dilution in groundwater is more difficult using a three-dimensional model because the contaminated area at the water table is bigger than the source area. Therefore, the water dilution volume has to be evaluated since the concentrations are not spatially constant within this area.

Model III
Model III consists of a subset module of a large model named CatchRisk [44]. This model evaluates the contaminant transport from a point source to the groundwater. In Model III, the unsaturated zone is designed like a reactor where the mass flux is governed only by water advection and the other transport mechanisms are neglected; the transport equation is not solved as in the Model I. On the contrary, Model III assumes that the dissolved concentration is constant along the depth.
The decay of the source concentration is taken into account. The mass of contaminant within the compartment decreases due to both leaching through the unsaturated zone and degradation of the contaminant in water phase. The dissolved concentration over the time can be obtained solving a mass balance. The variation of the chemical mass over the time is equal to the rate of mass depletion due to leaching and degradation. The solution of this mass balance is shown in Table 1-Model III. Furthermore, the model considers the residual non-aqueous phase concentration in the source. If the mass of contaminant is higher than the equilibrium mass, a separate phase is present. At the initial time, C tot and C w0 are described by Equations (11) and (13), respectively, and the total contaminant mass is equal to Sustainability 2020, 12, 2949 8 of 24 The mass flux J, Equation (25), leaving the unsaturated zone is constant on time until the contaminant mass is equal to the mass equivalent to the saturation conditions M sat , Equation (26); When M tot is equal to M sat , C w can be calculated through the solution shown in Table 1-Model III.

Model IV
Model IV is composed of three terms: α dep (t), α leach and LDF [45]. α dep (t) considers the source depletion caused by leaching and biodegradation. α leach describes the transport phenomena taking into account the mechanisms of advection, dispersion-diffusion, sorption and biodegradation. LDF evaluates the dilution of the contaminant in groundwater; however, this mechanism has not been taken into account.
The term α dep (t) is defined by an exponential law, obtained from a mass balance that includes the mass losses caused by leaching and biodegradation, reported in Table 1-Model IV. Furthermore, the model includes the presence of a residual phase in the source. In this case, the dissolved concentration in the source over the time is equal to and t* is the time when the initial source concentration reaches the saturation conditions, defined as The term α leach is obtained solving the advection-diffusion differential equation for the steady state: The dispersion-diffusion coefficient D z is calculated as where D w e is calculated as: Unlike other models, the seepage velocity is not obtained from the effective infiltration, but from the time required by the contaminant to reach the underlying water table t leach . The latter is calculated as a linear function of the time required to the infiltrating water to reach the water table t w , by means of the retardation coefficient of the contaminant R: Sustainability 2020, 12, 2949 9 of 24 R is calculated according to the following equation: The time t w is calculated using the Green and Ampt equation [57]: where H w is the ponding depth of water surface and h cr the wetting front suction head. The set of initial and boundary conditions to be given to the differential equation (30) is reported in Table 1-Model IV.

Model V
PESTAN (PESTicide ANalytical) is a software developed by the U.S. EPA (Environmental Protection Agency of United States) to assess the transport of organic contaminants, mainly pesticides [58]. The theoretical model underlying the software PESTAN was published in [46]. It differs from the others mainly due to the conceptualization of the source. The source is considered a mass of granular solid contaminant located on the ground surface which dissolves in a "slug" of contaminated water infiltrating into the soil. The thickness of the slug z 0 is the equivalent depth of the pore water required to dissolve the total mass of the solid chemical in water.
In order to compare the model with the others, the reference scheme presented above was adopted and z 0 was assumed equal to the equivalent depth of the pore water contained in the source (Equation (36)).
The differential equation underlying the model (Equation (37)) takes into account the following mechanisms: advection and dispersion in the vertical dimension, sorption and biochemical degradation described according to a first order kinetic.
where R is given in Equation (34). The term D is evaluated through the relationship of Biggar and Nielsen [59] D = D e w +2.93v 1.11 (38) where D and D w e are expressed in cm 2 day −1 and v in cm day −1 .
The initial and boundary condition for the Equation (37) and the proposed analytical solution are given in Table 1-Model V.

Model VI
Model VI was first implemented in the risk assessment software RISC 4 [47], then upgraded to the newer version RISC 5 , which incorporates the approach developed in [60]. In this model, the source depletes with time due to the simultaneous effect of leaching and volatilization of VOCs. Unlike the previous model, biological degradation is neglected. The presence of a residual phase in the source is considered as described above. The leachate concentration in the source decays exponentially with time. The exponential law is obtained from a mass balance where the variation over the time of the contaminant mass is equal to the rate of mass depletion due to leaching and volatilization. The source depletion law is The coefficient µ 3 is described by the following equation: where D eff is the effective diffusion coefficient in soil estimated using the Millington-Quirk relationships [61] and L d is the distance from the ground surface to the center of the source. Model VI is ruled by the same differential equation of Model V (Equation (37)). The dispersion coefficient is assumed as a linear function of the seepage velocity (Equation (42)) because the mechanism of molecular diffusion is neglected, where α l is calculated using the Gelhar relationship [62]: ln α l = −2.727 + 0.584 ln z m otherwise (44) being z m the distance from the source to the observation location.
The set of initial and boundary condition and the solution to be given are reported in Table 1-Model VI.

Methodology
The analysis has been carried out by comparing the solutions of the models with reference to the same scenario ( Figure 1) and varying soils, contaminants and thickness of the source L 1 . As one can expect, the models do not respond the same way to the variation of the parameters because of the different assumptions about the source and transport mechanisms.
The values of the geometrical parameters are reported in Table 3. Table 3. Values of the geometrical parameters.

Parameter * Value
Three types of soils, differing mainly in water content and hydraulic conductivity, have been chosen: sand, loam and clay. Furthermore, five organic contaminants have been selected according to their different mobility and biodegradability characteristics in order to explore a wide range of fate and transport mechanisms.
The soil parameter values are derived from [63,64] and listed in Table 4. Chemical properties of the contaminants have been derived by USEPA 1996 [65], except for the first order degradation constant λ, taken from the handbook of Howard [66]. The C w0 has been chosen according to the literature focusing on contaminated sites [43,[67][68][69][70]. The physico-chemical properties of the selected contaminants are reported in Table 5.

Predicted Concentration over the Time According to the Different Models
In order to evaluate the differences among the solutions of Models I-VI, the dimensionless concentrations C w /C w0 have been calculated at a certain spatial location (z = 5 m) as a function of time, fixing a single contaminant, Trichloroethylene (TCE), soil type (sand), source thickness (L 1 = 2 m). The results of the ratio between C w /C w0 are sketched in Figure 2.
As shown in Figure 2, the value of the initial dimensionless concentration and the shape of the curve are influenced by source hypotheses and transport mechanisms. According to these results, models belonging to group 1 (Model I and II) show a constant concentration over the time (with the exception of an initial transient phase for Model II); models belonging to group 2 (Model III and IV) show a decreasing exponential trend whereas models belonging to group 3 (Model V and VI) exhibit a temporal trend represented by a bell-shaped curve.
Model I is the most basic, as it neglects fate and transport mechanisms over time. Therefore, the dimensionless concentration C w /C w0 is constant. Conversely, Model II solves the transport equation, and it considers a constant source; therefore, at the initial time t = 0 the value of C w /C w0 is zero and after a transient phase reaches a constant value. Although both models showed a similar trend over the time, Model II depends on the spatial coordinates x, y, and z. The models belonging to group 1 can be used when the reduction processes of the source, like leaching and volatilization losses or biodegradation, are negligible, so the source can be modelled as constant. Pascuzzi et al. [71] used the model I for evaluating the risk for human health assuming a steady redeposition on the soil of pollutants in the selected site. The C W /C W0 was equal to 3.33 × 10 −2 , while according to our results the C W /C W0 calculated through the application of the SAM coefficient ranges between 0.167 and 0.444. These different values are related to the differences in the adopted reference schemes, which cause a stronger attenuation process in the scheme presented by Pascuzzi et al. They consider a source thickness L 1 of 0.2 m and a distance between the bottom of the source and the water table L f equal to 6m, while in our case L 1 ranges between 1 m and 4 m and L f is 5 m.

Predicted Concentration over the Time According to the Different Models
In order to evaluate the differences among the solutions of Models I-VI, the dimensionless concentrations Cw/Cw0 have been calculated at a certain spatial location (z = 5 m) as a function of time, fixing a single contaminant, Trichloroethylene (TCE), soil type (sand), source thickness (L1 = 2 m). The results of the ratio between Cw/Cw0 are sketched in Figure 2.
As shown in Figure 2, the value of the initial dimensionless concentration and the shape of the curve are influenced by source hypotheses and transport mechanisms. According to these results, models belonging to group 1 (Model I and II) show a constant concentration over the time (with the exception of an initial transient phase for Model II); models belonging to group 2 (Model III and IV) show a decreasing exponential trend whereas models belonging to group 3 (Model V and VI) exhibit a temporal trend represented by a bell-shaped curve.  Di Gianfilippo et al. [51] used Model I to reproduce the transport of leaching substances from alkaline waste material assuming a constant source. They carried out a Monte Carlo analysis to evaluate a wide range of scenario conditions, varying the input parameters within selected intervals. The SAM value was estimated using a value of L 1 varying from 0.1 m to 1 m and a value of L f from 5 to 10 m, thus obtaining a value of C W /C W0 ranging between 9.9 × 10 −3 and 0.167. These concentrations are slightly lower than those calculated in this work because of the different geometrical conditions, as in the previous case.
The solutions of the group 2 consist of exponential laws, describing the depletion processes of the source. Model III assumes instantaneous mixing in the unsaturated zone, and it does not consider transport mechanisms. Therefore, C w /C w0 is spatially constant along the vertical axis, and it decreases over time. In Model IV the exponential law is multiplied by α leach . In this scenario, the value of α leach is approximately equal to 1, because of the selected type of soil and the contaminant. In consequence, the exponential part is dominant. The depletion coefficients µ 1 and µ 2 are equal to 3.469 × 10 −9 s −1 and 4.800 × 10 −9 s −1 , respectively. They are indicators of the depletion rate of the contaminant and differ only for the retardation coefficient. Therefore, in these conditions, the results of the two models are very similar. For group 2, the source depletion processes (bio-chemical degradation and leaching losses) are predominant, while the transport mechanisms are less relevant.
Models V and VI can be considered as the more advanced ones as they consider both the source depletion and a transient transport equation. This causes the distinctive bell-shape, yielding a dimensionless concentration peak lower than the other models I-IV. The solutions contain two parts: a time-dependent exponential law and the space-and time-dependent solution of the transport equation. For Model V the first part takes into account only the biochemical degradation, while for Model VI, it considers leaching and volatilization losses. The second part of the solution is different due to the initial and boundary conditions of the transport equation. For these reasons, the peaks of the solutions differ between the Models V and VI. In this scenario, for Model V the maximum of the normalized concentration is 4.338 × 10 −2 at 8.17 y, while in Model VI, the maximum of the normalized concentration corresponds to 4.920 × 10 −2 at 8.08 y. These models allow to describe the fate and transport mechanisms more thoroughly than the other models. In this case, the predicted concentrations are approximately an order of magnitude lower than those obtained by using the other models.
For these reasons, in some applications, the models of these groups are preferable. Di Sante et al. [72] used Model VI to assess the sanitary and environmental risk of a disused industrial plant contaminated with polycyclic aromatic hydrocarbons, heavy hydrocarbons and polychlorinated biphenyls. Although the Italian guidelines [73] consider adequate a steady-state model to reproduce leaching phenomena, the authors suggest the use of a transient model for contaminants with high experimental leaching value, because they are more suitable to reproduce real conditions. They evaluated concentrations of aromatic fraction of heavy hydrocarbons and benzo(a)pyrene after 20 years roughly equal to 0.40 mg/L e 1 × 10 −6 mg/L, respectively. Their results consider not only the leaching process but also the uniform dilution in groundwater, which is not considered by this paper. Furthermore, their results cannot be expressed as a dimensionless quantity. For these reasons, a direct comparison with the results of this work is not possible.
The different output concentrations of the models lead to different amounts of contaminant reaching the water table. In order to determine and compare the mass value reaching the groundwater for the different models, the mass flux (J) from unsaturated zone to groundwater is integrated over the time. J is governed only by water phase advection, and it is evaluated as suggested in [43]: where A uz-gw is the exchange area between the two zones. For all the models but Model II, A uz-gw is equal to A, and the dissolved concentration C w is constant over the area. In the paper describing Model II [43], A uz-gw is bigger than A, and the concentration is variable in space; hence, it is necessary to integrate C(x,y,t) over space and time.
For the considered contaminant, soil type and source thickness, the mass reaching the water table, considering a period of twenty years and forty years, has been calculated starting from an initial source concentration C w0 (100 µg/dm 3 ); see Table 6. As can be noted, the third group of models returns a lower contaminant mass than the others. After 20 years, the exponential models belonging to group 2 exhibit the highest values, while extending the period of time to 40 years Model I exhibits the maximum value of J.

Predicted Concentration of Different Contaminants over Time
The physico-chemical characteristics of the contaminant, e.g., the volatility and the degradation rate [52,74], strongly affect its mobility through the unsaturated zone as well as its fate in time. The selected models are solved using the contaminants properties listed in Table 5, fixing the thickness of the source L 1 = 2 m and a sand soil.
The transport time ranges between some orders of magnitude according to the contaminant mobility. Hence, in order to overcome this discrepancy and compare the results consistently, the dimensionless concentration has been plotted versus a dimensionless time τ, calculated as the ratio between the physical time and a reference time t ref . This criterion has not been used for Model I because the results are not depended by time and contaminant properties. In this case, t ref has not been considered, and the dimensionless concentration is plotted versus physical time. The value of the t ref cannot be evaluated the same way for all the models. For models belonging to the groups 2 and 3, the dimensionless concentration goes to zero after a certain period of time, depending on contaminants and medium properties. Therefore, the value of t ref has been set equal to the time when C w /C w0 reaches the residual value of 0.1%. For model II, the value of t ref has been calculated with a different criterion, as in this case C w /C w0 tends to a horizontal asymptote: It is set equal to the time when the variation of the concentration is negligible (<0.01 µg/dm 3 ). As one can perceive, in this case the value of t ref depends on the temporal discretization, which has been set equal to 1 month for Benzene, Vinyl chloride and Trichloroethylene and 1 year for Benzo(a)pyrene and Aldrin.
The values of t ref are listed in Table 7. It can be observed that Model III has the highest t ref for all the contaminants. More generally, the models belonging to group 2 exhibit a higher t ref than the other models for the contaminants characterized by a greater mobility (Benzene, Vinyl Chloride, and Trichloroethylene). For Benzo(a)pyrene and Aldrin, the reference time t ref of Model V is not reported because C w is negligible. For the Model II (Figure 3b), it is possible to observe a marked distinction in the trend concentration between contaminants characterized by a less mobility (Benzo(a)pyrene and Aldrin) and the others (Benzene, Vinyl chloride, Trichloroethylene). The first group of contaminants reaches a steady-state concentration higher than the others: the C w /C w0 ratio is equal to 0.443 for Benzo(a)pyrene and 0.358 for Aldrin. The second group reaches lower values, the C w /C w0 ratio is about equal to 0.090 for the three contaminants. These results are in accordance with the application of the model presented by Locatelli et al. [75], considering the different hydrogeological setting of the simulated site. Locatelli et al. [75] used the model II to simulate leaching of benzene and toluene in a simulated contamination scenario; at 5 meters, they obtain a C w /C w0 ratio corresponding to 0.01 for benzene and 0.04 for toluene. The lower value compared to the obtained value in this paper can be attributed to the different hydrogeological setting and the different assessment of some parameters, like the degradation constant.
The C w /C w0 ratios for Benzo(a)pyrene and Aldrin in the Model II are higher than in Model I (C w /C w0 = 0.286). This proves that in some circumstances, the Model I is not the most conservative approach, as it is typically considered (ASTM, 2000). This result can be related to the fact that in Model I the SAM is assumed simplistically constant for all the contaminants.
In Model III, as can be seen from Figure 3c, the results collapse on the same curve as consequence of the use of the dimensionless time τ. This occurs in Model IV to some extent. In fact, as previously explained, Model IV consists of two terms: α dep and α leach . The former is characterized by scale invariance while the latter is approximately equal to one for all the contaminants. The minimum value of α leach corresponding to 0.930 has been obtained for Benzo(a)pyrene. approach, as it is typically considered (ASTM, 2000). This result can be related to the fact that in Model I the SAM is assumed simplistically constant for all the contaminants.
In Model III, as can be seen from Figure 3c, the results collapse on the same curve as consequence of the use of the dimensionless time . This occurs in Model IV to some extent. In fact, as previously explained, Model IV consists of two terms: dep and leach. The former is characterized by scale invariance while the latter is approximately equal to one for all the contaminants. The minimum value of leach corresponding to 0.930 has been obtained for Benzo(a)pyrene. In Model V the concentrations of contaminants with low mobility (Benzo(a)pyrene and Aldrin) is closed to zero during the whole considered range of time (10 4 -10 5 years) because of the exponential part of the solution. Therefore, the results of this model (Figure 3e) are shown only for the three contaminants Benzene, Vinyl Chloride and Trichloroethylene. In Figure 3f. sketching the results of Model VI, it is possible to identify two groups of contaminants, characterized by low and high mobility, respectively. The concentration of contaminants belonging to the first group reach a maximum value earlier than the second group.

Comparison between Different Degradation Constants
The contaminant degradation into the environment, λ, simulated through a first order kinetics, takes into account both chemical mechanisms (hydrolysis, redox reductions and photodegradations) and biodegradation. Several mechanisms are involved, and each of them can be affected by many environmental variables, making the assessments of these coefficients highly uncertain [76,77].
The λ values can vary by more than one order of magnitude for each contaminant in a natural porous media [15,78]. For these reasons, a sensitivity analysis based on the variation of λ has been carried out. The comparison is developed with respect to the Trichloroethylene compound, in the same conditions of the previous sections: depth equal to z = 5 m, source thickness L 1 equal to 2 m and the unsaturated zone composed by sand. In Howard [66], the values of λ for Trichloroethylene in groundwater range from 4.852 × 10 −9 s −1 to 2.499 × 10 −8 s −1 . Three values have then been chosen for the comparison: the first λ 1 is the minimum value (4.852 × 10 −9 s −1 ), the second λ 2 is the mean of logarithms of the two extremes (1.101 × 10 −9 s −1 ), and the third λ 3 is the maximum value (2.499 × 10 −8 s −1 ). The results of the comparison are shown in Figure 4.
As already highlighted, Model I is independent of contaminant properties; therefore, the ratio C w /C w0 is constant, equal to 0.286 for all the values of λ. In Model II (Figure 4b), little differences related to the decay coefficients can be noted. The obtained C w /C w0 , equal to 0.0901, 0.0884 and 0.0848 for λ 1 , λ 2 and λ 3 , respectively, present similar values, reaching the steady state conditions in a small period of time.
(1.101 × 10 s ), and the third λ3 is the maximum value (2.499 × 10 s ). The results of the comparison are shown in Figure 4.
As already highlighted, Model I is independent of contaminant properties; therefore, the ratio Cw/Cw0 is constant, equal to 0.286 for all the values of λ. In Model II (Figure 4b), little differences related to the decay coefficients can be noted. The obtained Cw/Cw0, equal to 0.0901, 0.0884 and 0.0848 for λ1, λ2 and λ3, respectively, present similar values, reaching the steady state conditions in a small period of time.  In the Models III and IV (Figure 4c,d), the λ coefficients play the same role, affecting the depletion coefficients µ 1 and µ 2 , see Table 1.
The variation of λ affects Model V more than Model VI (Figure 4e,f). In particular, it influences the C w /C w0 peak (maximum concentration) as well as the time when it occurs. In both models, the higher is the value of λ, the earlier and lower is the peak. In Model V, C w /C w0 peaks correspond to 0.0438 at 8 y, 0.0094 at 7.83 y, 0.00033 at 7.42 y for the three values of λ considered, respectively, whereas in Model VI the C w /C w0 are equal to 0.0492 at 8.08 y, 0.0247 at 7.92 y and 0.0053 at 7.83 y.

Comparison between Different Soils
The physical properties of soil represent a further element that influences the fate and transport mechanisms of the contaminants. Simulations with three different soils (sand, loam, clay) have been carried out, in order to evaluate the corresponding differences in the models. Again, the comparison has been developed with respect to the Trichloroethylene compound in the same conditions depicted in the previous section. The three types of soils are characterized by the same organic carbon fraction, whereas they differ in grain size, water content, hydraulic conductivity and infiltration rate. The soil parameters values are listed in Table 4.
It is worth noting that the analyzed models consider the sorption process only through the soil-water partition coefficient, which is a synthetic parameter considering the characteristics of the soil, the contaminant and their interaction. This parameter should be assessed experimentally with reference to the specific soil-contaminant complex. In this analysis, it has been chosen to evaluate the partition coefficient referring only to the contaminant characteristics. For this reason, the obtained results for the different types of soils are affected only by the difference between hydrological characteristics of the soils, while the differences related to the sorption capacity are neglected. This assumption has resulted in a lower partition coefficient value of the clays than the experimental values in these soils and, consequently, a higher value of the simulated contaminant concentration. Clay soils have, in fact, a higher sorption behavior than the other soils because of the high surface area and negative charged sites of the soils particles.
As previously seen, Model I is independent by the soil parameters (Figure 5a). Model II (Figure 5b) shows a similar behavior for sand and loam, reaching the steady state C w /C w0 ratios corresponding to 0.0902 and 0.0836, at 1.67 y and 2.42 y respectively. A lower mobility could be observed for clay. In this case the steady concentration, C w /C w0 equal to 0.0538, is reached at t = 7.67 y. Hence, the leaching velocity has a greater impact compared to the sorption capacity.
partition coefficient referring only to the contaminant characteristics. For this reason, the obtained results for the different types of soils are affected only by the difference between hydrological characteristics of the soils, while the differences related to the sorption capacity are neglected. This assumption has resulted in a lower partition coefficient value of the clays than the experimental values in these soils and, consequently, a higher value of the simulated contaminant concentration. Clay soils have, in fact, a higher sorption behavior than the other soils because of the high surface area and negative charged sites of the soils particles.
As previously seen, Model I is independent by the soil parameters (Figure 5a). Model II ( Figure  5b) shows a similar behavior for sand and loam, reaching the steady state Cw/Cw0 ratios corresponding to 0.0902 and 0.0836, at 1.67 y and 2.42 y respectively. A lower mobility could be observed for clay. In this case the steady concentration, Cw/Cw0 equal to 0.0538, is reached at t = 7.67 y. Hence, the leaching velocity has a greater impact compared to the sorption capacity. Models V and VI are most influenced by soil parameters, which affect both the time at which the maximum concentration occurs at z = 5 m and the related values. The transport mechanism lasts very long when Model V is applied for clay. Hence the concentration is so extremely low and therefore undetectable. The C w /C w0 peaks are equal to 4.380 × 10 −2 at 8 y and 8.651 × 10 −3 at 18.67 y in sand and loam, respectively. In Model VI the peaks are of the same order of magnitude for sand and loam, whereas it is negligible for clay. In this case, the C w /C w0 ratios are equal to 4.920 × 10 −2 at 8.08 y and 1.509 × 10 −2 at 19.75 y for sand and loam, respectively.

Comparison between Different Source Thickness
As can be seen in Table 1, the thickness of the source L 1 determines the soil attenuation of the contaminant in the Model I, and the time needed to deplete the source in the models composed of an exponential part (Models III to VI). Model II is not influenced by the change of source thickness because it does not take into account the source depletion and soil attenuation. The simulations have been carried out taking into account the scheme depicted in Figure 1, with three different source depths: L 1 = 1 m, 2 m and 4 m, using sand as type of soil and Trichloroethylene as contaminant. Again, the dimensionless concentration C w /C w0 has been calculated at the depth z = 5 m. The concentrations calculated with Model I (Figure 6a)

Conclusions
The use of fate and transport models may give a considerable support to tackle with environmental problems involving chemical leaching in the vadose zone. Unlike numerical based approaches, analytical models are used especially in the risk assessment of contaminants due to their straightforward application and simplicity. Six analytical models, reproducing contaminant fate and transport from a source located in the unsaturated zone to the groundwater were described and compared with the aim of giving some insights about their use.
On the basis of their characteristics three groups were identified. The first one is suited for contamination scenarios where the source may be assumed constant, e.g. in the case of continuous chemicals leakage. The second group is particularly appropriate when the depletion processes are predominant and the travel distance of the chemical is low, e.g. a localized spill of biodegradable organic contaminant near the water table. The third group describes more comprehensively all the processes occurring during the transport.
Simulations being carried out, showed that models belonging to the third group yielded lowest values of contaminant concentration with respect to the other models, giving the less conservative results. When a more conservative approach is needed, Models I to IV should be preferred instead.
Future research could be addressed to either enhance the existing analytical models or to propose new advanced ones, while preserving their suitability for engineering applications, in order to obtain a more realistic representation of contaminant fate and transport in the unsaturated zone. In the exponential Models III and IV (Figure 6c,d), the source thickness has a different effect. The depletion coefficients µ 1 and µ 2 ( Table 1)  In the Model V, an increased thickness of the source determines delayed and higher peaks. The C w /C w0 ratios are 2.25 × 10 −2 (t = 7.92 y), 4.38 × 10 −2 (t = 8.00 y) and 8.24 × 10 −2 (t = 8.25 y) for L 1 equal to 1 m, 2 m, 4 m respectively.

Conclusions
The use of fate and transport models may give a considerable support to tackle with environmental problems involving chemical leaching in the vadose zone. Unlike numerical based approaches, analytical models are used especially in the risk assessment of contaminants due to their straightforward application and simplicity. Six analytical models, reproducing contaminant fate and transport from a source located in the unsaturated zone to the groundwater were described and compared with the aim of giving some insights about their use.
On the basis of their characteristics three groups were identified. The first one is suited for contamination scenarios where the source may be assumed constant, e.g. in the case of continuous chemicals leakage. The second group is particularly appropriate when the depletion processes are predominant and the travel distance of the chemical is low, e.g. a localized spill of biodegradable organic contaminant near the water table. The third group describes more comprehensively all the processes occurring during the transport.
Simulations being carried out, showed that models belonging to the third group yielded lowest values of contaminant concentration with respect to the other models, giving the less conservative results. When a more conservative approach is needed, Models I to IV should be preferred instead.
Future research could be addressed to either enhance the existing analytical models or to propose new advanced ones, while preserving their suitability for engineering applications, in order to obtain a more realistic representation of contaminant fate and transport in the unsaturated zone.