Journal of Hydrology 598 (2021) 126301

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Unsaturated flow effects on solute transport in porous media
Luwen Zhuang a, b, c, *, Amir Raoof c, Mojtaba G. Mahmoodlu c, d, Sara Biekart c, Riemer de Witte c,
Lubna Badi c, Martinus Th. van Genuchten c, e, Kairong Lin a, b
aSchool of Civil Engineering, Sun Yat-sen University, Zhuhai 519082, China
Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai), Zhuhai 519082, Chinac
Department of Earth Sciences, Utrecht University, P.O. Box 80021, 3508 TA Utrecht, NetherlandsdDepartment of Watershed and Rangeland Management, Gonbad Kavous University, Irane
Center for Environmental Studies, São Paulo State University, UNESP, Rio Claro, SP 13506-900, Brazilb

A R T I C L E I N F O

A B S T R A C T

This manuscript was handled by Jiri Simunek,
Editor-in-Chief

A major contaminant transport process in soils is hydrodynamic dispersion by affecting the spreading and arrivalof surface-applied pollutants at underlying groundwater reservoirs. When a soil is unsaturated, hydrodynamic
dispersion is very much affected by soil water saturation. Centimeter- and decimeter-scale column experiments
were carried out to explore the effects of fluid saturation and particle size on the unsaturated solute dispersivity.
Measured in-situ breakthrough curves were analyzed in terms of both classical advection–dispersion and dualporosity (mobile-immobile) type transport equations. A clear non-monotonic relationship was found between
the dispersivity and soil water saturation. The extent of non-monotonicity was more pronounced for a relatively
coarse-textured sand compared to a finer sand. This finding has been reported rarely before; it explains some of
the inconsistencies of saturation-dispersivity relationships in the literature.

Keywords:
Dispersion
Dispersivity
Solute transport
Unsaturated zone

1. Introduction
Contaminants released at or near the earth surface may travel
through the soil unsaturated zone to arrive at and pollute underlying
groundwater resources. Knowledge of the processes governing the fate
and transport of pollutants in the unsaturated zone is a major aspect ofrisk assessments and remediation of contaminated aquifers (Bear andAlexander, 2010). Many studies have been carried out to investigate
solute transport processes in the saturated zone (e.g., Fitch and Jia,
1996; Chiogna et al., 2010; Rolle et al., 2010; Gai et al., 2011). Soil
texture and travel distance are well-known for influencing the dispersion
process significantly. An additional complexity in the unsaturated(vadose) zone is the effect of soil water content (or fluid saturation) on
solute dispersion (Bear, 1988), an issue that remains relatively illdefined (Scheidegger, 1961; Yule and Gardner, 1978; Bolt, 1979; De
Smedt et al., 1986; Gelhar, 1986). Several studies (e.g., Kirda et al.,
1973; De Smedt et al., 1986; Maraqa et al., 1997; Matsubayashi et al.,
1997; Devkota et al., 1998; Padilla et al., 1999; Kumahor et al., 2015)
have shown an increase in the solute dispersivity for unsaturated soils
compared to fully saturated conditions. Still other studies (Yule andGardner, 1978; Costa and Prunty, 2006) found that dispersion remains

constant, regardless of the saturation level. A few recent laboratory-scale
experimental studies (Bunsri et al., 2008; Toride et al., 2003; Kar­adimitriou et al., 2016) have shown a non-monotonic relationship be­
tween the solute dispersivity and soil water content, with thedispersivity increasing to a maximum value (referred to here as the
maximum dispersivity) at some intermediate saturation and thendecreasing with further desaturation. Raoof and Hassanizadeh (2013)
could explain this non-monotonic dispersivity behavior based on a nu­
merical study of pore-scale fluid flow and solute transport processes.
The laboratory studies referred to above were mostly conductedusing repacked soil columns. More uncertainties are expected for fieldscale undisturbed soils because of large-scale soil heterogeneity, soil
layering, the presence of fractures and macropores, spatial variations in
the flow velocity within the transport domain, time-dependent bound­ary conditions, root water uptake, and other complexities (e.g., Van­
derborght et al., 2000, 2001; Javaux et al., 2006; Vanderborght andVereecken, 2007; Koestel et al., 2009). Carefully constructed laboratory
tracer experiments permit a far more precise assessment of the effects ofsoil moisture content on the dispersivity, as compared to field-scalesituations. The discrepancies in the literature referred to above werelikely due in part also to different soil types being used in the various* Corresponding author at: School of Civil Engineering, Sun Yat-sen University, Zhuhai 519082, China.
E-mail address: zhuanglw3@mail.sysu.edu.cn (L. Zhuang).
https://doi.org/10.1016/j.jhydrol.2021.126301
Received 15 September 2020; Received in revised form 4 April 2021; Accepted 5 April 2021
Available online 8 April 2021
0022-1694/© 2021 Elsevier B.V. All rights reserved.

L. Zhuang et al.

Journal of Hydrology 598 (2021) 126301

