J. Hydrol. Hydromech., 65, 2017, 3, 209–221
DOI: 10.1515/johh-2017-0004

The role of heterogeneous lithology in a glaciofluvial deposit on unsaturated
preferential flow – a numerical study
Erij Ben Slimene1, Laurent Lassabatere1*, Jiří Šimůnek2, Thierry Winiarski1, Remy Gourdon3
1

Université de Lyon, Laboratoire d’Ecologie des Hydrosystèmes Naturels et Anthropisés, LEHNA, UMR 5023 CNRS, ENTPE,
UCB-Lyon-1, Rue Maurice Audin, 69518 Vaulx-en-Velin, France.
2
Department of Environmental Sciences, University of California Riverside, Riverside, CA, USA.
E-mail: Jiri.Simunek@ucr.edu
3
Université de Lyon, INSA-Lyon, LGCIE-DEEP, 20 Av. Albert Einstein, F-69621 Villeurbanne cedex, France.
*
Corresponding author. E-mail: laurent.lassabatere@entpe.fr

Abstract: An understanding of preferential flow in the vadose zone is crucial for the prediction of the fate of pollutants.
Infiltration basins, developed to mitigate the adverse effects of impervious surfaces in urban areas, are established above
strongly heterogeneous and highly permeable deposits and thus are prone to preferential flow and enhanced pollutant
transport. This study numerically investigates the establishment of preferential flow in an infiltration basin in the Lyon
suburbs (France) established over a highly heterogeneous glaciofluvial deposit covering much of the Lyon region. An investigation of the soil transect (13.5 m long and 2.5 m deep) provided full characterization of lithology and hydraulicproperties of present lithofacies. Numerical modeling with the HYDRUS-2D model of water flow in the transect wasused to identify the effects of individual lithofacies that constitute the deposit. Multiple scenarios that considered different levels of heterogeneity were evaluated. Preferential flow was studied for several values of infiltration rates applied
after a long dry period. The numerical study shows that the high contrast in hydraulic properties of different lithofacies
triggers the establishment of preferential flow (capillary barriers and funneled flow). Preferential flow develops mainly
for low water fluxes imposed at the surface. The role of individual lithofacies in triggering preferential flow depends ontheir shapes (layering versus inclusions) and their sizes. While lenses and inclusions produce preferential flow pathways,
the presence of the surface layer has no effect on the development of preferential flow and it only affects the effective
hydraulic conductivity of the heterogeneous transect.
Keywords: Preferential flow; Heterogeneity; Lithology; Numerical modeling; Soil hydraulic properties; Boundary conditions; HYDRUS-2D.
INTRODUCTION
Field and laboratory studies have demonstrated that preferential flow in the vadose zone is more a rule than the exception(Köhne et al., 2009a). Preferential flow can be seen as unevenand rapid movement of water and solutes through certain zonesof the porous media, while any remaining water is quasistagnant. Preferential flow can be caused by multiple physicalor chemical factors such as hydrophobicity, macropores, orhydraulic instability. Similar types of flow have also been related to heterogeneous lithology, specifically in heterogeneous
soils and deposits that are formed by lithofacies with contrasting properties. Preferential flow is considered to be an
important process since it is expected to enhance the transportof nutrients, pollutants, and bacterial pathogens (e.g., Jarvis,
2007; Köhne et al., 2009a, 2009b) by bypassing soil reactive
particles and thus limiting the storage, filtering, and bufferingfunctions of soils (e.g., Clothier et al., 1992). It is crucial tounderstand flow processes that trigger the establishment ofpreferential flow and to identify hydraulic conditions prone topreferential flow.
Infiltration basins are highly suitable for studying the establishment of preferential flow. Infiltration basins are ecotechnological structures designed to mitigate the impact ofincreased stormwater runoff from urban areas by infiltrating
runoff water directly into the ground. To ensure their infiltration function, these structures are usually established over highly permeable and strongly heterogeneous soils. Preferential
flow can thus easily develop in these soils and impact pollutant

transport. When instrumented to monitor flow and transport,these sites can serve as a field laboratory for the study of preferential flow and its impact on solute transport.
Modeling preferential flow in infiltration basins is a special
case of modeling preferential flow in the vadose zone. Modeling of preferential flow has been conducted for many different
purposes in order to improve the understanding of the water
cycle and the fate of pollutants in soils (e.g., Allaire et al.,2009; Köhne et al., 2009a, 2009b). Many different models, suchas the dual-porosity and dual-permeability non-equilibrium
models (Šimůnek et al., 2003), have been developed to simulate
preferential flow and transport during various hydrologicalprocesses, including water infiltration (e.g. Lassabatere et al.,2014). Modeling approaches at the macroscopic scale are mostoften based on the use of the Richards equation and can consider the spatial variability of hydraulic and transport properties ofthe soil (e.g., Birkholzer and Tsang, 1997). For strongly heterogeneous deposits and soils, different lithofacies with different
shapes of inclusions and layers can be defined in a numericaltransport domain, and flow can be modeled using theknowledge of the hydraulic properties of the lithofacies (e.g.,
Coutinho et al., 2015). However, to the authors’ knowledge, nostudies have clearly identified the role of individual lithofacies,
inclusions, or layers in heterogeneous soils on preferential flow.
Also, their roles have never been studied as functions of theinfiltration rate at the soil surface or as functions of the degreeof saturation, although it is well known that the degree of saturation is a key factor for the establishment of preferential flow,such as funneled flow, induced by lithological heterogeneities.

UnauthenticatedDownload Date | 7/25/17 6:13 AM

209

Erij Ben Slimene, Laurent Lassabatere, Jiří Šimůnek, Thierry Winiarski, Remy Gourdon

In this paper, we model the establishment of preferential
flow during infiltration in an infiltration basin involving astrongly heterogeneous deposit made of lithofacies with contrasting hydraulic properties and different shapes (layers andinclusions). Preferential flow is modeled at the scale of a soiltransect 13.5 m long and 2.5 m deep, accounting for lithologicalheterogeneity and considering a fixed flow rate imposed at thesurface. The establishment of preferential flow is depicted forboth transient and steady state conditions. This work focuses on
the influence of individual lithofacies and their shapes (inclusions versus layers) on flow heterogeneity at several flow rates
imposed at the soil surface, thus mimicking rainfall events withdifferent rainfall intensities.

  h n 
