Estimation of Conservative Contaminant Travel Time through Vadose Zone Based on Transient and Steady Flow Approaches

Estimation of contaminant travel time through the vadose zone is needed for assessing groundwater vulnerability to pollution, planning monitoring and remediation activities or predicting the effect of land use change or climate change on groundwater quality. The travel time can be obtained from numerical simulations of transient flow and transport in the unsaturated soil profile, which typically require a large amount of data and considerable computational effort. Alternatively, one can use simpler analytical methods based on the assumptions of steady water flow and purely advective transport. In this study, we compared travel times obtained with transient and steady-state approaches for several scenarios. Transient simulations were carried out using the HYDRUS-1D computer program for two types of homogeneous soil profiles (sand and clay loam), two types of land cover (bare soil and grass) and two values of dispersion constant. It was shown that the presence of root zone and the dispersion constant significantly affect the results. We also computed the travel times using six simplified methods proposed in the literature. None of these methods was in good agreement with transient simulations for all scenarios and the discrepancies were particularly large for the case of clay loam with grass cover.


Introduction
Reliable determination of contaminant travel time in an unsaturated zone (also known as an unsaturated time lag) is required to assess groundwater vulnerability to pollution, to plan monitoring and remediation activities and to predict the effect of changes in climate or land use practices on groundwater quality [1][2][3].While the susceptibility of aquifers to pollution can be quantified using different approaches, like the index methods (e.g., DRASTIC [4,5]) or the grey water footprint method (GWF) [6,7], methods which explicitly take into account the time scale of contaminant transport are often preferable [3].For the purposes of travel time estimation, it is often assumed that the dissolved pollutant behaves as a conservative substance; for example, it does not undergo any chemical or biological reactions and is not sorbed on the soil skeleton.The resulting travel times are considered as safer (shorter) estimates, compared to the case of non-conservative contaminants, which are subject to retardation or attenuation.Over the years, a number of approaches have been developed to estimate the time lag in the vadose zone, ranging from simple analytical formulas to comprehensive numerical models.
Analytical equations were developed for the case of steady downward water flow, purely advective transport and constant water content in the soil profile, or at least in each soil layer [1,8,9].For example, in Polish hydrogeological practice, one often applies the formulas proposed by Bindemann (as described in [10]), Witczak and Żurek [11] and Macioszczyk [12].They require only basic data, such as the thickness of soil layer, average recharge flux and average water content.It should be noted, however, that accurate estimation of the recharge rate still represents one of the major challenges in hydrogeology [13].Additionally, there is significant uncertainty in estimation of the average water content, which generally varies with both time and depth.Even for steady flow in homogeneous soil, the water content between the groundwater table and the soil surface is not uniform.This leads to large discrepancies between some of the simple formulas, as reported for example in [12,14].Sousa et al. [15] recommended calculating the travel time by numerical integration using a variable water content profile.The profile can be obtained either from the numerical solution of the steady unsaturated flow equation for a specific value of downward flux (recharge), or from the assumption of hydrostatic equilibrium.This method requires knowledge of soil hydraulic functions, describing relationships between the water pressure, water content (or saturation) and hydraulic conductivity.
An alternative approach is to carry out transient simulations of flow and transport processes between the soil surface and the water table, which are based on the Richards' equation combined with an equation describing the advective-dispersive solute transport.In this way, one obtains estimations of both water and solute fluxes at the water table on condition that detailed weather, soil and vegetation data are available as the model input.This approach is gaining increasing popularity; however, the results were shown to depend on many factors, including the time resolution of the weather data [16][17][18].Recently, Vero et al. [19] presented a comprehensive framework for estimating travel time based on simulations with the HYDRUS-1D computer code [20].
In spite of the growing application of sophisticated numerical models, simple formulas remain a valuable tool, especially in preliminary or large-scale studies, where the detailed data are not available.There is a need for more detailed comparison of the results of simple methods with numerical simulations based on the Richards' equation.In the present study, we focus on such a comparison for two soil profiles (sand and clay loam), which can be considered typical for post-glacial areas in Northern Poland.Numerical simulations are performed using the HYDRUS-1D code, which is a well-established tool in modeling recharge rates and travel times [19,[21][22][23].We evaluate the role of the hydrodynamic dispersion, which is represented in the transient transport equation, but not accounted for in the simple methods.We also investigate the role of the root zone.Travel times obtained from the transient simulations are compared to six simpler methods available in the literature.To the best of our knowledge, a systematic comparison of steady-state and transient solutions has not been presented in the literature so far.Our study shows that the agreement between simple methods and numerical simulations of time-variable flow heavily depends on the method and the considered scenario.None of the methods was capable of closely reproducing the transient simulation results in all cases.The discrepancies were larger for the clay loam than for the sand.