studies. Soil texture and soil structure affect water and air distributions
at a given saturation, which in turn affects prevailing fluid velocity
variations and solute mixing, resulting in different solute spreading anddispersion phenomena. Additionally, scale effects are important for so­lute dispersion (e.g., Dagan, 1986; Butters et al., 1989; Bromly and Hinz,
2004; Bromly et al., 2007; Mayer et al., 2008); however, their impact
under unsaturated conditions remains unclear.
Another unsettled issue is the selection of an appropriate macro-scalemodel, which could be the classical equilibrium Advection-Dispersion
Equation (ADE), a dual-porosity type Mobile-Immobile Model (MIM),
or some other formulation. Some studies suggested the use of MIM typemodels (e.g., De Smet et al., 1981; Beven and Young, 1988; Geiger andDurnford, 2000), while others (e.g., Maraqa et al., 1997) show thesuitability of the ADE to model transport without requiring immobile
liquid zones. Still, the various transport models do need relationships for
macroscopic coefficients related to solute dispersion and possible
nonequilibrium mass transfer processes. Similar to the constitutive re­lationships for fluid flow, transport parameters are known to depend on
pore structure as well as on actual air and water phase distributions. A
lack of available constitutive data is often cited as a primary barrier toobtaining acceptable predictions (Toride et al., 2003).
Although advection and dispersion are generally considered to be themost important transport processes, there is ample evidence that solute
diffusion into immobile or dead-end zones also affect solute residence
times and resulting concentration distributions. Criteria are needed forchoosing a particular macro-scale model that can be used for different
soil types and fluid saturation scenarios, including variably-saturated
conditions. Unfortunately, only a few well-controlled laboratory-scalehydrodynamic dispersion experiments exist for unsaturated flow con­ditions, especially experiments that consider a broad range of satura­
tions. One major reason for this relates difficulties to establish uniform
flow conditions at relatively low fluid saturations.
In this paper we present results of a systematic numerical andexperimental study on solute transport under saturated and unsaturated
conditions. Our aim is to investigate several factors influencing unsat­
urated transport, such as soil texture and the effects of scale. We per­formed a series of solute transport experiments for different sandy soils
having a wide range of fluid saturations. Moreover, we employed two
experimental setups with different column lengths to explore the scaledependency of unsaturated solute dispersion. Breakthrough curves(BTCs) were measured in situ at different points along the flow path. The
equilibrium ADE and nonequilibrium-based MIM models were used toanalyze the experimental data and to obtain solute dispersivities atdifferent fluid saturations. The resulting dependency of the solute dis­persivity on saturation is shown and discussed.

two sands are listed in Table 1. The HYPROP evaporation device (UMS
AG, Germany) was used to obtain soil water retention data of the twosands, and to find the corresponding van Genuchten (1980) soil hy­
draulic parameters. The resulting curves are shown in Fig. 1, while theoptimized hydraulic parameter values are listed in Table 1. Comparedwith sand S1, the coarser sand S2 had a lower air entry pressure as re­flected by a larger α value, and a larger value of residual water
saturation.
2.2. Experimental apparatus and procedures
We employed a 37-cm long plexiglass column with an inner diameterof 10 cm to perform the unsaturated flow and solute transport experi­ments. The setup is shown in Fig. 2. Since the dispersivity is sensitive topossible nonuniformities of flow across the inlet boundary, we placed apaper filter on top of the column to uniformly distribute water over theentire inlet surface area. A 5-mm thick hydrophilic filter was used at thebottom of the column to serve as a capillary barrier preventing air
penetration into the sample. A vacuum pump furthermore was used at
the outlet to precisely control the outlet water pressure head.
The columns were packed by pouring dry sand into the water-filledcolumns. During packing, we regularly tapped the sands and scratched
their surface to avoid layering. To ensure complete saturation, CO2 wasinjected overnight before saturating the columns. Deionized anddegassed water was used in all experiments. Several sensors were used
for the required measurements. Three micro-tensiometers (Rhizo In­
struments, Netherlands) were installed at depths of 10.5, 18.0, and 25.5cm, respectively, to measure pressure heads along the samples. The
tensiometers consisted of a ceramic cup having dimensions of 1 cm inlength and 4 mm in diameter, as well as a small pressure transducer. We
further inserted probes (5-TE sensors, Decagon Devices, USA) at thesame depths to obtain the electrical conductivity, dielectric permittivity,
and temperature simultaneously at those locations. This allowed us toalso estimate water saturation of a given sample from the dielectric
permittivity using the Topp empirical formula (Topp et al., 1980). Wheninjecting a tracer solution, the solute concentration is linearly related tothe electrical conductivity at a given value of water saturation (Toride
et al., 2003). All data were collected using a CR1000 data logger(Campbell Scientific, UK).
Solute transport experiments were carried out under both saturated
and unsaturated flow conditions. The saturated flow experiments
involved two different flow rates. After establishing steady-state flowusing solute-free water, a three-way valve was used to start injecting a
CaCl2 solution (0.08 mol/L) for a certain pore volume, and thenreturning to solute-free water to create a solute pulse. Concentration
breakthrough curves (BTCs) were measured at three depths using the 5-