S e (h ) = 1 +     hg  

MATERIALS AND METHODSStudied infiltration basin
The Django Reinhardt infiltration basin is located at Chassieu, in the suburbs of Lyon, France. It covers 13,775 m2 andreceives stormwater from an urban industrialized watershed of
185 ha (Winiarski et al., 2006). It was directly set up above theheterogeneous glaciofluvial deposit by mixing the top 50–80
cm of the soil (Goutaland et al., 2008; Winiarski et al., 2006).
The vadose zone is approximately 13 m deep (an average depth
to groundwater). Two trenches were dug in the infiltration
basin, with excavation being conducted by a mechanical shovel.
The excavation was followed by a sedimentological study
that provided the detailed description of the soil profile beneaththe soil surface. The sedimentological study (Fig. 1a) showedthat the glaciofluvial deposit is composed of four main lithofacies: (i) the upper layer that is a mixture of the soil matrix andgravel, (ii) a mixture of the soil matrix and gravel that has abimodal particle size distribution (PSD) and occupies most of
the deposit below the top 50–80 cm layer, (iii) large lenses ofsand, and (iv) smaller lenses of matrix-free gravel (Goutaland etal., 2008). These lithofacies are referred to as Gravel with matrix (or the upper layer), Bimodal Gravel, Sand, and Gravel,
respectively. These four materials were characterized withregards to their water retention and hydraulic conductivity
curves (Goutaland et al., 2013; Lassabatere et al., 2010) using
water infiltration experiments and the BEST method (Lassabatere et al., 2006). The lithofacies exhibit contrasting properties (Fig. 1b) and can be sorted with regards to their hydraulicconductivity (a decreasing order) and water retention (an increasing order): Gravel > Sand > Gravel with matrix > Bimodal
gravel. As expected, the lithofacies that exhibit the largest
values of the saturated hydraulic conductivity exhibit at thesame time the lowest water retention. The contrast in hydraulicproperties between Gravel and Sand, on the one hand, and the
predominant lithofacies Bimodal Gravel, on the other hand, is
expected to favor funneled flow (Walter et al., 2000).
Numerical model
Modeling of water flow was carried out using the HYDRUS2D model (Šimůnek and van Genuchten, 2008; Šimůnek et al.,2008, 2016) that numerically solves the generalized Richards
equation, which is based on the combination of the DarcyBuckingham law and the mass conservation equation (Hillel,1998).

∂θ= ∇.( K (θ ) ∇ H )
∂t

(1a)

H =z+h(1b)210

where θ (L3 L–3) is the volumetric water content, x (L) and z (L)are the spatial coordinates, h (L) is the water pressure head, t(T) is time, and K (L T–1) is the unsaturated hydraulic conductivity tensor, which depends on the water content, θ. The solution of equations (1) requires knowledge of the water retention
and hydraulic conductivity functions, h(θ) and K(θ), respectively, in addition to the knowledge of boundary and initial conditions. In this study, we used the water retention and hydraulicconductivity functions as estimated by Goutaland et al. (2013),
considering the van Genuchten relation (van Genuchten, 1980)
in conjunction with the Mualem capillary model (Mualem,1976):
−m

(

(2a)

)

m
1


K ( S e ) = K s S el  1 − 1 − S e m 



(2b)

1
n

(2c)

θ −θr
θs −θr

(2d)

m =1−

Se =

2

where Se (–) is the degree of saturation, hg (L), n (–), m (–), andl (–) are hydraulic shape parameters, θr (L3 L–3) and θs (L3 L–3)
are the residual and saturated volumetric water contents, respectively, and Ks (L3 T–1) is the saturated hydraulic conductivity. The lithofacies were assigned the hydraulic properties asshown in Fig. 1b.
The transect was represented by a two-dimensional flow
domain 13.5 m long and 2.5 m deep. The flow domain was
discretized using unstructured triangular elements with 18,500
nodes (Fig. 1a). The lower boundary was assigned a free drainage boundary condition and vertical sides were assigned no flux
conditions. The entire transect was initially set at a water pressure head of –0.1 m to ensure almost full saturation and thenallowed to drain for 168 h. The water pressure heads at 168 h
were then used as initial conditions for the infiltration experiment. The transect was subjected to the target infiltration rate at
the soil surface. Numerically, the upper boundary condition was
assigned an atmospheric boundary condition with no evaporation and constant precipitation as often considered when modeling water infiltration in soils (Šimůnek et al., 2008). Several
values of precipitation rates were imposed between 0.0025 and0.04 m/h, corresponding to a range of final hydrological states
between “very dry” and close to saturation. These values werelower than saturated hydraulic conductivities of lithofacies,
ensuring full infiltration of water into the profile and no runoff
at the surface. These values were selected in agreement withincoming water fluxes measured in the retention pond locatedupstream of the infiltration basin (data not shown).
The role of lithofacies on flow heterogeneity was investigated for the following scenarios with different levels of heterogeneity (Fig. 1c). A control run (Scenario 1) denoted as "Hom"assumed a homogeneous soil profile consisting of a predominant lithofacies, i.e., bimodal gravel. The following five scenarios considered different levels of heterogeneity by defining
additional soil layers and inclusions in the homogeneous transect. Scenario 2, denoted as "Layer," additionally considered asurface layer of gravel with matrix. Scenario 3, denoted as"All," additionally considered all lithofacies and heterogeneities.
UnauthenticatedDownload Date | 7/25/17 6:13 AM

The role of heterogeneous lithology in a glaciofluvial deposit on unsaturated preferential flow – a numerical study

Fig. 1. Graphical representation of evaluated scenarios: (a) Observed heterogeneity in a two-dimensional transect, (b) Soil hydrauliccharacteristic curves, (c) Evaluated scenarios with different types of heterogeneities.