Vadose Zone Profiles
We considered two homogeneous soil profiles, composed of sand and clay loam, respectively.The depth to the groundwater table is 6 m in each case.We followed an approach commonly applied in the vadose zone hydrology and assumed that the relationships between the water pressure h (negative in the unsaturated zone), the volumetric water content θ and the hydraulic conductivity (permeability) k of each soil were described by the van Genuchten model [24]: where S e is the effective (normalized) water saturation, θ r is the residual water content, θ s is the water content at fully saturated (or field saturated) conditions, α is a parameter related to the average pore size, n g and m g are parameters related to the pore size distribution (m g = 1 − 1/n g ), k s is the hydraulic conductivity at saturation, and k r is the relative hydraulic conductivity.For simplicity, we neglected hysteresis, soil deformation and air trapping.While Equations (1a) and (1b) are well established in the literature, their limitations should be noted, in particular with regard to θ r , which is often defined as the water content which cannot be removed from soil in field conditions.According to the above equations, h tends to negative infinity and k tends to 0 as the water content approaches θ r .On the other hand, it has been shown that the limit value of the negative pressure potential is equivalent to about 10 7 cm for oven dry state (θ = 0) [25].Moreover, water flow is possible in thin films even in very dry soils [26].Thus, θ r does not have a clear physical interpretation and is treated rather like a fitting parameter of the model, while the model is supposed to be used in the range of θ much larger than θ r .Parameters of the van Genuchten model for each soil were taken as average values, reported by Carsel and Parrish [27] and implemented in the HYDRUS-1D data base.They are listed in Table 1.In order to keep the present study concise, we focused on homogeneous soils and did not consider layered soil profiles, which commonly occur in field condition.All methods described in this paper can be also applied to layered soils, but it can be expected that the influence of heterogeneous structure depends on the position and relative thickness of individual layers, which makes it a topic for a separate study.Similarly, we did not take into account the possible effects of preferential flow, caused by the presence of macropores or other factors.The preferential flow has a large impact on the travel time of contaminants [28]; however, it was not represented in the analytical equations for travel time considered in the present study.