2. Materials and methods
2.1. Materials
Two different sands were used in the experiments. A relatively finesand, S1, was obtained from a mining site in Belgium, and a coarse sand,S2, from a riverbed in the Netherlands. Before their use the sands were
washed using deionized water to remove clay particles. Properties of theTable 1
Properties of sands S1 and S2 used in the experiments.
Properties

S1

S2

Particle size (mm)Mean particle diameter, d50 (mm)Uniformity coefficient, d60/d10
van Genuchten parameter, α (cm− 1)
van Genuchten parameter, n (-)
Residual water saturation, Sr (-)
Saturated conductivity, Ks (cm min− 1)
Average porosity, φ (-)

0.1–0.5
0.20
1.7
0.022
10.0
0.17
1.02
0.39

0.3–0.8
0.50
1.2
0.062
10.0
0.29
4.80
0.37

Fig. 1. Water retention data for fine sand S1 and coarse sand S2.
2

L. Zhuang et al.

Journal of Hydrology 598 (2021) 126301

et al. (2017). For the smaller sample we used gamma ray transmission tomeasure water saturations and solute concentrations, thus avoiding theinsertion of any physical sensors into the small sample. This avoided any
disturbance of the soil, including possible gaps between the surface ofthe sensors and the soil particles. Measurement of the attenuation of
gamma photons has been used often to determine the soil bulk density
and water content, but not solute concentrations. For this study, wecalibrated the attenuation to be able to also measure solute concentra­
tions. Different solution concentrations lead to different attenuation
intensities. The required relationship was based on calibration prior to
the solute transport experiments. Detailed information about the gammatransmission method is given in the Supporting Information Text S1.
The transport experiments in the S1-sandbox under saturated flow
conditions were carried out using two different water flow rates. For
unsaturated flow conditions, a continuous flux was applied to the sandsample in the sandbox starting from primary drainage. After packing thesaturated sand samples, we applied a flow rate slightly smaller than themeasured saturated hydraulic conductivity (obtained with the constanthead method) to the top of the sample. The hanging water column wasthen kept at the same level as the bottom of the sand sample. Readings ofthe two tensiometers were monitored continuously. When readings ofthe two tensiometers were identical, water saturation of the sample wasmeasured using the gamma ray system. The three-way valve subse­
quently was switched from solute-free water to a CaCl2 solution for acertain number of pore volumes, and then switched back to solute-freewater in order to create a solute pulse. The resulting concentration
breakthrough curves (BTCs) were measured every 30 s using the gamma
system.
In all we conducted eight experiments for sand S1 and ten experi­ments for sand S2 under variably-saturated flow conditions. The BTCs of
the longer columns were observed at three depths for each experiment.
For the small sample, we carried out eight experiments under bothsaturated and unsaturated flow conditions. The BTCs in the small col­umn were measured at 1.5 cm from the inlet. Each experiment was
carried out at least twice to validate the results.

Fig. 2. Schematic view of the 37-cm long column.

TE sensors. For the transport experiments during unsaturated flow, unitgradient flow conditions were established first before injecting the so­lute pulse at a given water saturation. This was done to obtain uniform
water contents (i.e., no gradients in the capillary pressure) in the sam­ples, which later allowed us to relate the observed solute dispersivity toa certain well-defined saturation value. To do so, starting from saturated
flow conditions, the water inflow rate at the top and water pressure headat the bottom of the column were decreased gradually in order tominimize local hysteresis effects. The top of the column was open to air
to keep the air pressure fixed at atmospheric level. During this process,we monitored the three tensiometers until they reached the samenegative unsaturated pressure head, thus ensuring that unit gradient
flow was achieved.
The solute dispersivity is known to be length-scale dependent. To
investigate this effect at the laboratory scale, we performed, in addition
to the experiments using 37-cm long sand columns, also transport ex­
periments using a much smaller sample. The smaller sample had di­mensions of 3 cm (height) by 3 cm (length) by 2 cm (width), filled with
S1 sand. A schematic view of the small sandbox system is shown in
Fig. 3. More details about the sample container can be found in Zhuang

2.3. Transport models
Since the solutes were injected uniformly across the entire inlet crosssection, and the lateral boundaries were no-flow, overall macroscopic
flow and solute transport in both setups can be considered to be one
dimensional. The BTCs hence could be simulated using one-dimensional
continuum scale modeling. Two different formulations were considered
for this purpose: The classical Advection Dispersion Equation (ADE) anda dual-porosity type Mobile-Immobile model (MIM). We used the
STANMOD software (Toride et al., 1995; Šimůnek et al., 1999) for all
forward and inverse calculations.
For non-reactive solute transport such as in this study, the ADE isgiven by

∂C
∂2 C
∂C
=D 2− v
∂t
∂x
∂x

(1)

where C represents the solute concentration, v denotes the pore watervelocity, D is the hydrodynamic longitudinal dispersion coefficient, t is
time, and x is the spatial coordinate. The dispersivity λ in the ADE
formulation is defined as λ = D/v, thus implying that diffusion in ourstudy had a negligible effect on longitudinal transport.
Nonequilibrium transport often exists during both saturated andunsaturated flow, leading to early arrival and tailing in observed BTCs(van Genuchten et al., 1977; De Smedt et al., 1986). One approach toaccount for this is to modify the ADE model to assume the presence ofstagnant water in relatively small or dead-end pore spaces. MIM models
consider total water saturation (S) to be made up of two regions: mobile
water saturation (Sm) and immobile water saturation (Sim), with soluteexchange between the two regions simulated as a first-order mass