Scenario 4, denoted as "Incls," additionally considered bothsandy and matrix-free gravel inclusions. This scenario differsfrom Scenario 3 by not considering the surface layer.
Additionally, Scenarios 5 and 6, denoted as "Sand" and"Gravel," alternatively considered either the sandy or matrixfree gravel inclusions, respectively. A comparison of results
obtained for the homogeneous scenario (Scenario 1) with thosefor the heterogeneous scenarios (Scenarios 2 through 6) allows
for the discrimination of separate effects of the surface layer
and individual inclusions (sand lens, gravel lenses, and alllenses simultaneously) and of the overall effect of allheterogeneities together.
Analysis of the numerical data
The arrival of the water front at the lower boundary of thesoil profile and the spatial variability of water fluxes at thesteady state were compared between scenarios. To simplify thecomparison, we defined both several characteristic times toquantify the impact of lithological heterogeneities on the arrivalof the water front at the lower boundary, and dispersion indicators to characterize the spatial heterogeneity of the water flux
field at steady state. The following characteristic times were
defined (see the Result section for graphical illustration of thesetimes): a) breakthrough time (tbr), when the flow rate at thebottom boundary starts increasing (a transition between thedrainage phase and the arrival of the wetting front at the lower
boundary), b) wetting front arrival time (twf), when the flow rate
reaches the average value between initial and final flow rates,
and c) stabilization time (tstab), when the flow rate reaches aplateau approximately corresponding to the value of the flowrate imposed at the surface (i.e., the beginning of the steadystate). The value of twf can be considered as the average arrival
time of the wetting front at the lower boundary, and the difference between stabilization and breakthrough times (Δt =(tstab–tbr)/2) can be viewed as an indicator for time dispersion
around the averaged arrival time. These characteristic times canbe viewed as measures of water residence time in the systems.
The characteristic times allow for the quantification of theinfluence of different lithological heterogeneities on the flow
patterns and their variation depending on applied flow rates.
In addition, the flow resistance of a system may also be a
key parameter, and the effects of inclusions and layering may
be important when the hydrological response of such deposits
needs to be quantified. The flow resistance was derived below
from the analysis of pressure head profiles at steady states. To
quantify the effective hydraulic conductivity at the scale of atransect, we computed the hydraulic gradient over the entire
transect from the average water pressure heads at the upper andlower boundaries. We then derived the effective hydraulicconductivity using the Buckingham–Darcy law:

K eff =

( L + (hqLz =0

− hz = L ))(3)

where Keff (L T–1) is the effective hydraulic conductivity, L isthe depth of the profile (L), q (L T–1) is the infiltration rate
imposed at the surface, and hz = 0 (L) and hz = L (L) refer to theaverage water pressure heads at the surface and the lower
boundary, respectively. The values of effective hydraulic conductivities were then compared between scenarios.
UnauthenticatedDownload Date | 7/25/17 6:13 AM

211

Erij Ben Slimene, Laurent Lassabatere, Jiří Šimůnek, Thierry Winiarski, Remy Gourdon

Fig. 2. Time evolution (at 0, 25, 50, and 75 h) of water contents (top) and water fluxes (bottom) for the six evaluated scenarios at the lowest
flow rate (0.0025 m/h). Capital letters in the figure are referred to in the text.

RESULTS
Preferential flow at the lowest flow rate (0.0025 m/h)
Simulated volumetric water contents and water fluxes arepresented in Fig. 2 as a function of time for the first flow rate of
0.0025 m/h. Wetting fronts for the homogeneous and layered
profiles are typically piston-like and horizontal. The water flux
field is uniform with similar fluxes in the entire transect above
the wetting front. In contrast, for all the other scenarios, wetting
fronts have singular and complex shapes reflecting the locationsof the various inclusions that water attempts to bypass. Water
contents are lower in inclusions (e.g., letter A in Fig. 2), andwater accumulates above inclusions (e.g., letter B in Fig. 2),
which is a common capillary barrier effect. Water fluxes showthat flow pathways are distorted and mostly develop between
inclusions, leading to funneled flow. Locally, along preferential
flow pathways, the water flux can increase more than ten times

212

in comparison to the value imposed at the soil surface (e.g.,
letter C in Fig. 2., bottom).
Water fluxes at the lower boundary are depicted in Fig. 3aand the definition of related characteristic times is given in Fig.
3b. During the transient state, water fluxes at the lower boundary increase from a very low initial value (the end of the preceding draining phase) to the value of the flow rate imposed at thesoil surface (Fig. 3a). The increase does not occur at the sametime and depends on the scenario. There are two distinct responses. There is an earlier arrival of the wetting front for all
heterogeneous scenarios that include the sand lens, and a latearrival of the wetting front for all other scenarios (i.e., without
the sand lens) (Fig. 3a). This is reflected by lower values ofbreakthrough and arrival times (Figs. 3c and 3e) for scenarioswith the sand lens. The stabilization time and time dispersion
are increased by any kind of lithological heterogeneity, except
for the scenario with the surface layer that has no clear impact
UnauthenticatedDownload Date | 7/25/17 6:13 AM

The role of heterogeneous lithology in a glaciofluvial deposit on unsaturated preferential flow – a numerical study

Fig. 3. (a) Water fluxes at the bottom boundary, (b) definition of characteristic times, (c) breakthrough times tbr, (d) stabilisation times tstab,(e) wetting front arrival times twf, and (f) time dispersions Δt for the six evaluated scenarios at the lowest flow rate (0.0025 m/h).