Numerical Modeling
The reference estimation of the contaminant travel time was based on numerical solutions of the differential equations describing the unsteady water flow and contaminant transport in the vadose zone, which allows obtaining estimates of both the groundwater recharge flux and contaminant travel time.For this purpose, we used the HYDRUS-1D computer program [20].The water flow was described by the Richards equation, which couples Darcy's law with the mass conservation equation as follows: where t is the time, z is the vertical coordinate (oriented upwards), and S(h) is the sink function representing water uptake by plant roots.In order to solve Equation ( 2), the HYDRUS-1D used the finite element discretization in space to solve the Richards' equationand fully implicit discretization in time.
In our case, simulations were carried out for 5 years (sand) and 60 years (clay loam).Each profile was uniformly discretized with a node spacing of 0.6 cm (1001 nodes in total).Distributions of initial pressure head along the profiles were calculated in preliminary simulations for a period of 5 years.The upper boundary condition reflected variable atmospheric conditions (evaporation and precipitation fluxes, with free surface runoff), while the lower boundary condition was assigned as a constant pressure head (water table level, h = 0).The weather data, such as the daily minimum and maximum temperatures, amount of precipitation, solar radiation, wind speed and air humidity observations from 2011 to 2015, were obtained from a meteorological station of Gda ńsk University of Technology in Sopot.For simulation of flow in the clay loam profile, the 5-year data series was repeated 12 times to obtain a 60-year period.The annual precipitation varied between 515 and 577 mm/year, with an average of 550 mm/year.Potential evapotranspiration was estimated with the Penman-Monteith FAO equation [29].To simulate vegetative cover (grass), the actual root water uptake was estimated with the Feddes' model [30] based on the default parameters for grass implemented in the HYDRUS-1D.The thickness of the grass root zone was set to 0.5 m [31].
Migration of the conservative solute was simulated with the advective-dispersive transport equation implemented in HYDRUS-1D (here shown in a simplified form, applied in the present study): where c is the solute concentration, D is the hydrodynamic dispersion coefficient and q is the Darcy velocity.The above equation includes dispersion, which is not considered in the analytical formulas for travel time described below.For the 1D transport, D = α L |q| +D m *, where α L is the longitudinal dispersivity (or dispersion constant) and D m * is the effective molecular diffusion.
In our simulations, D m * was set to the chloride diffusion coefficient in bulk water (D m = 1.2 cm 2 /day) multiplied by a tortuosity factor depending on the water content, as implemented in the HYDRUS-1D.
We considered two cases with different dispersivity values, α L = 0.6 m and 0.06 m, respectively.The first case represents the typical choice of dispersion constant, equal to 10% of the length of the transport domain [16].In the second case, the dispersion was reduced, in order to improve consistency with the analytical methods, which were developed for purely advective cases.Note that it was not possible to entirely remove the dispersive term, due to the instability and errors of the numerical algorithm.The algorithm used the Galerkin finite element method for spatial discretization of the transport equation, while the Crank-Nicholson scheme was chosen for time discretization.Initially, the concentration was assumed to be 0 in the whole profile.The boundary condition at the soil surface was assigned as a third-type condition (specific contaminant flux, depending on the infiltration), with a constant concentration of contaminant in the infiltrating water c = 1 mg/cm 3 .We chose this kind of boundary condition, because it is more physically based than specifying a constant concentration at the boundary [20] and also consistent with the assumption of piston-flow model underlying the simple analytical methods.At the bottom of the profile, the zero concentration gradient condition was assumed.
The resulting time of vertical transport was reported as the time of arrival of the concentration equal to 0.01 mg/cm 3 at the bottom boundary of the model (watertable).We also reported the arrival time for the concentration of 0.99 mg/cm 3 .The difference between these two times depends on the assumed dispersion coefficient.

