Modeling Water Infiltration and Solute Transfer in a Heterogeneous Vadose Zone as a Function of Entering Flow Rates

Due to its rapid movement, preferential flow (PF) in the vadose zone allows much faster contaminant transport, which may have a significant impact on ground-water quality. PF can occur in heterogeneous vadose zones and it strongly depends on hydric and hydraulic conditions like entering flow rates at surface. This study deals with the modeling of the establishment of PF, and related solute transfer during the infiltration phase in a strongly heterogeneous glaciofluvial deposit. This deposit is made of four contrasting lithofacies (sand, gravel, bimodal gravel and matrix-free gravel) and lies underneath an urban infiltration basin (Lyon, France). Previous studies have been carried out on this site and linked the regionalization of soil pollution with the lithological heterogeneity. But none of them clearly demonstrated how heterogeneity could impact flow and solute transfer and may explain such a regionalization. In this study, we model flow and solute transfer at the trench scale for both uniform and heterogeneous profiles in order to characterize the effect of lithological heterogeneity. In addition, such a modeling was performed for two different entering flow rates to depict the influence of condition at surface on PF. A key result is that heterogeneity clearly impacts unsaturated flow and solute transfer. Numerical modeling permitted pointing out the existence of PF paths associated with the sedimentary heterogeneity of the glaciofluvial deposit. For lower surface fluxes, the sand lens and matrix-free gravel were the sources of capillary barrier effects, leading to a funneled flow and a groundwater recharge characterized by earlier and more dispersed wetting fronts. Such a flow pattern enhances solutes transfer and reduces solute retention by soil. Thus, the effect of heterogeneity on solute transfer is significant, especially for the most reactive solutes. E. Ben Slimene et al.