in these characteristic times (Figs. 3d and 3f). It can beconcluded that the surface layer has no significant effect on
characteristic times, while the gravel inclusions increase thedispersion in arrival times, and the sand inclusion at the sametime decreases and disperses arrival times.
The water flux field at steady state is shown for the scenario
with all lithofacies in Fig. 4a. It is clear that the sand lens andthe gravel lenses induce deviations in flow and behave as capillary barriers at this low flow rate. This has a clear impact on theflux spatial heterogeneity and horizontal distributions of water
fluxes at all depths (Figs. 4b–g). Water fluxes are quite homogeneous for depths of 0 and 0.5 m for all scenarios. For depths
larger than 1.5 m, the effects of various inclusions become
more important. The sand lens clearly produces both a decreasein water fluxes inside the sand lens and an increase of them at
its right side due to flow diversion (e.g., letter A in Fig. 4f). The
gravel inclusions produce clear oscillations in water fluxes atthe same depth (e.g., letters C to E in Figs. 4d–g). In the middle, the effects of one gravel lens and the sand lens overlap,
leading to the superposition of respective effects (see letter B in
Figs. 4e–g).
Coefficients of variation CVq (the standard deviation divided
by the mean) were calculated for all water flux distributions
both at different depths and for different scenarios (Figs. 4b–g)
and then plotted as a function of depth in Fig. 4h. The coefficients of variation were also averaged to obtain a coefficient ofvariation at the transect scale and to define a global indicator of
the flux spatial heterogeneity for different scenarios (Fig. 4i).
Clearly, the highest spatial heterogeneity of water fluxes was
obtained for the heterogeneous transect containing all lithofacies (Fig. 4i). The sand and gravel lenses cause a spatial variability in the water fluxes to a lesser extent. The surface layer
has no effect on the spatial variability of the water fluxes,
which are as homogeneous over the length of the transect as in
the homogeneous transect.
The surface layer in our system does not seem to have an
important effect either on water residence time in the system orwater flux distributions. The sand lens triggers the main preferential flow pathway at its side, resulting in an increase in thespatial variability of water fluxes and at the same time, both adecrease in mean residence time and an increase in time dispersion. While the gravel inclusions have a substantial impact on
the spatial heterogeneity of water fluxes, they only have a limited impact on water residence time in the system. The geometry of inclusions seems to be a key factor. Small inclusions,randomly distributed in the soil profile, will produce flow
pathways that will bypass some inclusions before encounteringothers. Water will then be successively accelerated and slowed
down along the flow pathways. This will result in similar residence times and a larger time dispersion than in the homogeneous profile, as observed in the simulation results presented inFig. 3.
UnauthenticatedDownload Date | 7/25/17 6:13 AM

213

Erij Ben Slimene, Laurent Lassabatere, Jiří Šimůnek, Thierry Winiarski, Remy Gourdon

Fig. 4. (a) Water fluxes for the steady state flow conditions for a scenario with full heterogeneity at the lowest flow rate (0.0025 m/h), (b–g)
water flux profiles at different depths (0, 0.5, 1.0, 1.5, 2.0, 2.5 m, respectively), (h) corresponding coefficients of variation CVq at different
depths, (i) profile averaged values of coefficients of variation CVq (i).

Preferential flow at the highest flow rate (0.04 m/h)
In comparison to the simulations with the lowest flow rate,
the differences in the movement of the wetting fronts between
the homogeneous and heterogeneous scenarios at the highest
evaluated flow rate are smaller during the transient state. In all
the scenarios, wetting fronts are more horizontal, and water
infiltration resembles more a piston-like process. For instance,
at 4 and 6 h, wetting fronts are almost horizontal and indicate amore piston-like and less funneled flow than observed at thelower flow rate (e.g., letter A in Fig. 5 versus letter D in Fig. 2).
However, in the scenarios that consider the sand lens, wetting
fronts are curved below the sand lens, indicating that water
moves faster within the sand lens than in the rest of the profile(e.g., letter B in Fig. 5). This indicates a preferential water214

infiltration within the body of the sand lens. Similar conclusions can be obtained from water flux fields (Fig. 5c).
Wetting fronts arrive at the lower boundary earlier in all
scenarios that consider the sand lens than in the scenarios without the sand lens, irrespective of what other heterogeneities areincluded (Fig. 5b). The wetting fronts arrive later in scenarioswith the homogeneous or layered soil profiles. The scenariothat considers only the gravel inclusions produces intermediatearrival time. Similar trends can be observed with respect tocharacteristic times. All the scenarios containing the sand lenshave shorter wetting front arrival times and breakthrough times(Table 1). The scenarios with homogeneous and layered soil
profiles exhibit longer characteristic times, while the scenario
with the gravel inclusions produces intermediate values. The
earlier arrival of wetting fronts at the lower boundary in scenarios
UnauthenticatedDownload Date | 7/25/17 6:13 AM

The role of heterogeneous lithology in a glaciofluvial deposit on unsaturated preferential flow – a numerical study
Table 1. Characteristic times quantifying the impact of lithological heterogeneities on the arrival of the water front at the lower boundary,i.e., tbr - breakthrough time, tstab - stabilization time, twf - wetting front arrival time, Δt = (tstab–tbr)/2, and the coefficients of variation forwater flux distributions at the lower boundary, CVq, at the lowest and the highest imposed flow rates.
VariableHomtbr (h)tstab (h)twf (h)Δt (h)CVq (–)

36.4
85.4
59.4
24.5
0.02%

tbr (h)tstab (h)twf (h)Δt (h)CVq (–)

6.70
11.10
59.4
2.20
0.01%

Layer

All
Incls
Sand
The lowest flow rate (0.0025 m/h)
35.4
24.4
21.9
20.9
81.4
96.9
102.2
198.7
57.4
52.4
53.4
54.4
23
36.3
40.5
38.5
0.37%
66.4%
66.0%
51.6%
The highest flow rate (0.04 m/h)
6.65
5.41
4.88
4.89
11.15
10.10
9.80
10.85
57.4
52.4
53.4
54.4
2.25
2.35
2.46
2.98
0.08%
12.0%
8.97%
3.16%

Gravel
33.4
198.2
60.4
34.7
46.7%
5.92
10.35
60.4
2.22
9.79%

Fig. 5. Time evolution (at 0, 25, 50, and 75 h) for the six scenarios of water contents (a), water fluxes at the lower boundary (b), andvelocity fields at steady state for the homogeneous and fully heterogeneous scenarios at the highest flow rate (0.04 m/h). Capital letters in
the figure are referred to in the text.

with the sand lens is attributed to the enhancement of downward flow within the sand lens. All the scenarios exhibit similarvalues with respect to time dispersion (Table 1), indicating thatthe wetting fronts simply moved backward or forward along atime axis with no change in dispersion around the averagearrival time. Although the time interval during which water
fluxes increase appears to be narrower in scenarios withthe homogeneous or layered soil profiles (Fig. 5b), this is not
UnauthenticatedDownload Date | 7/25/17 6:13 AM

215

Erij Ben Slimene, Laurent Lassabatere, Jiří Šimůnek, Thierry Winiarski, Remy Gourdon