Analytical Equations for Vertical Travel Time
The time lag of dissolved contaminants in the vadose zone is commonly estimated using analytical equations based on some simplifications of the flow and transport processes occurring in real-life conditions.It is assumed that the conservative solute is transported only by advection, that is, it moves at the same velocity as the infiltrating water.Moreover, the steady water flow is considered, with the Darcy velocity q uniform along the soil profile and equal to the recharge rate R.However, the actual advection velocity varies, because of the variations in water content.The steady state profile of water pressure and water content can be obtained by solving the following equation, which is a simplification of Equation ( 2): At any elevation z, the advective velocity v a and the increment of travel time dt a over the depth increment dz are given by: Integrating Equation ( 5) over the whole depth of the profile provides the total travel time t a .This approach was suggested by Sousa et al. [15].Another possibility, also discussed in [15], is to use the vertical water content profile corresponding to the hydrostatic equilibrium, which can be easily obtained from the water retention curve, without the need to solve Equation (4).While the hydrostatic equilibrium represents no-flow conditions, it can be considered a reasonable approximation of the water content variability above the groundwater table in many real-life situations.In the following section, we present results obtained for both types of water content profiles.The steady-state solution was obtained from HYDRUS-1D, by imposing a constant water flux equal to recharge as the top boundary condition.For the sake of consistency, we applied the average recharge fluxes calculated from the transient flow simulations described in Section 2.2.
In the simplest methods for travel time estimation, the water content was assumed to be uniform within each homogeneous soil layer and equal to θ av .In such a case, the advective velocity v a and the corresponding time of travel t a through a soil layer of thickness L can be calculated as [11,32,33]: For any given recharge flux R, the accuracy of estimation crucially depends on the choice of θ.Some estimations can be obtained from the literature or field measurements.In Poland, a rather detailed guidance is available in the works of Witczak and Żurek [11] and Duda et al. [34].One should note that these authors suggest the use of the total volumetric water content in Equation ( 6), which means that the whole volume of water in pores is considered to participate in the advective transport of contaminant.This is consistent with the transport Equation (3) implemented in HYDRUS-1D and also described in other works [32,33].However, some authors [8] proposed to use the mobile (or effective) water content in Equation ( 6), which is defined as the part of water participating in flow and advective transport.The latter approach seems to be more consistent with the concept of effective porosity n e used in contaminant transport models for the saturated zone [35].The effective porosity, applied to calculate the advective velocity in the saturated conditions, is considered smaller than the total open porosity, and the difference is particularly significant in fine-textured media like loams or clays.
In the present study, Equation ( 6) was evaluated for the volumetric water content θ av equal to 0.07 and 0.1 in sand and 0.24 to 0.32 in clay loam, representing the expected range of variability [1,34].
If one assumes steady flow, the value of average water content in the soil profile θ av is related to the amount of recharge.For the same soil, larger recharge fluxes correspond to larger θ av values.A simple relationship can be obtained assuming that the water pressure gradients are negligible and the flow is driven mostly by the gravity potential.In such a case, Equation (4) was reduced to: Such a simplification seems more suitable for coarser sands, rather than for fine-textured soils.Nevertheless, it allows us to calculate a constant value of water content corresponding to the specific recharge flux R for any type of hydraulic conductivity functions.For the van Genuchten function given by Equation (1b), one has to solve a nonlinear equation; however, a closed-form solution can be easily obtained if one uses a power law (Brooks-Corey type function): where b is the soil-dependent parameter.Using the above expression in combination with Equation ( 7), Charbeneau and Daniel [9] (cited in [33]) obtained: In the present study, b was set equal to 4.19 for sand and 9.45 for clay loam, based on the average parameters of the Brooks-Corey functions for these soil textural classes reported in [33].
If we assume b = 3, θ − θ r as equivalent to the mobile water content and θ s − θ r as equivalent to the effective porosity n e (bearing in mind the ambiguity of θ r mentioned above), the result is the Bindemann formula, as presented by Szestakow and Witczak [10]: In the present study, the Bindemann formula was evaluated for sand with n e = 0.2 [34] and n e = θ s − θ r = 0.385.For clay loam, we chose n e = 0.1 [34] and n e = θ s − θ r = 0.315.
The Bindemann formula was further modified by Macioszczyk [12], who suggested the replacement of n e in Equation ( 10) by the total volumetric water content: Macioszczyk [12] argued that Equation (11) predicts travel times which are more realistic compared to either Equation (6) or Equation (10).However, there seems to be no physical basis for such a modification.We evaluated Equation (11) for the same range of water content values as Equation ( 6) above.