Fig. 3. Schematic view of the 3-cm long sandbox system.
3

L. Zhuang et al.

Journal of Hydrology 598 (2021) 126301

transfer process. The MIM model for non-reactive transport can be
described as (Coats and Smith, 1964; van Genuchten et al., 1977)
Sm

∂Cm
∂Cm
∂2 Cm
∂Cm
+ Sim
= Sm Dm 2 − Sm vm
∂t
∂t
∂x
∂x

φSim

∂Cim
= ω(Cm − Cim )∂t

The MIM model showed a slightly better description of the nonsymmetric structure and tailing of the BTCs (R2 values can be found in
Tables S1-S12).

(2)

3.2. Dispersivities of the ADE and MIM models(3)

3.2.1. ADE modelFig. 6 shows estimated dispersivity values as a function of water
saturation using the ADE model. The results clearly indicate nonmonotonic relationships between the solute dispersivity, λ, and satura­tion, S. The relatively fine sand (S1) column showed an increase indispersivity as saturation decreased from 1.0 to 0.5, but then λ slightly
decreased when S became less than 0.5. This behavior was consistent for
all depths. Overall, λ showed a slight non-monotonic behavior with amaximum value (the maximum dispersivity) at an intermediate satu­
ration value of 0.5. We note that only one set of the BTCs was measured
for saturations less than 0.5 for sand S1. This because of the very low
permeability of the fine sand at low saturations, and hence the extremely
long times needed to obtain complete BTCs. Coarse sand S2 showed far
more visible non-monotonic behavior and a better defined maximum
dispersivity (Fig. 6b). These results confirm that non-monotonicity exist
in the dispersivity versus saturation, with the extent of nonmonotonicity depending upon particle size (being less pronounced for
our finer sand). Moreover, the S value corresponding to the maximum
dispersivity shifted to a relatively large value of around 0.6 for thecoarse sand. Our results hence indicate that the changes in λ with fluid
desaturation depend upon soil texture, thus explaining some of thediscrepancies about the unsaturated dispersivity reported in the litera­ture. For instance, Toride et al. (2003) and Bunsri et al. (2008) showed
non-monotonic relationships between λ and fluid saturation for
0.25–0.50 mm and 0.15–0.5 mm sand particle sizes, which are slightly
coarser than our sand S1. By comparison, Padilla et al. (1999) reported amore monotonic dispersivity-saturation relationship for their sandranging from 0.05 to 0.43 mm (finer than our S1 sand). We emphasize
here that the results of repacked laboratory column experiments may
differ markedly from those of undisturbed soils and/or field situations.
Many other uncertainties for field-scale undisturbed soils exist, such aslarge-scale heterogeneities, soil layering, fractures and macropores,
spatial moisture variations along the transport domain, near-surfacehydrophobicity, root water uptake, time-dependent soil surface
boundary conditions, and other complexities. For these reasons one
must be careful when comparing laboratory results for steady-state flow
conditions using from repacked coarse-textured soils directly with fieldscale observations.
At a given saturation level, the dispersivity changed with distance

where the subscripts m and im refer to the mobile and immobile regions,
respectively, φ is the porosity, ω is the mass transfer coefficient between
the mobile and immobile regions, and vm and Dm are the pore watervelocity and dispersion coefficient of the mobile liquid phase, respec­tively. We will use the parameter β for the ratio of mobile saturation tototal saturation (i.e., β = Sm/S), which implies that vm in Eq. (2) is related
to v through vm = v/β, while further using λm (=Dm/vm) for the dis­persivity in the MIM model.
We specified third-type inlet boundary conditions for both models,thus considering all concentration data to be volume-average variables(van Genuchten and Parker, 1984). The dispersion coefficient and thepore water velocity were the optimized parameters in the ADE model,
while for the MIM model we also needed to optimize simultaneously themass transfer coefficient ω and the mobile fluid ratio β (in addition to vm
and Dm). Optimizations were carried out in terms of the total resident
concentration CT given by(4)

CT = Cm Sm + Cim Sim
3. Results and discussion
3.1. Breakthrough curves

We first present a selection of the BTCs that were measured at
different locations and for different fluid saturations, as well theiranalysis in terms of the ADE and MIM models. Fitted values of thetransport parameters are listed in Tables S1 to S12 of the Supporting
Information. We note here that application of different macroscopic
models (ADE and MIM in our case) will lead to different parametervalues depending upon the formulation of the macroscopic equations.
Fig. 4 shows BTCs obtained under saturated flow conditions. Bothfine sand S1 and coarse sand S2 provided symmetrical BTCs at all
depths. Given the symmetry, the ADE model was sufficient to describe
the data (the black solid lines in Fig. 4). By comparison, observed BTCs
at intermediate saturation values, such as shown in Fig. 5 for S = 0.44,
were less symmetrical with some tailing at both the higher and lower
concentrations. This suggests the use of the MIM model in addition to
the ADE. Results are shown in Fig. 5 as solid and dash lines, respectively.