Introduction
Large cities are mostly situated in areas close to water resources in order to meet the water needs of their populations.Alluvial soils harbor large aquifers that are used to supply water.The Rhone-Alpes region is a good illustration of it.However, the increase of soil sealing has led to the development of the best management practices such as infiltration basins aimed at infiltrating storm water in order to reduce the amount of collected and treated water in usual systems.Yet, these infiltration basins are mainly settled over highly permeable geologic formations so as to ensure water infiltration and thus better pollutant retention.Most of these formations are strongly heterogeneous, since they are made of different materials with contrasting sedimentological properties (e.g.particle size distribution) and transfer properties [1]- [3].
Soil heterogeneity is responsible for the difficulty in predicting mass movement.For example, it often results in faster mass movement than would be expected from the soil matrix properties [4] [5].This movement is associated with PF that denotes water and solute transport through a small portion of a soil volume receiving input over its entire inlet boundary.Transport along PF paths can cause water and chemicals to move more rapidly and reach greater depths than predicted by models, i.e. the Richards equation.PF can occur in macropores [6], as finger flow [7], as funnel flow [8] or in water repellent soils [9].There have been several studies of PF in soils during the last 30 years [10] [11].As a result, PF is now routinely included in models that predict water, solute and particle transport in soils [12].
Most modelling approaches for flow and solute transfer rely on the use of darcean approach, i.e. assuming that flow is ruled by the Darcy Buckingham Law [13].In such a case, unsaturated water flow in the vadose zone is modeled considering the Richards equation, written as follows, considering no source/sink terms: where h stands for the pressure head, θ is the volumetric water content, x i , x j are the spatial coordinates, t is time, and K ij are components of the unsaturated hydraulic conductivity tensor.We assumed the soil in each layer to be isotropic, with both entries K xx and K zz , equal to the unsaturated hydraulic conductivity function, K(h).Soil water retention and hydraulic curves are characterized by the van Genuchten-Mualem formulation [14]: where θ r and θ s denote the residual and saturated water contents, respectively; K s and K r are the saturated and relative hydraulic conductivities, respectively; h g = 1/α is the scale parameter of water pressure heads, n is a pore-size distribution index, l is a pore-connectivity parameter, and m = 1 − 1/n.Solute transfer is often modeled considering the Convection-Dispersion Equation (CDE).This equation assumes that solutes are transported through convection and dispersion in the liquid and may be retained in the soil through sorption onto solid particles [13]: where C stands for the liquid concentrations and S for the sorbed concentrations, D ij is the dispersion tensor, described as well by the longitudinal and transversal dispersion coefficients, D L and D T , respectively, which account for both molecular diffusion and mechanical dispersion: where λ L and λ T are the medium longitudinal and transversal dispersivities, τ is the medium tortuosity, and D 0 is the solute molecular diffusion coefficient.Concerning pollutant retention, sorption is assumed to be instantaneous, kinetically limited mechanisms being the focus of further research.The total sorbed concentration, S, is described as a function of the liquid concentration, C, with the following general expression [15]: where K d , β and η are empirical coefficients that can be estimated from experimental sorption isotherms [16].
For some pollutants, and especially for organic pollutants, sorption isotherms may be linear and the estimation of K d may be the sole parameter required for modeling solute transfer.These flow and solute equations can be adapted to account for preferential flow and solute transfer [17] [18].Preferential flow (PF) can be seen as uneven and rapid movement of water and solutes through restricted zones of the porous media (e.g.cracks, wormholes, etc.), with slow or even negligible movement of water and solutes elsewhere [17].In this case, differential equations can be changed to consider the concomitancy of fast flow along preferential pathways and pockets of stagnant water.Solutes undergo convection and dispersion along the main stream with solute diffusion between the "mobile" water and stagnant water pockets [19].These equations are applied at porous medium scale (i.e. at the centimeter-meter scale); the porous medium is described as a set of representative elementary volumes (REV).At larger scales, the ground may exhibit an additional spatial heterogeneity of transfer properties, specifically for deposits that are constituted of several lithofacies.When processes need to be described at larger scales (decameter or larger), the lithological heterogeneity of the ground or the deposit must be accounted for, and the ground has to be considered as an assemblage of several materials, each requiring a proper description of its transfer properties [3].The combination of these two types of heterogeneity, including the intrinsic heterogeneity related to each material (potential macropores, cracks, etc.) and the lithological heterogeneity resulting from the architecture of the deposit, is a real challenge.Studying the impact of the lithological heterogeneity is a prerequisite to the understanding of flow and mass transfer in heterogeneous vadose zones.
This paper addresses flow and solute modeling during the infiltration phase in the vadose zone of the heterogeneous deposit of Lyon area.In particular, we want to pinpoint numerically the influence of lithological heterogeneity on the establishment of PF and mass transfer.We also want to identify the hydric and hydraulic conditions (mostly entering flow rates) that are favorable to PF and mass transfer.For this purpose, the proposed numerical model is based on several studies that propose a complete sedimentological description of the deposit architecture and lithofacies along with their hydraulic properties [1]- [3].Considering these data, water infiltration and solute transfer were modeled for two imposed flux values.Preliminary studies focused on numerical options and mesh optimization with regards to both accuracy and computation time.Data derived from these tests were used in the proposed study.Then, a sensitivity analysis addressed the role of the deposit architecture and boundary conditions (entering flow rate) in the establishment of PF.For flow, the time required reaching steady state over the whole soil profile, wetting fronts, water fluxes during transient state and water fluxes at steady state were accurately characterized and compared to those of a uniform section.Solute transfer was investigated at steady state for a tracer and two reactive solutes with two different distribution coefficients (parameter K d , in Equation ( 7)).Then, the influence of lithological heterogeneity on entering flow rate was discussed.