Results
Examples of vertical profiles of the volumetric water content obtained from HYDRUS-1D simulations of transient flow are shown in Figures 1 and 2. It can be seen that the volumetric water content is fairly uniform in sand, in particular at depths larger than 2 m, where the influence of evaporation and transpiration disappears.The range of computed values is consistent with that reported in [11], i.e., 0.07 for gravelly sands to 0.1 for fine sands while it is somewhat smaller than the values given by [34] (from 0.08 for gravelly sands to 0.17 for fine sands).In clay loam, the computed variations of the water content are significantly larger and reach the depth of 4 m in the case without vegetation.The obtained values are generally consistent with those reported for loams (Polish: gliny) by either Witczak and Żurek [11] (from 0.26 to 0.32) or Duda et al. [34] (from 0.24 to 0.32).The average annual values of recharge (Table 2) computed by HYDRUS-1D show very significant influence of the vegetative cover, which reduces the recharge by a factor of about 2 for sand and about 4 for clay loam.The results for sand correspond very well with the graph of Dyck and Chardebellas [36], as reproduced in [37], based on lysimeter experiments in Central Europe.The values computed for clay loam are between the lines for sandy loam and loam in the same graph.The obtained recharge/precipitation ratios for sand are generally higher than those given in the Polish literature.For example, Duda et al. [34] suggest, in the case of typical vegetation, the ratio from 0.14 for fine sand to 0.25 for gravelly sand, which should be multiplied by 1.2 for bare soil conditions.For loams, the same authors give ratios ranging from 0.08 to 0.12, which are in between the simulation results for vegetated and bare soil cases.The travel times of conservative solute obtained from numerical simulations of transient flow and transport are shown in Table 3.The early appearance of contaminant at the water table (c = 0.01 mg/dm 3 ) is strongly influenced by the dispersion coefficient, especially in sand where the travel time is 5 to 8 times shorter for large dispersion case compared to the small dispersion case.For clay loam, the differences are smaller, but still very significant, with the difference in travel time by a factor of 2. On the other hand, the differences in arrival time of the high concentration (c = 0.99 mg/cm 3 ) are much less significant and do not exceed 35%.It seems that a careful examination of the role of dispersion is necessary in studies related to travel time estimations, especially if the contaminant is considered harmful even in small concentrations.Travel time estimates based on the assumption of advective transport and steady-state or hydrostatic θ profiles are given in Table 4.As expected, the hydrostatic profile assumption results in shorter travel times; however, the differences are not very big.The largest relative difference between the steady flow and hydrostatic cases occurs for bare sand, the smallest one for vegetated clay loam.For sand, the hydrostatic case seems to be in a slightly better agreement with transient simulations than the steady flow case, but the latter one is also reasonably close to the transient results.However, for clay loam, both types of profiles lead to much longer travel times than the transient simulations.This could be explained by Figure 2, which shows that the water content profiles in clay loam are more variable in time than the profiles in sand.Table 5 presents travel times calculated with simple analytical formulas.For Equations ( 6), ( 10) and ( 11), we reported a range of values, corresponding to the assumed range of variability of θ av or n e , as explained in Section 2.3.It can be seen that the results vary greatly between the formulas.For sand, Equation ( 6) seems to be in relatively good agreement with the results from transient simulations (small dispersion case).For clay loam, the estimated travel time was significantly longer than the one obtained from HYDRUS-1D, especially in the case of vegetation.It was possible to obtain good match between the two approaches if θ av was chosen in the range from 0.11 to 0.16 for the bare soil case and in the range from 0.15 to 0.22 for the vegetation case.However, these values are smaller than the typical ones reported in the literature and they are not consistent with the values calculated in numerical simulations (Figure 2).The method developed by Charbeneau and Daniel [9] gave travel times within the range predicted by Equation ( 6).In contrast, both Bindemann's and Macioszczyk's formulas in all scenarios gave travel times much shorter than the ones computed from HYDRUS-1D and other methods.This discrepancy has been shown previously by several authors [8,10].For sand, the equation of Macioszczyk leads to the shortest travel times, because the volumetric water content is smaller than the effective porosity used in the Bindemann formula.For clay loam, if one uses small values of the effective porosity, as commonly reported in the literature, the Bindemann formula predicts shorter time lag than Macioszczyk formula.