Fig. 6. Comparison of flow patterns in the vicinity of the gravel and sand lenses at the two extreme flow rates (0.0025 m/h (left) versus 0.04
m/h (right)). The color scales in (a) and (b) are the same as in Fig. 2 and 5, respectively.

reflected in the values of time dispersion Δt, which does notreflect this observation (Table 1).Scenarios with the surface soil
layer and both sand and gravel inclusions showed that theseheterogeneities have only a slight impact on the water flux field
at steady state at the largest evaluated flux (e.g., Fig. 5c andFig. 6). The sand lens has lost its effect of a capillary barrier atthis flux rate, as revealed by the fact that flow no longer avoids
the sand lens (Fig. 6, letters C and D versus A and B). In thiscase, the water pressure head and the corresponding hydraulicconductivity of the sand is high enough, compared to the surrounding bimodal gravel, to prevent the capillary barrier effectand to allow water to directly enter the sand lens. The sand lens,rather than diverting flow as occurred at low fluxes, preferentially conducts water through its body at higher fluxes (evenmore visible in the left part of the lens). On the other hand, thegravel lenses continue to act as capillary barriers, but to asmaller extent than at the lowest flow rate (Fig. 6, letters C andD versus A and B). This change in the flow pattern also has an
impact on the horizontal distribution of water fluxes at the scaleof the transect (Fig. 5c). Distributions of water fluxes are relatively uniform in the scenario with only the sand lens and in thehomogeneous and layered scenarios. The sand lens no longer
distorts flow pathways. The scenarios containing the gravel
lenses exhibit more water flux variations. This is confirmed by
the values of CVq (Table 1). This last parameter indicates thatthe flow field appears more heterogeneous for all scenarios that216

consider gravel inclusions (Table 1, scenarios All, Incls, andGravel, The highest flow rate). The numerical simulations performed at the highest flow rate indicate that sand loses its capillary barrier effect and that gravel remains the sole source of
spatial heterogeneity. The spatial heterogeneity of flow issmaller at the higher flow rate than at the smaller flow rate at
the spatial scale of these simulations, with values of the coefficients of variation around 10% at a flux of 0.04 m/h, versus
60% at a flux of 0.0025 m/h (Table 1).
It can be concluded that at the highest flow rate, the capillarybarrier effect disappears for the sand lens and decreases for thegravel lenses. An increase in the applied surface flow rate lessens the capillary barrier effect because it increases the pressure
heads in the soil profile and reverses the difference in corresponding hydraulic conductivities of the lenses and surroundingsoil, which actually favors the moistening of inclusions. Consequently, flow becomes much more homogeneous and water
residence times are less dispersed at the scale of simulations.
The sand lens conducts a part of infiltrating water through its
body, leading to a lower residence time in any scenario that
contains the sand lens. In contrast, a similar decrease in residence time was due to preferential funneled flow at the lowest
flow rate. The results of these simulations demonstrate thatflow heterogeneity is boosted at the lowest flow rate and thatthe lithological heterogeneity loses its potential for disturbing
flow at the highest flow rate.
UnauthenticatedDownload Date | 7/25/17 6:13 AM

The role of heterogeneous lithology in a glaciofluvial deposit on unsaturated preferential flow – a numerical study

Fig. 7. Time enhancement (Rtwf) as a function of the lithological heterogeneity (a) and the flow rate (b), dimensionless time dispersion (RΔt)
as a function of the lithological heterogeneity (c) and the flow rate (d), and the coefficient of variation (CVq) for water fluxes as a functionof the lithological heterogeneity (e) and the flow rate (f).

Preferential flow as a function of flow rate – at intermediate
fluxes
The effect of lithological heterogeneities on flow is further
investigated as a function of flow rates. Fig. 7 depicts characteristic times in the form of dimensionless indicators: (i) timeenhancement (Rtwf), which is defined as the arrival time of a
wetting front (twf) for a particular scenario divided by a corresponding value for the homogeneous scenario, (ii) the relative
time dispersion (RΔt), which is defined as the time dispersion(Δt) divided by the arrival time of the wetting front (twf), and(iii) the coefficient of variation for the spatial heterogeneity ofwater fluxes at the transect scale (CVq). These parameters areshown as a function of lithological heterogeneities (Fig. 7, leftside) and flow rates imposed at the soil surface (Fig. 7, right side).
With respect to the role of lithofacies, the main trends become apparent when results for all flow rates are simultaneously evaluated. The surface layer alone has an insignificant effecton water fluxes, which is similar to the homogeneous system(Figs. 7a, c, e), indicating that flow homogeneity is similar inboth systems. The gravel lenses only have a small effect ontime enhancement (Fig. 7a) and relative time dispersion (Fig. 7c),
UnauthenticatedDownload Date | 7/25/17 6:13 AM

217

Erij Ben Slimene, Laurent Lassabatere, Jiří Šimůnek, Thierry Winiarski, Remy Gourdon