Fig. 4. Saturated flow conditions: Observed and simulated BTCs for fine sand S1 and coarse sand S2.
4

L. Zhuang et al.

Journal of Hydrology 598 (2021) 126301

Fig. 5. Unsaturated flow conditions: Observed and simulated BTCs, at S = 0.44, for fine sand S1 and coarse sand S2.

Fig. 6. ADE model results: dispersivity, λ, versus water saturation, S. at different depths for fine sand S1 (a), and coarse sand S2 (b).

due to its scale dependency. Our results showed an increase in the dis­persivity up to a distance of 18.0 cm, where it reached an asymptotic
value for sand S1. The accuracy of measurements close to the bottom (i.
e., at 25.5 cm) may have been affected by hysteresis in the unsaturated
hydraulic properties. The coarser sand, S2, was found to be less sensitive
to the length scale. To further explore this effect, we used a nondestructive method (to prevent sample disturbance) to collectdispersivity data for fine sand, S1.
3.2.2. MIM model
Changes in the dispersivity obtained with the MIM model (i.e., forthe mobile phase dispersivity, λm) as a function of water saturation forboth sands are shown in Fig. 7a. The transport parameters were obtainedusing the BTCs at the central location (i.e., at 18.0 cm) since those are

Fig. 7. Changes in solute transport parameters with saturation using the MIM model: variations in a) solute dispersivity of the mobile phase, λm, b) the mobile phasefraction, β = Sm/S, and c) the mass transfer equivalent time between the mobile and immobile phases, 1/ω.
5

L. Zhuang et al.

Journal of Hydrology 598 (2021) 126301

affected the least by the inlet and outlet boundaries. The relationshipbetween λm and water saturation for coarse sand S2 showed a nonmonotonic relation, similar as for the ADE model. However, the MIM
values of λm were six times smaller than the ADE λ values.
Mass exchange between the mobile and immobile liquid phases isinfluenced by the relative amount of these phases and the mass transfer
coefficient associated with this exchange. Fig. 7b presents the ratio (β)
between mobile water saturation (Sm) and total saturation (S). The
mobile water fraction was at its maximum during saturated conditions,but decreased as more air occupied the soil pores at lower saturation
values (Griffioen et al., 1998). For fine sand S1, β values decreasedslightly and then remained nearly constant at a value of 0.96 as S further
decreased. For coarse sand S2, β values decreased significantly at first,
but then stayed at around 0.7 with further desaturation. Differences in
the pore size and flow pathways between the two sands are the maincauses of these differences in β. Sand S1 is a relatively fine to medium
sand with relatively small pores and many connected flow pathways at agiven saturation, while sand S2 is coarser with much larger pores andfewer connected pathways. Sand S2 with its fewer connected flow
pathways was expected to show more immobile water.
Fig. 7c shows the mass transfer equivalent time, 1/ω, as a function ofsaturation under unsaturated flow conditions. For fine sand S1, when S> 0.6, the transfer time was very small, implying rapid solute exchange,
while the transfer time increased dramatically by two orders of magni­
tude or more when S < 0.6, which qualitatively is consisted with con­clusions by Griffioen et al. (1998). For coarse sand S2, the transfer time
was relatively small over the entire saturation range, with a slight in­crease only when S < 0.3. Mass transfer hence required far less time forcoarser sand S2 compared to finer sand S1.
We calculated the overall effective dispersivity also in terms of the
MIM parameters. This parameter can be defined as follows (Valocchi,
1985; Toride et al., 2003):
S2 Sm φvm
λMIM = λm + im 2
ωS