Summary and Conclusions
The numerical simulations of time-variable flow and transport showed that the groundwater recharge flux strongly depends on the presence of root zone, for either sand or clay loam, although the difference is larger in the case of clay loam.Consequently, the travel time of pollutant also depends on the presence of root zone.The assumed dispersion constant has significant influence on the arrival time of contaminant at the water table.While this conclusion may seem obvious, it seems to be important in view of the widespread calculation of travel time based on the assumption of purely advective flow.
The methods using the steady flow approximation showed mixed performance, even though it was assumed that the exact value of average groundwater recharge is known for each soil profile.In the case of sand profiles, the simplest method given by Equation ( 6) used in combination with literature-based range of water content values gave results comparable to the transient flow analysis.Taking into account vertical variability of the water content distribution [15] or the relationship between the water content and the recharge rate [9] did not offer any significant improvement over Equation (6).On the other hand, the methods of Bindemann and Macioszczyk led to severely underestimated travel times.This point may have important consequences, as both methods are widespread among hydrogeology practitioners in Poland.
In the case of clay loam, most of the steady-state methods predicted travel times significantly longer than those resulting from the simulations of unsteady flow.It seems that the transient factor may be more important in this case.Again, the Bindemann and Macioszczyk methods gave travel times much shorter than the HYDRUS-1D computations, by a factor of at least 2 or 3.
Care should be taken if simple analytical formulas are to be used to estimate unsaturated zone travel time.In view of the growing computer capacities and availability of simulation software and parameter data, it seems advisable that numerical simulations are used for at least partial comparisons with the analytical results.One advantage of numerical models of flow and transport in unsaturated zones is that both water and contaminant fluxes reaching the groundwater table can be estimated.Additionally, one dimensional model for unsaturated zone processes can be integrated with 3D models for flow and transport in the saturated zone, such as MODFLOW or MT3D, which is currently undertaken by the authors [18,38,39].
Numerical calculations presented in this paper were based on the so-called single porosity model, which assumes that the whole volume of pores occupied by water contributes in the same proportion to the advective flow.In soil science, the phenomenon of preferential or non-equilibrium flow has been long recognized [28,40], with rapid water flow in larger pores bypassing smaller pores.This leads to faster movement of solutes than that predicted by the standard transport Equation (3).However, the largest pores become active only in conditions close to full saturation, while drained state flow occurs mostly in smaller pores.The concept of non-equilibrium flow potentially allows reconciling the small values of effective porosity characterizing transport in saturated fine-textured soils with relatively large values of the water content occurring in such soils during unsaturated flow.Further research on this topic is warranted.Additionally, the relationship between steady-state and transient methods for soil profiles composed of distinct layers remains to be investigated.

Figure 1 .
Figure 1.Volumetric water content profiles obtained from HYDRUS-1D simulations for bare sand (a) and sand with grass cover (b).Initial profiles shown in black, dark blue, green, light blue and red correspond to the end of the first, second, third and fourth year of the simulation, respectively.

Figure 2 .
Figure 2. Volumetric water content profiles obtained from HYDRUS-1D simulations for bare clay loam (a) and clay loam with grass cover (b).Initial profiles shown in black, dark blue, green, light blue and red correspond to the end of the first, second, third and fourth year of the simulation, respectively.

Table 1 .
Parameters of the van Genuchten hydraulic model for each type of soil.

Table 2 .
Annual recharge and recharge/precipitation ratios obtained from numerical simulations.

Table 3 .
Estimates of travel time (days) from numerical simulations of transient flow.

Table 4 .
Estimates of travel time (days) for steady flow and hydrostatic θ profiles.

Table 5 .
Estimates of travel time (days) from analytical equations 1 .