but a significant effect on flow spatial heterogeneity (Fig. 7e).
Despite their strong effect at the local scale (inducing funneledflow), they have no clear effect at the transect scale. This is theresult of a random distribution of inclusions that distort flow
pathways more than they boost downward preferential flow. In
contrast, the sand lens has a significant effect on time enhancement (Fig. 7a) and relative time dispersion (Fig. 7c). It
also significantly affects the spatial heterogeneity of flow, but
to a lesser extent than the gravel lenses (Fig. 7e). The sandinclusion induces a more structured preferential flow pattern,
resulting in a more significant effect on water residence time at
the transect scale. In other words, the geometry and the characteristic length of inclusions are the key parameters at the transect scale.
The impact of each lithological heterogeneity also depends
on the flow rate. This can be demonstrated by comparing theeffects of the sand and gravel lenses at the low and high flow
rates of 0.0025 and 0.04 m/h, respectively. The gravel lenses
lose their capillary barrier effect with an increase in the flow
rate. The sand lens even transforms from a capillary barrier to amore conductive porous media. This shift in the hydraulicbehavior is reached at a threshold flux value of 0.03 m/h. For
the sand lens, there is no visible trend in time enhancement(Fig. 7b, Sand). The sand lens boosts the downward movementof the wetting front by creating preferential flow pathways atflow rates lower than 0.03 m/h or by conducting water throughits body at flow rates higher than 0.03 m/h. However, in allcases, the effect on the wetting front arrival is roughly the same
(Fig. 7b, Sand). For the gravel lenses, the enhancement of thedownward movement of the wetting fronts is more important atlarger flow rates (Fig. 7b, Gravel). Both sand and gravel inclusions lose their effect on time dispersion with increasing flow
rates, probably as a result of diminishing capillary barrier effects (Fig. 7d, Gravel or Sand). The same conclusion can bemade for spatial flow variability, which decreases when theflow rate is increased (Fig. 7f). An increase in the flow rate
reduces spatial flow variability to the level observed for the
scenario with the homogeneous profile (Fig. 7f, heterogeneous
scenarios versus Hom). Flow clearly becomes much more homogeneous when the flow rate is increased. Given a gradualdecrease of relative time dispersion (Fig. 7d) and spatial flux
heterogeneity (Fig. 7f) with an increase in the flow rate, it canbe concluded that the shift in the flow regime (homogeneous
versus heterogeneous) is continuous.
The homogenization of flow with the flow rate may be counterintuitive. It is usually assumed that higher flow rates and thus
higher water pressure heads in heterogeneous soils will activate
flow through zones with higher saturated hydraulic conductivities (Birkholzer and Tsang, 1997; Vogel et al., 2000). This isindeed often the case for soils with macropores, fractures, orcracks (Gerke and van Genuchten, 1996). In such cases, preferential flow is usually observed under conditions with highsurface pressure heads, such as under ponded conditions at thesoil surface (e.g., Gupta et al., 1999). In the studied deposit, thelithofacies with very high saturated hydraulic conductivities, inparticular the gravel lenses, form inclusions that do not constitute a continuous network. As a result, water cannot flow in awell-connected network, bypassing the soil matrix and moving
quickly from the soil surface to the bottom of the profile, which
explains why flow remains mostly homogeneous under saturated conditions and high flow rates. In contrast, the main sourceof preferential flow in the studied deposit is linked to the capillary barrier effects and is thus more apparent at low flow rates
under unsaturated conditions. Under these conditions, gravel
lenses act as capillary barriers, thus triggering funneled flow.

218

However, their impact remains local, as a result of the small
characteristic lengths of these inclusions. The sand lens also
induces funneled flow. However, since its characteristic lengthis larger, the sand lens induces a flow deviation large enough tolead to funneled flow without a real bypass of the soil matrix.
The combination of different lithofacies does not alter thecontribution of individual lithofacies. The spatial flow heterogeneity seems to be induced mainly by gravel at the local scale.
Indeed, the values of CVq for the scenarios involving both gravel and sand lenses have similar CVq values as the scenario containing only the gravel inclusions (Fig. 7e, All and Incls versus
Gravel). In contrast, the scenarios ‘Incls’ and ‘All’ exhibit
similar values of time enhancement as the scenario with thesand lens alone (Fig. 7a). In other words, the sand lens is themain lithofacies contributing to the enhancement of the downward movement of the wetting fronts at this spatial scale, evenwhen the gravel lenses are included. The effects of the sand andgravel lens on time dispersion are additive. Indeed, the relative
time dispersion is higher for the scenario that considers all thelenses compared to the scenario that considers only the sandlens (Fig. 7c, All and Incls versus Gravel).
Finally, the addition of the surface layer does not change either spatial flux variability, time enhancement, or time dispersion (Fig. 7a, Fig. 7c, Fig. 7e, All versus Incls). It is interesting
to note that the addition of a surface layer does not change flow
patterns in lower horizons and does not affect flow homogeneity. Under other circumstances, several authors found a significant effect on flow homogeneity from the included soil layers(e.g., Lamy et al., 2013; Lassabatere et al., 2004). These authors
showed that the inclusion of a geotextile into heterogeneous
soils or soils containing macropores could interrupt preferential
flow. However, these layers were placed in the middle of highly
conductive zones (e.g., macropores) and had a lower hydraulicconductivity than the conductive zones. In the studied deposit,
the surface layer is located above the zones where preferential
flow develops and does not have necessarily a lower hydraulicconductivity, depending on the value of flow rate imposed atsurface (Fig. 1).
Influence of lithological heterogeneities on flow resistance
The results presented above give some insight on the effectsof heterogeneous lithology on flow patterns, heterogeneity ofwater fluxes, and wetting front arrival times. The analysis ofnumerically simulated pressure head profiles provides an insight on the effects of lithological heterogeneities on flow
resistance. In the case of a homogenous soil profile, the steadystate water pressure head will be the same over the entire
profile (Fig. 8a), resulting in a unit hydraulic gradient (Figs.
8c–d, Hom). Equation (3), along with the hypothesis of
uniformity of water pressure heads between boundaries( h z = 0 ≈ h z = L ), shows that Keff is equal to q. The steady-statewater pressure head hsst will be equal to the value that satisfies
the equation K(hsst) = q, where K(h) refers to the hydraulicconductivity function defined by equations (2), and q is theapplied surface flux. These assumptions are valid provided thatthe infiltration rate imposed at the soil surface does not exceed
the saturated hydraulic conductivity, which is the case in thisstudy. Otherwise, the steady-state infiltration rate will be equal
to the saturated hydraulic conductivity and water runoff or
ponding will occur at the soil surface When the transects are nolonger homogeneous, the water pressure heads at boundaries
will no longer be uniform (Fig. 8b). The sand and gravelinclusions cause a slight increase in the water pressure head at

UnauthenticatedDownload Date | 7/25/17 6:13 AM

The role of heterogeneous lithology in a glaciofluvial deposit on unsaturated preferential flow – a numerical study

Fig. 8. The steady-state pressure head profiles for uniform (a) and heterogeneous (b) scenarios, the hydraulic gradient (i) as a function of
scenarios (c) or flow rates (d), and the effective hydraulic conductivity of the heterogeneous scenarios divided by that of the homogeneous
scenario as a function of scenarios (e) or flow rates (f).