as such, as well by solute mixing between the mobile and immobile
phases as predicated by the fraction mobile liquid β (=Sm/S) and themass transfer coefficient ω), while in the ADE all of the mixing processes
are lumped into a single dispersivity parameter λ. The λMIM values forsand S1 were slightly larger than those of λ over the entire saturation
range, especially at the lower saturations, while the values of λMIM and λwere almost the same for sand S2.
3.3. Dispersivity at different scales
The above experiments were carried out using 37-cm long columns.
We additionally used a very short 3-cm column to perform transport
experiments for fine sand, S1. The resulting parameter values for the
short column are listed in Table S13. Selected BTCs under saturated andunsaturated flow conditions are shown in Fig. 9. The BTCs exhibitedsome scattering in the concentration values, most likely due to local flow
variations (the measurement area normal to the gamma ray flux was
only 6 mm in diameter) and subsequent dynamic changes in the air–­water interfaces during gamma ray passage. We acknowledge here that
BTCs from very short columns inherently are affected more by experi­mental imperfections associated with packing and the implementationof boundary conditions.
Fig. 10a shows a plot of normalized dispersivity values (relative tosaturated λ values) obtained using both the small 3-cm sample, as well as
the larger 37-cm columns. Shown are data for the point closest to theinlet of the longer column (i.e., at 10.5 cm), which possibly may havebeen affected still slightly by the input boundary. The inset provides aview with a smaller scale of the y axis. As saturation decreased, λ values
obtained with the small sample for fine sand S1 increased to around 0.6
cm at S = 0.5, and then decreased slightly, hence again showing a nonmonotonic relationship. By comparison, the dispersivities obtained forcoarse sand S2 were much more non-monotonic. Still, all of the samples
with their different soil textures and sample lengths confirmed a satu­
ration dependency of the solute dispersivity, with the dependency being
far more pronounced and non-monotonic for the coarser sand (Fig. 10a).
For the finer sand, for which we performed experiments with two
different sample lengths, saturation had much more effect on the dis­persivity of the longer column, but with the non-monotonicity being far
less as compared to the coarser sand. Fluid saturation at the turning
point (the critical saturation, Sc) was slightly lower for fine sand S1 (Sc =0.5) compared to the coarse sand S2 (Sc = 0.6).
Interestingly, the results in Fig. 10a show a similar trend as found inprevious studies where 3D pore-scale modeling was used to explore the(5)where λMIM is the overall effective MIM dispersivity, which accounts fordispersion in the mobile phase as well includes the effects of solute ex­change between the mobile and immobile phases, and where Sm and Sim
and are mobile and immobile saturations respectively.
The values of λMIM, λ, and λm for the two sands are shown in Fig. 8.
Values of λm were always much smaller than those of λ over the entire
saturation range. This was to be expected since dispersion in the MIM
formulation is represented using three parameters (i.e., by dispersivity

Fig. 8. Comparison between the effective overall dispersivity λMIM as calculated with Eq. (5), λ in the ADE model, and λm in the MIM model for fine sand S1 (a), andcoarse sand S2 (b).
6

L. Zhuang et al.

Journal of Hydrology 598 (2021) 126301

Fig. 9. Saturated and unsaturated flow conditions: Observed and simulated BTCs for fine sand S1 using the 3-cm sandbox.

Fig. 10. a) Comparison of normalized dispersivity values, λ, using the ADE model for both sands, as well as for different sample sizes for sand S1, and b) relationships
between the normalized dispersivity λ and saturation as reported by Raoof and Hassanizadeh (2013) based on 3D pore-scale network modeling of samples withdifferent pore sizes. Values of λ in a) for the larger columns were for the BTCs at 10.5 cm (closest to the inlet of the columns); solid lines show the fitted relationships
to the data.

saturation dependency of the dispersivity using samples with different
pore size distributions (Fig. 10b, Raoof and Hassanizadeh, 2013).
Consistent with that study, the experimental observations can be
explained by the underlying pore scale processes. For unsaturated con­ditions, λ is affected by the amount of mixing at the level of individual
pores, which becomes limited due to the presence of air in the larger
pores. Furthermore, at the level of sample size, dispersivity is affected bypore connectivities and whether in an unsaturated soil a cluster ofconnected saturated pores exist to form a flow pathway that percolates
though the entire sample. Sand S2 having larger pores as reflected also
by its capillary pressure–saturation curve showed the much higher dis­persivity values. For this sand, air during drainage enters larger poresand effectively disconnects them from the otherwise still connectedsaturated pore system that forms fast flow pathways through the soil.
Eliminating these fast pathways causes velocity variations to becomeless among different parts of the sample and, therefore, causes λ todecrease. Raoof and Hassanizadeh (2013) showed that this dispersivity
turning point (λc) occurs at higher saturation values for soils having
larger pore sizes, consistent with our observations (Fig. 10). For a soil
with smaller pore sizes (e.g., sand S1) the effect of the pathways andtheir elimination by the invading air phase become less distinct, which isalso supported by the pore-scale network studies. For such a sample, theturning point of the dispersivity occurs at lower saturations, whilechanges in λ are affected less by saturation compared with samples
having larger pore sizes.
To describe the relationship between the normalized dispersivity andsaturation explicitly, we propose an empirical formula with a Gaussian
shape as follows( )2−

λ(S) = ae

S− b
c

(6)

in which a, b, and c are fitting parameters, and λ(S) is the normalized
dispersivity. We applied Eq. (6) to the three sets of λ-S data. The
resulting curves are shown as solid lines in Fig. 10a, while the parametervalues are given in Table 2. The values of a and b were close to theTable 2
Values of the parameters a, b and c in Eq. (6) for sands S1 and S2.

7

Experiment

A

b

c

S1-sandbox
S1-column
S2-column

1.50
4.91
37.2

0.41
0.45
0.49

0.81
0.39
0.24

Journal of Hydrology 598 (2021) 126301

L. Zhuang et al.

maximum dispersivity and the critical saturation, respectively. Knowing
the maximum dispersivity and the critical saturation, it is then possible
to predict λ over the entire saturation range (especially saturation levels
larger than residual saturation). Relationships such as those given by Eq.
(6) may well become useful for field-scale modeling of solute transport
in the unsaturated zone.