Study Site
East Lyon is located in a filled area of the Rhone Valley settled on an ancient formation dating from the last Glacial Maximum.The study site (GPS coordinates: 45.7360˚N/4.9570˚E) is an infiltration basin situated beneath this deposit and called Django Reinhardt-DjR (Figure 1) and has already been the subject of several studies [3].This basin, covering 13,775 m 2 , receives runoff storm waters from an urban and industrialized watershed of 185 ha in surface area [1] [2] [20]- [22].The depth of the unsaturated zone under the infiltration basin is 13 m.This site is instrumented by the field observatory for urban water management (http://www.graie.org/othu/).To characterize the lithological properties of the deposit, a trench of a length of 13.5 m and a height of 2.5 m has been mechanically excavated and is depicted in Figure 2(a) [2].This profile was the subject of a detailed lithographic description [1] [2].The studied trench was composed of four structural units made of four main materials.Referring to particle size distribution (cumulative fraction and frequency), we can distinguish in Figure 2(c) the following materials: (1) a poorly-sorted to moderately sorted medium sand with a grain size of 0.40 ± 0.21 mm, (2) a bimodal mixture of fine to medium sand and medium to coarse gravel, (3) a mixture of sandy matrix and gravel with a gravel fraction of 85%, and (4) a poorly sorted to moderately well sorted matrix-free gravel.These materials are referred to as sand, bimodal gravel, gravel, and matrix-free gravel, respectively.The gravel forms the upper layer, around 50 -70 cm thick, and results from the blending of all original lithofacies by mechanic engines while monitoring the basin and removing the upper layer.The bimodal gravel makes up most of the deposit (Figure 2).Matrix-free gravel and sand form inclusions into the bimodal gravel (Figure 2).Briefly, the architecture of the deposit is characterized by stratification with two blocks, namely the upper layer and the bimodal gravel with its inclusions (sand lens and free-matrix gravel).

Modeling: Numerical Options
Modeling was performed using the HYDRUS program [15].This code is a finite element model for simulating the movement of water, heat, and multiple solutes in variably saturated media.The program numerically solves the Richards equation as based on the Darcy-Buckingham law and mass conservation [23] for saturated-unsaturated water flow and Fickian-based advection dispersion equations for heat and solute transport.Additional models were also implemented to account for PF [23].In this study, HYDRUS was used to solve Equations (1-7) for the specific scenarios detailed below.
The trench was represented by a 2D flow domain of 13.5 m in length and 2.5 m in depth meshed through triangular elements with 200,000 as maximum number of nodes.All materials were geometrically represented according to the sedimentary description of the trench.Figure 2(b) illustrates the numerical domain that mimics the real deposit (Figure 2(a)).In order to spotlight the effect of heterogeneity, a uniform profile was built considering only the bimodal gravel that constitutes the major part of the deposit.
For flow, boundary conditions were set at no flux on the sidewalls, free drainage at the lower boundary and atmospheric boundary on the top.For the atmospheric boundary, the evaporation was fixed at zero and precipitation rate was chosen so as to mimic the targeted values of entering flow rates, i.e. 1.5 mm•h −1 and 15 mm•h −1 .These values correspond to medium rainfall intensities.For initial conditions, a first calculation is computed.The two profiles (uniform and heterogeneous) were assigned a water pressure head of hi = −0.01m to mimic hydric conditions close to saturation.Drainage and water redistribution were then modeled by imposing a zero flux at surface for one week (168 h).Water pressure heads obtained at the end of this drainage phase were considered as typical values encountered before rainfall events.These water pressure heads were then implemented as initial conditions into the modeling of the rainfall events.For hydraulic properties, each material was assigned the values proposed by [1] [2] who determined these values by infiltrometry (BEST method) [24] and by Arya and Paris method [25] for the coarser material (matrix-free gravel).Hydraulic parameters are listed in Table 1.
For solutes, we considered a third-type boundary.This kind of boundary condition controls solute concentrations and the quantity of solute that enters or exits the domain by mass balance considerations.For initial condition, we considered a zero value over the whole numerical domain.Solutes were injected at 1 L•g −1 which corresponds to the order of concentrations of major elements.Model solutes were considered with one tracer with no reactivity (K d = 0) and two reactive solutes with adsorption coefficients of K d = 10 −4 L•g −1 and K d = 10 −3 L•g −1 .These values are in the range between non-reactive solutes and values obtained for reactive solutes like copper in the same lithofacies [3].The comparison of the results for reactive and non-reactive solutes permits the identification of the separate effects convection, dispersion and retention.Finally, in agreement with previous studies [3], the longitudinal dispersivities of the materials were set at 2 cm for allcoarse materials and at 2 mm for the sand; and transversal dispersivities were fixed at one tenth of longitudinal dispersivities.ɬθr and θs (m 3 •m −3 ) denote the residual and saturated water contents, respectively; Ks (m•h −1 ) the saturated hydraulic conductivity; α (m −1 ) the scale parameter for water pressure head, n (−) the pore size distribution index and l the pore connectivity parameter (l = 0.5).

Effect of Heterogeneity on Water Flow
The evolution of volumetric water content is presented as a function of time for the two entering flow rates imposed on the top of the two profiles described above.In Figure 3, the evolution of water content is illustrated for two contrasted profiles: uniform (left) and heterogeneous (right).
In case of low velocity (Figure 3(a)), it is clear that in the uniform profile, water moves down to the low boundary condition along horizontal isolines, revealing a piston type infiltration.Wetting front progression can be assimilated to plug flow.In opposite, in heterogeneous sections, the geometry of the wetting fronts is influenced by inclusions.Indeed, water contents remain low in the sand lens and matrix-free gravel due to their lower water retention capacities.This results in a drastic drop in hydraulic conductivity and leads to flow deviations at their interfaces with the predominant bimodal gravel.There is a significant regionalization of flow and a patent barrier effect.This is usually explained by the fact that the fine-grained material overlay a coarse-grained material.The matrix potential at the textural interface is then so low that water cannot enter the coarser soil [26] [27].PF caused by capillary barriers is often referred to as funneled flow [28].In the studied case, such a flow pattern results also from the drainage phase.At the beginning of drainage, the whole section is saturated.Under this condition, the conductivity of bimodal gravel is low compared to other materials, and thus drains more slowly.Conversely, sand drains faster, leading to lower water contents (Figure 3(a)).Once unsaturated, sand lens will present a hydraulic conductivity lower than that of bimodal gravel and will therefore limit the infiltration, promoting water retention above and deviating flowaway.Matrix-free gravel also has the same capillary barrier effect, but this effect is more localized.This induces funneled flow.
In case of high entering flow rate (Figure 3(b)), the wetting front progression is more rapid.But, the same pattern is observed with a plug flow in the homogeneous profile and a funneled flow in the heterogeneous section.Yet, PF becomes less marked than for the lower entering flow rate.In particular, after 10 h infiltration (Figure 3(b)), water starts infiltrating the sand lens but keeps on avoiding inclusions.Under these conditions, sand, which used to form a capillary barrier, exhibits a higher hydraulic conductivity, becoming more and more permeable.Sand lens loses its capillary barrier effect.Conversely, the matrix-free gravel keeps its capillary effect and remains impermeable, thus forcing water to avoid it.
Clearly, these results show that the lithological heterogeneity induces a modification of flow pattern by leading to funneled flow.Such a flow pattern results from the more intense drainage of the sand lens and the matrix-free gravel during the drainage phase, which results in lower water content and hydraulic conductivity during the subsequent infiltration phase.This impact is enhanced for the lower entering flow rate, since this keeps the sand and the matrix-free gravel in a drier state, thus favoring the capillary barrier effect.
The lithological heterogeneity impacts also fluxes at the lower boundary.For both entering flow rates, there is earlier arrival of wetting fronts in the heterogeneous profile (Figure 4).In addition, for the lower infiltration rate, more dispersion is observed around the wetting front due to PF.These results show that the lithological heterogeneity impacts water fluxes at any depth of the soil profile and may thus impact groundwater recharge.After enough time, flow rates reach the value of the entering flow rate.We consider that steady state is reached since the free drainage boundary flux stabilizes and reaches the surface value and since the volume of water in the system stabilizes.When the steady state is established, the flow field differs significantly between the uniform and the heterogeneous profiles (Figure 5).The highest velocities are encountered at the center of funneled flow pathways, i.e. in the vicinity of interfaces between the sand lens and the bimodal gravel between the matrix-free gravel lenses and the bimodal gravel.For the lower entering flow rate, streamlines avoid completely the inclusions (i.e.all sand and gravel lenses), and the magnitude of velocity exceeds that of the uniform profile by a factor of more than 10.For the higher entering flow rates, flow regionalization is less marked, some of the streamlines pass through the sand lens and the difference between profiles in terms of velocity magnitude decreases.Clearly, steady state flow differs between the uniform and the heterogeneous profiles.The impact of such a difference on solute transfer is investigated in the following part.

Effect of Heterogeneity on Solute Transfer
Lithological heterogeneity markedly impacts the tracer breakthrough curves (BTC).For instance, at low entering flow rate (Figure 6(a)), we can notice the significant impact of PF on the tracer BTC with a faster breakthrough (around 75 h of outdistance) and a more important tail, in opposite to the symmetric and narrow peak for tracer BTC in the uniform profile.This asymmetry of tracer BTCs, typical for PF, reveals that some solutes arrive earlier (between 100 h and 150 h) and others stay longer in the section (exit between 300h and 450h).The concomitancy of PF pathways and stagnant zones explains such a pattern.The solutes that will pass through the PF pathways will have a short residence time and will exit early.On the contrary, the solutes that undergo diffusion between fast flow and stagnant zone will be delayed and exit later.Similar trends were obtained at the high entering flow rate with a more asymmetric tracer BTC in the heterogeneous profile (Figure 6(b)).Yet, the discrepancy between the BTCs is less marked.This results from the fact that flow pattern is much more heterogeneous at low than at high flow rate.With more numerous and faster PF pathways, flow field is strongly heterogeneous (Figure 5) and impacts drastically tracer BTCs.
Reactive solutes are transferred similarly to the tracer apart from the fact that they are trapped in the soil and consequently, their residence times are longer.Thus, reactive solutes are delayed in comparison with the tracer (Figures 6(c)-(  above the sand lens and matrix-free gravel lenses, due to slower flow in these zones.At high flow rate, heterogeneity effect remains but is less important.In fact, as shown in curves are less spread and asymmetric.The sand lens loses its effect on flow, and solutes remain mostly trapped above gravel lenses.This is corroborated by Figure 7 where the most reactive solute is represented at two intermediate times of its concentration evolution.Clearly, the PF explains the regionalization of retained solutes and the modeled regionalization is in agreement with field observations [1].

Discussion
Obtaining more patent PF for the lowest input flow rates seems counterintuitive.Indeed, it is usually considered that PF takes place at high pressure when highly conductive areas are activated due to high water pressure heads.This conceptual model is convenient to some configurations such as macropores or fractured soils.In this case, the soil is constituted of a matrix with a macropore and crack network.The activation of the crack and macropore network occurs when water pressure head is high enough to fill the network with water.This generates a PF.
In the studied case, the water pressure heads are not high enough to activate the most draining materials (i.e. the sand and the matrix-free gravel).The averaged water pressure heads are small enough to induce lower hydraulic conductivity for the sand and the matrix-free gravel in comparison with the main lithofacies (gravel).The inclusions of sand and gravel act as less permeable inclusions.This leads to the development of the capillary barriers observed.In fact, following the drainage phase, the most permeable materials (i.e. the sand and the gravel-free matrix) desaturate completely due to their low water retention capacity.The surrounding matrix becomes thus more conductive and serves as a transfer medium.The in clusions only induce a local deflection of flows.When these inclusions are slightly larger (e.g., sand lens), the deviation is large enough to induce a funneled flow without a true by pass of the matrix.
PF generates a typical pattern for solute transfer.It is obvious that PF allows some of the solutes to pass through the whole section quite quickly.This reduces the chance of these solutes to be trapped by the soil and disfavors their retention.This assumption has been demonstrated at several scales and using field observations [11] or laboratory column experiments [29] [30].A more homogenous flow means more pores accessible to flow and solutes.The increase in access to pores warranties the access of solutes to reactive surfaces of the pore walls.In addition, the increase in the water volume involved by flow means a decrease in pore velocity and thus an increase of the contact time between pollutants and the reactive surfaces of particles.These hypotheses are in full agreement with previous studies that pointed out an improvement of retention in relation with flow homogeneity for pollutants [18]- [29].Clearly, the homogeneity of flow must be targeted for an optimum retention of pollutant by the soil.

Conclusions
In this paper, we studied the effect of lithological heterogeneity on PF with regards to input flow rates and solute reactivity.For this purpose, a numerical study was performed on the basis of previous researches that offer a sedimentological description of the subsoil underneath the studied site.Numerical modeling permitted highlighting PF occurrence.According to these conditions, the time required reaching steady state over the whole soil profile, wetting fronts, water contents evolution during transient state and velocity vector distributions at steady state were accurately characterized and compared to those of a uniform section.For solutes, several distribution coefficients were also tested and compared to a tracer behavior.The numerical study revealed that lithological heterogeneity clearly favored PF, mostly at low input flow rates.Paradoxically, low input flow rate or dry hydric conditions are most conducive to the establishment of PF.Indeed, pressure head must be low enough to make the hydraulic conductivity of the most draining materials drop below that of the surrounding matrix.This sudden hydraulic conductivity decrease induces capillary barrier effects.Indeed, sand lens and matrix-free gravel lenses were the sources of capillary barrier effects, leading to a funneled flow and a strongly heterogeneous flow field.Lithological heterogeneity also affected the recharge of groundwater with earlier and more dispersed wetting fronts.Considering distribution coefficient, PF affects more the transfer of reactive solutes.
These results suggest that unsaturated flow favors PF, mostly for low water fluxes imposed at surface.That corroborates the conclusions made by other researches, which confirmed that PF was more likely to perform under unsaturated conditions, influencing relevantly pollutant transfer by enhancing retention far from PF zones [3]- [31].That can explain specific soil pollution pattern with pollutant accumulation at surface along with pollutant concentration at specific locations far from surface, as already observed for several heavy metals in a studied infiltration basin [20].These data are of great interest since they allow generating a conceptual and numerical model to better understand and predict the development of PF and its impact on pollutant transfer in highly heterogeneous deposits.These results can give indicative advice about the monitoring of infiltration basins.For instance, the fact that that PF is more important under unsaturated conditions and at low flow rates means that additional precautions must be considered for these conditions regarding soil and groundwater pollution underneath infiltration basins.

Figure 1 .
Figure 1.Geological settings of the region of Lyon (a) and location of the study site area (b) according to [3].

Figure 2 .
Figure 2. Lithological description (a); mesh and material distribution implemented in HYDRUS (b); and particle size distribution (cumulative fraction and frequency) of the four main lithofacies (c).

Figure 3 .
Figure 3.Time evolution of the volumetric water content, calculated with HYDRUS for the two imposed velocities: 1.5 mm•h −1 (a) and 15 mm•h −1 (b). 1024

Figure 4 .
Figure 4. Boundary flux at 2.5 m depth (groundwater recharge) for the imposed velocity of 1.5 mm•h −1 (a) and 15 mm•h −1 (b); data for the uniform (solid line) and heterogeneous profile (dashed line).(a) Low velocity; (b) High velocity.

Figure 5 .
Figure 5. Velocity vectors field at steady state time for the imposed velocity of 1.5 mm•h −1 (a) and 15 mm•h −1 (b), calculated with HYDRUS for uniform profile (left) and heterogeneous profile (right); steady state times are 300h and 40 h, respectively.
e) versus Figure 6(a)).Solute reactivity highlights the effect of PF on BTCs.The higher K d the more patent the impact of PF on BTCs (Figure 6(c) and Figure 6(e)).The solute resident concentrations are illustrated in Figure 7.There is a regionalization of pollution as observed in the heterogeneous section (Figure 7(a) and Figure 7(b)).The observed regionalization reveals zones where solutes are trapped preferentially, i.e.