the soil surface, resulting in hydraulic gradients larger than 1
for all flow rates imposed at the soil surface (Figs. 8c–d, Incls,
Sand and Gravel). This indicates that the gravel and sandinclusions slightly impede flow and produce an effective
hydraulic conductivity (profile-averaged) smaller than for theuniform profile (Figs. 8e–f, Incls, Sand and Gravel versus
Hom). For the layered system, the hydraulic gradient depends
strongly on the flow rate. For flow rates below a thresholdvalue of 0.035 m/h, the effective hydraulic conductivity of thelayered system is significantly lower than that of thehomogeneous system (Figs. 8e–f, Layer versus Hom),producing a hydraulic gradient higher than 1 (Figs. 8c–d, Layer
versus Hom). Otherwise, the effective hydraulic conductivity of
the layered system exceeds that of the homogeneous system(Fig. 8f, Layer and q = 0.04), leading to a hydraulic gradientbelow unity (Fig. 8d, Layer and q = 0.04) and a lower waterpressure head at the soil surface. This shift results from the factthat at low flow rates, i.e., at low water pressure heads, thesurface layer has a smaller hydraulic conductivity than the
predominant lithofacies, whereas at high flow rates close tosaturation, the situation is opposite (Fig. 1b). Finally, for the
fully heterogeneous soil profile, the effective hydraulicconductivity is mainly impacted by the surface layer. Indeed,similar values are found for the effective hydraulic conductivity
and hydraulic gradients between the fully heterogeneous
systems and the layered system (Figs. 8c–d and Figs. 8e–f,
Layer versus All). It can be concluded for studied systems thatthe surface layer is the predominant lithofacies that defines theeffective hydraulic behavior of the deposit and that the graveland sand inclusions only have a slight effect.
CONCLUSIONS
The results of this study support the following conclusions
concerning the role of lithological heterogeneities on
preferential flow in a strongly heterogeneous glaciofluvial
UnauthenticatedDownload Date | 7/25/17 6:13 AM

219

Erij Ben Slimene, Laurent Lassabatere, Jiří Šimůnek, Thierry Winiarski, Remy Gourdon

deposit. The soil stratification mainly affects the equivalentpermeability of the heterogeneous system. It increases theeffective hydraulic conductivity for conditions close tosaturation and lowers it for unsaturated conditions. It is
expected that under saturated conditions, the surface layer,
which is more conductive than underlying layers at saturation,
should not adversely affect water infiltration at the soil surface
or even increase water infiltration into the soil profile, thus
reducing the risk of water runoff or ponding at the soil surface.
Despite its effect on the effective hydraulic conductivity, thesurface layer has no effect on flow pathways and flow
heterogeneity.
Various sand or gravel inclusions only have a negligible
effect on the effective hydraulic conductivity of theheterogeneous systems. Therefore, they should not affectinfiltration of surface water or overall water fluxes. Inclusions
will only impact the flow patterns (homogeneous orheterogeneous flow patterns
with channeling). For all inclusions, the main mechanism
responsible for the establishment of preferential pathways isfunneled flow, which is more a result of the lower hydraulicconductivities of inclusions than of the surrounding soil for aparticular flow rate. The size of inclusions also has animportant influence. Large inclusions (e.g., sandy lenses) will
favor the development of well-organized and structured
channeling, impacting water flow across the system. However,smaller inclusions will induce only local preferential flow,which will have less impact on the mass transfer at larger
scales.
The degree of flow heterogeneity is also a function of thewater content or the flow rate imposed at the surface. Dry
conditions and low flow rates are the most conducive for theestablishment of preferential flow for the systems evaluated inthis study. In fact, these conditions cause the coarser lithofacies
to be drained and to have a significant drop in hydraulicconductivity, which will cause capillary barrier phenomena inthe system. Gravel inclusions have a hydraulic conductivity
lower than the predominant lithofacies under all circumstances.
However, this effect decreases for the sand lens when the flowrate increases. Therefore, the sand lens is transformed from aless permeable inclusion that diverts flow at lower fluxes to amore permeable inclusion that produces homogeneous flow at
higher fluxes.
This numerical study clearly shows that the geometry oflithofacies will determine their impact on hydrological behaviorof the transect. When the lithofacies constitute a continuous
layer, they may more greatly impact the effective hydraulicconductivity of heterogeneous deposits. Conversely, thelithofacies that form small inclusions will more greatly impact
flow heterogeneity, but are not expected to impact water fluxes
or the overall water balance. It should be noted that theseeffects were observed for a specific case of one transect in a
glaciofluvial deposit and for a set of specific hydraulicconditions. However, the methodology presented here could beapplied to other transects with the spatial variability oflithological heterogeneities on the scale of the deposit(Coutinho et al., 2015; Goutaland et al., 2013) and to othertypes of heterogeneous deposits and soils. This study could also
be expanded by investigating the dependency of the effects oflithological heterogeneities on rainfall event characteristics andon a succession of multiple rainfall events. Further research isneeded to understand the effects of such modification of waterflow on the fate of pollutants. However, inclusions are expected
to have a dramatic impact on solute transport due to theireffects on preferential flow (Winiarski et al., 2013). Indeed, it

220