Butters, G.L., Jury, W.A., Ernst, F.F., 1989. Field scale transport of bromide in an
unsaturated soil: 1. Experimental methodology and results. Water Resour. Res. 25,
1575–1581.
Chiogna, G., Eberhardt, C., Grathwohl, P., Cirpka, O.A., Rolle, M., 2010. Evidence ofcompound-dependent hydrodynamic and mechanical transverse dispersion bymultitracer laboratory experiments. Environ. Sci. Technol. 44, 688–693.
Coats, K.H., Smith, B.D., 1964. Dead-end pore volume and dispersion in Porous Media.
Soc. Pet. Eng. J. 4, 73–84.
Costa, J.L., Prunty, L., 2006. Solute transport in fine sandy loam soil under different flow
rates. Agric. Water Manag. 83, 111–118.
Dagan, G., 1986. Statistical theory of groundwater flow and transport: pore tolaboratory, laboratory to formation, and formation to regional scale. Water Resour.
Res. 22, 120S–134S.
De Smedt, F., Wauters, F., Sevilla, J., 1986. Study of tracer movement throughunsaturated sand. J. Hydrol. 85, 169–181.
Devkota, L., Matsubayashi, U., Takagi, F., 1998. A new form of dispersion coefficientmodel for the porous media. Annu. J. Hydraul. Eng. 42, 355–360.
Fitch, A., Jia, D.U., 1996. Solute transport in clay media: Effect of humic acid. Environ.
Sci. Technol. 30, 12–15.
Gai, K., Shi, B., Yan, X., Wang, D., 2011. Effect of dispersion on adsorption of atrazine byaqueous suspensions of fullerenes. Environ. Sci. Technol. 45, 5959–5965.
Geiger, S.L., Durnford, D.S., 2000. Infiltration in homogeneous sands and a mechanistic
model of unstable flow. Soil Sci. Soc. Am. J. 64, 460–469.
Gelhar, L.W., 1986. Stochastic subsurface hydrology from theory to applications. Water
Resour. Res. 22, 135S–145S.
Griffioen, J.W., Barry, D.A., Parlange, J.Y., 1998. Interpretation of two-region model
parameters. Water Resour. Res. 34 (3), 373–384.
Javaux, M., Vanderborght, J., Kasteel, R., Vanclooster, M., 2006. Three-dimensionalmodeling of the scale- and flow rate-dependency of dispersion in a heterogeneous
unsaturated sandy monolith. Vadose Zone J. 5 (2), 515–528.
Karadimitriou, N.K., Joekar-Niasar, V., Babaei, M., Shore, C.A., 2016. Critical role of theimmobile zone in Non-Fickian two-phase transport: a New Paradigm. Environ. Sci.
Technol. 50, 4384–4392.
Kirda, C., Nielsen, D.R., Biggar, J.W., 1973. Simultaneous transport of chloride and water
during infiltration. Soil Sci. Soc. Am. Proc. 37, 339–345.
Koestel, J., Vanderborght, J., Javaux, M., Kemna, A., Binley, A., Vereecken, H., 2009.
Noninvasive 3-D transport characterization in a sandy soil using ERT: 1.
investigating the validity of ERT-derived transport parameters. Vadose Zone J. 8 (3),711–722.
Kumahor, S.K., de Rooij, G.H., Schlüter, S., Vogel, H.-J., 2015. Water flow and solutetransport in unsaturated sand—A comprehensive experimental approach. Vadose
Zone J. 14.
Maraqa, M.A., Wallace, R.B., Voice, T.C., 1997. Effects of degree of water saturation ondispersivity and immobile water in sandy soil columns. J. Contam. Hydrol. 25,
199–218.
Matsubayashi, U., Devkota, L.P., Takagi, F., 1997. Characteristics of the dispersion
coefficient in miscible displacement through a glass beads medium. J. Hydrol. 192,
51–64.
Mayer, A., Sandman, T., Breidenbach, M., 2008. Effect of flow regime on physicalnonequilibrium transport in unsaturated porous media. Vadose Zone J. 7, 981–991.
Padilla, I.Y., Yeh, T.-C.-J., Conklin, M.H., 1999. The effect of water content on solutetransport in unsaturated porous media. Water Resour. Res. 35, 3303–3313.
Raoof, A., Hassanizadeh, S.M., 2013. Saturation-dependent solute dispersivity in porous
media: pore-scale processes. Water Resour. Res. 49, 1943–1951.
Rolle, M., Chiogna, G., Bauer, R., Griebler, C., Grathwohl, P., 2010. Isotopic fractionation
by transverse dispersion: flow-through microcosms and reactive transport modelingstudy. Environ. Sci. Technol. 44, 6167–6173.
Scheidegger, A.E., 1961. General theory of dispersion in porous media. J. Geophys. Res.
66, 3273–3278.
Simunek, J., van Genuchten, M.T., Sejna, M., Toride, N., Leij, F.J., 1999. The STANMOD
computer software for evaluating solute transport in porous media using analyticalsolutions of convection-dispersion equation. IGWMC-TPS-71, International Ground
Water Modeling Center, Colorado School of Mines, Golden, Colorado, 32p.
Smet, F. de, Wierenga, P.J., Beken, A., 1981. Theoretical and experimental study ofsolute movement through porous media with mobile and immobile water. VUB
Hydrol. 6, V.U.B., Brussel, Belgium.
Topp, G.C., Davis, J.L., Annan, A.P., 1980. Electromagnetic determination of soil watercontent: Measruements in coaxial transmission lines. Water Resour. Res. 16,
574–582.
Toride, N., Inoue, M., Leij, F.J., 2003. Hydrodynamic dispersion in an unsaturated dunesand. Soil Sci. Soc. Am. J. 67, 703–712.
Toride, N., Leij, F.J., van Genuchten, M.Th., 1995. The CXTFIT code for estimatingtransport parameters from laboratory or field tracer experiments. U.S. Salinity
Laboratory Res. Rep. 137, U.S. Salinity, Riverside, CA.
Valocchi, A.J., 1985. Validity of the local equilibrium assumption for modeling sorbing
solute transport through homogeneous soils. Water Resour. Res. 21 (6), 808–820.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
van Genuchten, M.T., Parker, J.C., 1984. Boundary conditions for displacement
experiments through short laboratory soil columns. Soil Sci. Soc. Am. J. 48,
703–708.
van Genuchten, M.T., Wierenga, P.J., O’Connor, G.A., 1977. Mass transfer studies insorbing porous media: III. Experimental evaluation with 2,4,5-T. Soil Sci. Soc. Am. J.
41, 278–285.

4. Conclusions
This study provides a comprehensive analysis of solute transport for
two different sandy soils under saturated and a wide range of unsatu­
rated flow conditions. Results showed a well-defined dependence of thesolute dispersivity on soil water saturation, with the resulting relation­
ship being non-monotonic, particularly for the coarser medium. The
established relationship between solute dispersivity and water satura­tion may prove to be useful for field-scale contaminant transport modelsof the unsaturated zone to provide a better basis for subsurface envi­
ronmental management and risk analyses. Our study involved carefully
packed homogeneous fine and coarse sands. Further studies are clearlyneeded for different natural soils, especially structured or macroporous
field soils, to test and/or improve the suggested relationship between thesolute dispersivity and water saturation in the unsaturated zone.
CRediT authorship contribution statement
Luwen Zhuang: Methodology, Software, Validation, Writing original draft. Amir Raoof: Conceptualization, Supervision, Writing review & editing, Funding acquisition. Mojtaba G. Mahmoodlu:
Methodology, Supervision. Sara Biekart: Investigation. Riemer deWitte: Investigation. Lubna Badi: Investigation. Martinus Th. van
Genuchten: Supervision, Writing - review & editing. Kairong Lin:
Writing - review & editing.
Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Acknowledgments
The authors acknowledge financial support by the Veni Talent
Scheme Award (No. 016.151.047), which is (partly) financed by the
Netherlands Organization for Scientific Research (NWO), and the Na­tional Natural Science Foundation of China (Grant No. 51822908 and
42007165). The authors would like to thank the Editor and two anon­
ymous reviewers for providing helpful comments to the manuscript.
Appendix A. Supplementary data
Supplementary data to this article can be found online at https://doi.
org/10.1016/j.jhydrol.2021.126301.
References
Bear, J., 1988. Dynamics of Fluids in Porous Media. Dover, New York.
Bear, J., Alexander, H.D.C., 2010. Modeling Groundwater Flow and Contaminant
Transport. Springer, Dordrecht.
Beven, K.J., Young, P.C., 1988. An aggregrated mixing zone model of solute transport
through porous media. J. Contam. Hydrol. 3, 129–143.
Bolt, G.H., 1979. Movement of solutes in soil: principles of adsorption/exchange
chromatography. Dev. Soil 5, 285–348.
Bromly, M., Hinz, C., 2004. Non-Fickian transport in homogeneous unsaturated repacked
sand. Water Resour. Res. 40, W07402.
Bromly, M., Hinz, C., Aylmore, L.A.G., 2007. Relation of dispersivity to properties ofhomogeneous saturated repacked soil columns. Eur. J. Soil Sci. 58, 293–301.
Bunsri, T., Sivakumar, M., Hagare, D., 2008. Influence of dispersion on transport oftracer through unsaturated porous media. J. Appl. Fluid Mech. 1, 37–44.

8

L. Zhuang et al.

Journal of Hydrology 598 (2021) 126301

Vanderborght, J., Timmerman, A., Feyen, J., 2000. Solute transport for steady-state andtransient flow in soils with and without macropores. Soil Sci. Soc. Am. J. 64 (4),1305–1317.
Vanderborght, J., Vanclooster, M., Timmerman, A., Seuntjens, P., Mallants, D., Kim, D.J.,
Jacques, D., Hubrechts, L., Gonzalez, C., Feyen, J., Diels, J., Deckers, J., 2001.
Overview of inert tracer experiments in key Belgian soil types: relation between
transport and soil morphological and hydraulic properties. Water Resour. Res. 37(12), 2873–2888.

Vanderborght, J., Vereecken, H., 2007. Review of dispersivities for transport modeling insoils. Vadose Zone J. 6, 29–52.
Yule, D.F., Gardner, W.R., 1978. Longitudinal and transverse dispersion coefficients inunsaturated plainfield sand. Water Resour. Res. 14, 582–588.
Zhuang, L., Bezerra Coelho, C.R., Hassanizadeh, S.M., van Genuchten, M.T., 2017.
Analysis of the hysteretic hydraulic properties of unsaturated soil. Vadose Zone J.
16, 9p.

9