has often been demonstrated that preferential flow will decrease
the contact between solute and reactive particles of the soils,
thus reducing the potential for pollutant sorption into the soil
(e.g. Köhne et al., 2009a, b; Lamy et al., 2013; Lassabatere etal., 2004, 2007). Inclusions may thus be a factor in enhancing
the transport of pollutants towards groundwater, depending on
their reactivity.
REFERENCES
Allaire, S.E., Roulier, S., Cessna, A.J., 2009. Quantifying
preferential flow in soils: A review of different techniques. J.
Hydrol., 378, 179–204.
Birkholzer, J., Tsang, C.-F., 1997. Solute channeling in
unsaturated heterogeneous porous media. Water Resour.
Res., 33, 2221–2238.
Clothier, B.E., Kirkham, M.B., McLean, J.E., 1992. In Situ
Measurement of the Effective Transport Volume for Solute
Moving Through Soil. Soil Sci. Soc. Am. J., 56, 733–736.
Coutinho, A.P., Lassabatere, L., Winiarski, T., Cabral, J.J.S.P.,
Antonino, A.C.D., Angulo-Jaramillo, R., 2015. Vadose Zone
Heterogeneity Effect on Unsaturated Water Flow Modelingat Meso-Scale. J. Water Resour. Prot., 7, 353–368.
Gerke, H.H., van Genuchten, M.T., 1996. Macroscopic
representation of structural geometry for simulating waterand solute movement in dual-porosity media. Adv. Water
Resour., 19, 343–357.
Goutaland, D., Winiarski, T., Dubé, J.-S., Bièvre, G.,
Buoncristiani, J.-F., Chouteau, M., Giroux, B., 2008.
Hydrostratigraphic characterization of glaciofluvial deposits
underlying an infiltration basin using ground penetratingradar. Vadose Zone J., 7, 194–207.
Goutaland, D., Winiarski, T., Lassabatere, L., Dubé, J.S.,
Angulo-Jaramillo, R., 2013. Sedimentary and hydrauliccharacterization of a heterogeneous glaciofluvial deposit:
Application to the modeling of unsaturated flow. Eng. Geol.,
166, 127–139. http://dx.doi.org/10.1016/j.enggeo.2013.09.006
Gupta, A., Destouni, G., Jensen, M.B., 1999. Modelling tritiumand phosphorus transport by preferential flow in structured
soil. J. Contam. Hydrol., 35, 389–407.
Hillel, D., 1998. Environmental Soil Physics: Fundamentals,
Applications, and Environmental Considerations. Academic
Press.
Jarvis, N., 2007. A review of non-equilibrium water flow andsolute transport in soil macropores: Principles, controlling
factors and consequences for water quality. Eur. J. Soil Sci.,
58, 523–546.
Köhne, J.M., Köhne, S., Šimůnek, J., 2009a. A review of model
applications for structured soils: a) Water flow and tracer
transport. J. Contam. Hydrol., 104, 4–35.
Köhne, J.M., Köhne, S., Šimůnek, J., 2009b. A review of model
applications for structured soils: b) Pesticide transport. J.
Contam. Hydrol., 104, 36–60.
Lamy, E., Lassabatere, L., Bechet, B., Andrieu, H., 2013. Effectof a nonwoven geotextile on solute and colloid transport inporous media under both saturated and unsaturated
conditions.
Geotext.
Geomembr.,
36,
55–65.
http://dx.doi.org/10.1016/j.geotexmem.2012.10.009
Lassabatere, L., Winiarski, T., Galvez-Cloutier, R., 2004.
Retention of three heavy metals (Zn, Pb, and Cd) in a
calcareous soil controlled by the modification of flow with
geotextiles. Environ. Sci. Technol., 38, 4215–4221. DOI:
10.1021/es035029s.
Lassabatere, L., Angulo-Jaramillo, R., Soria Ugalde, J.M.,
Cuenca, R., Braud, I., Haverkamp, R., 2006. Beerkan
UnauthenticatedDownload Date | 7/25/17 6:13 AM

The role of heterogeneous lithology in a glaciofluvial deposit on unsaturated preferential flow – a numerical study

estimation of soil transfer parameters through infiltration exexperiments-BEST. Soil Sci. Soc. Am. J., 70, 521–532.
Lassabatere, L., Spadini, L., Delolme, C., Février, L.,
Galvez-Cloutier, R., Winiarski, T., 2007. Concomitant
Zn-Cd and Pb retention in a carbonated fluvio-glacial
deposit under both static and dynamic conditions.
Chemosphere,69,
1499–1508.
http://dx.doi.org/10.1016/j.chemosphere.2007.04.053
Lassabatere, L., Angulo-Jaramillo, R., Goutaland, D., Letellier,L., Gaudet, J.P., Winiarski, T., Delolme, C., 2010. Effect of
the settlement of sediments on water infiltration in two urban
infiltration
basins.
Geoderma,156,
316–325.
http://dx.doi.org/10.1016/j.geoderma.2010.02.031
Lassabatere, L., Yilmaz, D., Peyrard, X., Peyneau, P.E., Lenoir,T., Šimůnek, J., Angulo-Jaramillo, R., 2014. New analyticalmodel for cumulative infiltration into dual-permeability
soils. Vadose Zone J., 13. DOI: 10.2136/vzj2013.10.0181
Mualem, Y., 1976. A new model for predicting the hydraulicconductivity of unsaturated porous media. Water Resour.
Res., 12, 513–522.
Šimůnek, J., van Genuchten, M.T., 2008. Modeling
nonequilibrium flow and transport processes using
HYDRUS. Vadose Zone J., 7, 782–797.
Šimůnek, J., Jarvis, N.J., van Genuchten, M.T., Gärdenäs, A.,
2003. Review and comparison of models for describing nonequilibrium and preferential flow and transport in the vadosezone. J. Hydrol., 272, 14–35.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2008.
Development and applications of the HYDRUS andSTANMOD software packages and related codes. Vadose
Zone J., 7, 587–600.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent
developments and applications of the HYDRUS computer
software packages. Vadose Zone J., 15, 7, pp. 25, DOI:
10.2136/vzj2016.04.0033
van Genuchten, M.T., 1980. A closed-form equation forpredicting the hydraulic conductivity of unsaturated soils.
Soil Sci. Soc. Am. J., 44, 892–898.
Vogel, T., Gerke, H.H., Zhang, R., van Genuchten, M.T., 2000.
Modeling flow and transport in a two-dimensional dualpermeability system with spatially variable hydraulicproperties. J. Hydrol., 238, 78–89.
Walter, M.T., Kim, J.-S., Steenhuis, T.S., Parlange, J.-Y.,
Heilig, A., Braddock, R.D., Selker, J.S., Boll, J., 2000.
Funneled flow mechanisms in a sloping layered soil:
Laboratory investigation. Water Resour. Res., 36, 841–849.
Winiarski, T., Bedell, J.-P., Delolme, C., Perrodin, Y., 2006.
The impact of stormwater on a soil profile in an infiltration
basin. Hydrogeol. J., 14, 1244–1251. DOI: 10.1007/s10040006-0073-9.
Winiarski, T., Lassabatere, L., Angulo-Jaramillo, R.,
Goutaland, D., 2013. Characterization of the heterogeneous
flow and pollutant transfer in the unsaturated zone in thefluvio-glacial deposit. Procedia Environ. Sci., 19, 955–964.
Received 9 August 2016
Accepted 9 November 2016

UnauthenticatedDownload Date | 7/25/17 6:13 AM

221