Advances in Water Resources 105 (2017) 82–95

Contents lists available at ScienceDirect

Advances in Water Resources
journal homepage: www.elsevier.com/locate/advwatres

Multiscale modelling of dual-porosity porous media; a computational
pore-scale study for ﬂow and solute transport
Enno T. de Vries a, Amir Raoof a,∗, Martinus Th. van Genuchten a,bab

Department of Earth Sciences, Utrecht University, Utrecht, The Netherlands
NIDES Interdisciplinary Center for Social Development, Federal University of Rio de Janeiro, Rio de Janeiro, Brazila r t i c l e

i n f oArticle history:
Received 8 December 2016
Revised 18 April 2017
Accepted 19 April 2017
Available online 20 April 2017
Keywords:
Pore network model
Aggregated soil
Dual-porosity
Solute transport
Mobile-immobile transport modela b s t r a c t
Many environmental and agricultural applications involve the transport of water and dissolved constituents through aggregated soil proﬁles, or porous media that are structured, fractured or macroporous in other ways. During the past several decades, various process-based macroscopic models havebeen used to simulate contaminant transport in such media. Many of these models consider advectivedispersive transport through relatively large inter-aggregate pore domains, while exchange with thesmaller intra-aggregate pores is assumed to be controlled by diffusion. Exchange of solute between thetwo domains is often represented using a ﬁrst-order mass transfer coeﬃcient, which is commonly obtained by ﬁtting to observed data. This study aims to understand and quantify the solute exchange termby applying a dual-porosity pore-scale network model to relatively large domains, and analysing the porescale results in terms of the classical dual-porosity (mobile-immobile) transport formulation.
We examined the effects of key parameters (notably aggregate porosity and aggregate permeability)
on the main dual-porosity model parameters, i.e., the mobile water fraction (φ m ) and the mass transfer
coeﬃcient (α ). Results were obtained for a wide range of aggregate porosities (between 0.082 and 0.700).
The effect of aggregate permeability was explored by varying pore throat sizes within the aggregates. Solute breakthrough curves (BTCs) obtained with the pore-scale network model at several locations along
the domain were analysed using analytical solutions of the dual-porosity model to obtain estimates of
φ m and α . An increase in aggregate porosity was found to decrease φ m and increase α , leading to considerable tailing in the BTCs. Changes in the aggregate pore throat size affected the relative ﬂow velocity
between the intra- and inter-aggregate domains. Higher ﬂow velocities within the aggregates caused achange in the transport regime from diffusion dominated to more advection dominated. This change increased the exchange rate of solutes between the mobile and immobile domains, with a related increasein the value of the mass transfer coeﬃcient and less tailing in the BTCs.
© 2017 The Authors. Published by Elsevier Ltd.
This is an open access article under the CC BY license. (http://creativecommons.org/licenses/by/4.0/)

1. Introduction
Soil and groundwater pollution by a broad range of industrial and agricultural contaminants is an ever-increasing problem
worldwide. One issue exacerbating effective management of thesubsurface is the problem of preferential ﬂow of surface applied chemicals such as fertilizers, pesticides, trace elements andpathogenic microorganisms. Much evidence exists that preferential
ﬂow through especially the vadose zone is contributing to surface and subsurface pollution problems (e.g., Flury et al., 1994;
Abbaspour et al., 2001; Hendrickx and Flury, 2001; Allaire et al.,2009, Vogel et al., 2010; Zhang et al., 2013; Mahmoodlu et al.,∗Corresponding author.
E-mail address: a.raoof@uu.nl (A. Raoof).

2013, 2014). In attempts to describe and quantify the basic processes leading to preferential ﬂow, a large number of dual-porosity,dual-permeability, multi-porosity and multi-permeability models
and related approaches have been developed and veriﬁed. Comprehensive reviews of alternative modelling approaches are provided
by NRC (2001), Šimůnek et al. (2003), Gerke (2006), Jarvis (2007),Šimunek and van Genuchten (2008) and Köhne et al. (2009).
A commonality of many preferential ﬂow models is the assumption that the liquid phase can be divided into relatively
mobile and immobile parts representing the macropore (orinter-aggregate or fracture) domain, and the micropore (or intraaggregate or soil matrix) domain, respectively, together withappropriate coupling terms to account for the exchange of water
and/or dissolved constituents between the two domains. The
simplest formulation arises when no ﬂow occurs in the micropore
domain, and a ﬁrst-order macropore/micropore solute exchange

http://dx.doi.org/10.1016/j.advwatres.2017.04.013
0309-1708/© 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. (http://creativecommons.org/licenses/by/4.0/)

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

term is used, to lead to mobile-immobile (MIM) type dual-porosity
models of the form (Coats and Smith, 1964; van Genuchten and
Wierenga, 1976):∂ 2 c̄m∂ c̄m∂ c̄m= φm θ D m−q− α (c̄m − c̄im )
∂t
∂x
∂ x2
∂ c̄
(1 − φm )θ im = α (c̄m − c̄im )
∂t

φm θ

(1)
(2)

where the subscripts m and im refer to the mobile and immobile
region, c̄ is the average concentration, θ is the volumetric watercontent, φ m is the mobile water fraction deﬁned as θ m /θ , Dm isthe dispersion coeﬃcient, q is the volumetric ﬂuid ﬂux, and α is aﬁrst-order mass transfer coeﬃcient.
The dual-porosity model given by Eqs. (1) and (2) involves several assumptions, the most important being that advective transport in the smaller intra-aggregate pores can be neglected. Thisimplies that the overall pore water velocity distribution within the
porous medium is approximated by the step function, with onepart of the medium having an average pore water velocity equal
to vm = q/θm , while water in the other part is completely stagnant.
This assumption is generally not met since the intra-aggregate (soilmatrix) region often has some non-zero permeability, even if small
compared to the inter-aggregate (or fracture) region.
Another simplifying assumption of Eqs. (1) and (2) is that soluteexchange between the mobile and immobile regions can be described using a quasi-empirical ﬁrst-order rate term proportional
to the difference between the average concentrations of the mobile and immobile regions. Various attempts have been made toobtain a more physical basis of the mass transfer coeﬃcient, α ,in terms of such parameters as the diffusion coeﬃcient and theshape and size of the aggregates or soil matrix. This has led toa number of analytical models that explicitly considered diffusion from the inter-aggregate region into immobile intra-aggregate
regions of various shapes (e.g., Rasmuson and Neretnieks, 1980;
Tang et al., 1981; Sudicky and Frind, 1982; van Genuchten et al.,
1984; van Genuchten, 1985). These and related studies also allowed derivation of approximate relationship for the mass transfer coeﬃcient, including through the use of Laplace transforms or
moment analyses, as exempliﬁed in studies by van Genuchten and
Dalton (1986), Bolt (1979), Barker (1985), van Genuchten and Dalton (1985), Parker and Valocchi (1986), Goltz and Roberts (1987),
and Hantush and Marino (1988).
Another approach for estimating the parameters φ m and α inEqs. (1) and (2) is by direct measurement (e.g., Clothier et al., 1992;
Jaynes and Shao, 1999), or by analyzing a large number of previously published data such as shown by Maraqa (2001). A latterstudy revealed an approximately linear relationship between themass transfer coeﬃcient and the residence time of the solute in
the transport domain (Pontedeiro et al., 2010).
Both of the above assumptions (i.e., negligible advection within
the aggregates, and the use of a ﬁrst-order exchange term) require
more research. This includes how best to account for aggregate
shape and size, which are known to vary widely and may involve
various spherical, blocky, columnar, and prismatic geometries,
or mixtures thereof (e.g., Tisdall and Oades, 1982; Hillel, 2003).
Imaging techniques can be of great value for determining the sizeand shapes of aggregates, and hence can reduce the dependency
on using idealized shapes (Bultreys et al., 2015, 2016). The masstransfer coeﬃcient is, in actuality, a more complex integrated
parameter whose value depends on many porous media characteristics including pore and aggregate geometry, the solute diffusion
coeﬃcient, the intra-aggregate permeability, the relative magnitude of the mobile and immobile region, as well as the dynamics
of the overall system such as the concentration gradients and theapplied ﬂow velocity. At relatively low pore water velocities, thetime scale of solute diffusion into aggregates may be comparable

83

with the transport time scale within the macropore domain, in
which case a larger fraction of the solutes may diffuse into or outof the aggregates, thus limiting any tailing in observed BTCs. At
the other extreme, at relatively high pore water velocities, macroand micro-pores may become essentially disconnected, leading
to negligible solute exchange and possibility of dual peaks in the
BTCs (Zhou et al., 2014).
The value of the mass transfer coeﬃcient, α , is commonly
obtained by ﬁtting macroscopic models to observed solute BTCs
(Toride et al., 1995). Since BTCs often show tailing, this approach
may be very time-consuming in terms of getting appropriate resolution in the data. Moreover, the BTC data then provide information only of the macroscopic concentrations at selected observation
points within the medium, or from column outﬂow experiments,
thus providing little insight into the internal concentration distributions and interactions between the mobile and immobile
zones. This makes it diﬃcult to estimate contribution of different
transport processes into the BTCs and extrapolate such BTC data
to other transport regimes. Methods are hence needed to provide
information in a systematic manner about the internal state ofaggregated media, including quantiﬁcation of the basic transport
processes operating at the microscopic level. An alternative would
be to obtain the BTCs using pore scale modelling (Raoof et al.,2010; Raoof and Hassanizadeh, 2013). For transport in unsaturated
non-aggregated media, Raoof and Hassanizadeh (2013) used pore
network modelling to obtain several BTCs in this manner, which
would be diﬃcult and time consuming to obtain experimentally.
Advanced 3D X-ray microtomography and related imaging techniques are now being increasingly used to obtain non-destructive
visualizations of pore structures (Allaire et al., 2009; Dal Ferro etal., 2013; Mangalassery et al., 2013; Zhou et al., 2013; Martínez etal., 2015). This includes studies of dynamic processes such as ﬂuid
ﬂow and structural dynamic processes (Cnudde and Boone, 2013;
Wildenschild and Sheppard, 2013). Soil aggregates for such studies
could be imaged and analysed individually to collect data on poremorphology and connectivity (Zhou et al., 2013; Dal Ferro et al.,2013), or direct ﬂuid ﬂow experiments could be performed on adual-porosity medium. Resulting information can then be utilized
to construct pore network structures needed for pore-scale ﬂow
and/or transport models.
Various pore-scale modelling approaches are now being pursued, with differences pertaining to the speciﬁc mathematical formulation such as the use of direct numerical solutions (Bijeljic etal., 2013; Fathi et al., 2017a, 2017b), Lattice Boltzmann methods(Jafari et al., 2011) or pore network modelling (Raoof et al., 2013),
including differences in the invoked modelling resolution and required computational time (Sahimi, 2011; Bultreys et al., 2016).
Multi-scale pore network modelling have been applied for simulating two-phase ﬂow (Bultreys et al., 2015; Jiang et al., 2013;
Prodanović et al., 2015) as well as solute transport (Bijeljic etal., 2013). and also using grain ﬁlling and pore ﬁlling method(Mehmani and Prodanović, 2014).
The aim of our study was to use a pore-scale network model tosimulate ﬂow and transport in dual-porosity domains containing alarge number of aggregates. The pore network modelling approachassumes that the porous medium continuum can be divided intopore elements made up of pore bodies representing the larger
voids in the medium, and pore throats representing narrow openings connecting the pore bodies (Raoof et al., 2013; Bultreys et al.,2016). Applying mass balance equations, ﬂow and transport is simulated within each individual pore. Averaging over a large numberof pores will then allow estimation of the macroscopic transport
properties for porous media containing unimodal pore sizes (Raoofet al., 2010, 2013), or for multi-scale media with multimodal poresize distributions (Bultreys et al., 2015; Mehmani and Prodanović,
2014; Mehmani et al., 2015). The obtained information at the84

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

pore-scale makes it possible to relate macroscopic transport properties to the underlying physical pore-scale processes and pore size
distributions, and their connectivities (Raoof et al., 2010; Raoofand Hassanizadeh, 2013). Using such multi-scale porous media
requires description of the connectivity between the macro- andmicro-porosity domains, which computationally can be particularly
challenging for heterogeneous multimodal pore scale domains.
Our objective is to develop a dual-porosity pore structure to
represent an aggregated medium, and to use this structure to simulate ﬂow and transport within the composite medium in orderto obtain pore-scale distributions of the solute concentration. The
resulting solute breakthrough curves are then compared with solutions of the macroscopic equations given by Eqs (1) and (2) to obtain estimates of the fraction of mobile water (φ m ), and the masstransfer coeﬃcient (α ). A large number of simulations will be carried out to obtain dependencies of φ m and α on such soil aggregate properties as porosity and internal permeability. Resulting insight may be useful for ﬁeld scale models to improve predictions
needed at the larger scales for evaluating alternative management
or contaminant remediation strategies, as well as industrial porous
media applications.
2. Theoretical development
In this section we describe the pore-scale network model thatwas used to simulate ﬂow and transport in the aggregated dualporosity pore structure. After a brief review of the pore-scale network modelling approach in general, details are provided on how
the pore-scale model was used to generate pore structures andconduct simulations for an aggregated medium.
2.1. Pore network model

certain size was assumed to contain a large number of micropores. Algorithms were developed to connect the collection of pores
located at the outer surface of the aggregates to their adjacentmacropores. We generated for this purpose several domains withdifferent aggregate fractions, with the aggregates themselves randomly placed within the medium. An example is shown in Fig. 1.
Once the pore structures were constructed, ﬂow and transport processes were simulated within the pore network domain using the
PoreFlow package developed by Raoof et al., (2013), which is capable of simulating saturated and variably-saturated ﬂow, as wellas multi-component reactive transport, within arbitrary pore structures.
Since the adopted pore-scale model has been explained in detail elsewhere (Raoof and Hassanizadeh, 2010, 2012; Raoof et al.,2010, 2013), only a brief overview is provided here of the mainequations governing ﬂow and transport within the porous medium.
Assuming laminar ﬂow, the volumetric ﬂow rate within each pore
was obtained using the Hagen-Poiseulle equation:

qi j = gi j p j − pi(3)where qij is the volumetric volume rate through the pore throat
between two adjacent connected pore bodies i and j, pi and pj arepressures of the two adjacent pore bodies, and gij is the conductance of the pore throat which, for a pore throat with a cylindricalcross section, is given by:

gi j =

π R4i j
8 μl i j(4)where Rij is the pore throat radius, lij is the throat length, and μ is
the ﬂuid dynamic viscosity.
Considering incompressible ﬂow, the continuity equation for
ﬂow is applied at pore junctions:zi

Pore network modelling considers a porous medium as a system of pore elements composed of pore bodies and pore throats.
While pore bodies, due to their larger size, mainly control theporosity of the medium, the smaller pore throats determine thehydraulic conductivity. The complete pore structure of an aggregated soil requires the use of two types of pore spaces with different sizes. The aggregated medium in our study was obtainedby ﬁrst generating a macropore domain. Since pore connectivity
is an important topological property of the medium, we used apore network model with random pore connectivity (Raoof & Hassanizadeh, 2012) to represent the randomness of real porous media. Once the macropore domain was generated, several aggregateswere superimposed within the overall domain. Each aggregate ofqi j = 0;

j = 1, 2, . . . , zi ,

(5)

j=1

where qij is the volumetric ﬂow rate within a pore throat frompore body i to pore body j, and zi is the pore coordination numberof pore i.
Combining Eqs. (3)–(5) for all pores results in a linear systemof equations with a sparse, symmetric and positive-deﬁnite coeﬃcient matrix, which is to be solved for the pore body pressures. Considering the sample as a representative elementary volume (REV), the average pore water velocity, v̄, can be calculated as(Raoof et al., 2013):v̄ =Qtot LVf(6)

Fig. 1. Construction of the aggregated domain: macroscopic pore domain (the left plot) as well as the aggregate domain (the middle plot) are created and superimposed toprovide the ultimate multi-scale domain (the right plot).

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

85

Table 1
Selected literature on aggregate properties: aggregate diameters, pore sizes, and porosities.
Reference

Aggregate diameter (mm)Pore size (μm)

Mean aggregate porosity [–]Imaging method

Czachor et al. (2015)
Dal Ferro et al. (2012)
Dal Ferro et al. (2013)
Khan et al. (2012)
Peth et al. (2008)
Rabbi et al. (2015)
Zhou et al. (2013)

3–4
1–2
5–6
1–5
3–5
2–4
1.5–2

2.5–100
0.074–100
12.5–200
1–85
5–50
5.2–60
10–500

0.35
0.37
0.11
0.40
0.13
0.048
0.14

X-ray microtomography and paraﬃn intrusion
X-ray microtomography and mercury intrusion
X-ray microtomography
Synchrotron-based nano-CT
X-ray microtomography
X-ray microtomography
X-ray microtomography

where Qtot is the total discharge rate through the pore network,
which can be determined at the inlet or outlet of the pore networkas the sum of all ﬂuxes, L is the length of the pore network, and
Vf is the total volume of the ﬂuid phase within the pore network.
The intrinsic permeability, κ , of the sample can be determinedusing Darcy’s law:κ=μQtot LA P

(7)

where μ is the dynamic viscosity, P is the pressure difference between the inlet and outlet pores, and A is the cross-sectional areaof the pore network.
Solute transport through the pore network was described forthe general case involving both advective transport and diffusion(Vasilyev et al., 2012). Calculations were done by considering eachpore element (i.e., a pore body or a pore throat) as the control
volume. We used a backward Euler method for the temporal discretization and ﬁrst-order upwind and central schemes for spatial
discretization of the advection and diffusion terms, respectively.
For a given pore body i, one can write the following mass balance
equation:

throat



Nin
zi


ci j − cid ciVi=qi j ci j − Qi ci +D0 Ai j
,
dt
li j
j=1

(8)

j=1

where Vi is the volume of the pore body, Qi is the total volumetricrate going out of the pore body, Aij is the cross-sectional area ofthe pore throat, D0 is the ionic or molecular diffusion coeﬃcient,ci is the concentration in the pore body, cij is the concentration inthroatthe pore throat, lij is again the length of the throat, while Nin
represents the number of pore throats ﬂowing into pore body i.
Similarly, the mass balance equation for solutes in a pore throat
ij may be written as:







c j − ci j
ci − ci jd ci jVi j= qi j c j − qi j ci j + D0 Ai j
+ D0 Ai j
dt
li j
li j



(9)

which assumes that pore body j is the upstream node. Eqs. (8) and(10) were solved using a fully implicit scheme available in the
PoreFlow software package (Raoof et al., 2013). Flux-averaged
breakthrough curves at selected points were obtained by averaging
concentrations over the network cross-section at successive times
as follows: xc (x, t ) =Ntci x, t
i
Ntx
Qi
i

(

)Qi



1
c0

i = 1, 2, . . . , Nt

(10)where c(x,t) is the normalized average concentration at location x
and time t, Ntx refers to the total number of pores at location x,
and c0 is the input concentration of the solute.
2.2. Application to aggregated media
Imaging techniques such as X-ray microtomography are being
increasingly used to visualize the inner structures of aggregates
and to provide information on pore size distribution, aggregate

porosity, and aggregate size and shape (Mangalassery et al., 2013;
Martínez et al., 2015; Zhou et al., 2013; Dal Ferro et al., 2013). We
note that studies using only X-ray microtomography often show
lower aggregate porosities as compared to studies using multiple
measuring techniques (Czachor et al., 2015; Dal Ferro et al., 2012,2013). The resolution of X-ray tomography has been shown to bea limiting factor since it tends to neglect the presence of belowresolution pores. Additional techniques such as mercury intrusion
hence are often used to obtain a more realistic value for the porosity of the aggregates (Dal Ferro et al., 2013). Focused ion beamscanning electron microscopy (FIB-SEM) can provide much higher
resolution images, as compared to X-ray tomography, for analysing
microporosity (Hemes et al., 2015). Using a combination of multiple techniques applicable to different scales will be most valuable
for multi-scale pore topology applications such as in this study.
Table 1 shows aggregate sizes, pore sizes and aggregate porosities
as obtained from the literature. We used these data to establish as
realistic ranges of properties as possible for the aggregates used inour calculations.
Transport processes within the aggregated media we considered were studied by ﬁrst creating a reference model. The reference model was then modiﬁed to consider i) domains with different aggregate porosities by changing pore body sizes, and ii)
domains with different aggregate permeabilities by changing porethroat sizes to explore the effects of ﬂow velocity differences between the macro- and aggregate domains. The reference model
had a physical size of 100 × 20 × 20 mm3 and contained a total of
10 0 0 aggregates. The mean distance between pore bodies in theaggregates was taken to be 0.1 of the separation distance between
the macropores, while the mean coordination number was taken
to be 5. The minimum and maximum aggregate diameters were
set to 1.0 and 3.0 mm, respectively. Pore body sizes and aggregate
throat sizes were assigned from a truncated log-normal distribution given by (Raoof & Hassanizadeh, 2012).

     2√
2exp − 12 ln RRmi /σ
 √
 √
f ( Ri , σ ) = √π σ 2 Ri er f ln RRmax/ 2σ 2 −er f ln RRmin/ 2σ 2mm(11)
where Rmin , Rmax and Rm are the minimum, maximum and meanof the distribution, and σ 2 is the variance of the distribution. The
macropore and micropore bodies each had their own size distribution. The macropore throat sizes were set as the smallest of theadjacent pore body sizes. The throat lengths were calculated bydetermining the length between the center points of two adjacentpore bodies of the throat, while subtracting the radii of the twopore bodies.
Table 2 provides the pore size statistics used for the calculations, with Cases I-3 and II-2 serving as the reference models.
Case I calculations were meant to show the effects of aggregate
porosity, and Case II calculations the effects of aggregate permeability. Based on the data in Table 1, the single aggregate porosity (Case I in Table 2) was varied between 0.082 and 0.700 in 6
steps by changing the aggregate pore body size distributions. The

86

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

Permeability
effects

Porosity effects

Table 2
Minimum (Rmin ), maximum (Rmax ) and mean (Rm ) pore body radii (Case I) and pore throat radii (Case II), as
well as standard deviations (σ ), of the pore-size distributions used in the calculations.

3. Results and discussion
Before presenting results of the various pore network calculations, we ﬁrst summarize the pore network that was used for thereference simulations.
3.1. Reference model
The structured macropore network consisted of 12,929 porethroats connected to each other at 6850 junctions serving as the

Table 3
Statistics of the generated aggregates.
Property

Value∗

Number of aggregates
Mean aggregate diameter (mm)Aggregate pore bodies in domain
Mean pore bodies per aggregate
Aggregate pore throats in domain
Mean pore throats per aggregate
Mean coordination number

10 0 02.141 (0.475)
965,985966 (578)∗
2,391,040
2391 (1424)∗
4.95∗

Values in parenthesis indicate standard deviations.

0.4

Normalized frequency

aggregate throat radii as well as the macropore sizes were kept
constant in these simulations to examine the effects of varying aggregate porosities on ﬂow and transport. For the Case II simulations (Table 2), the aggregate permeability was varied by changing the aggregate throat sizes in 4 steps, while the aggregate bodyand macropore body/throat sizes were kept constant. Simulation II1 had smaller throats, and simulations II-3 and II-4 larger throats,
compared to the reference model (Case II-2).
The presence of macropores and aggregates with distinctly different micropore sizes created various pore systems. To explore
the contribution of the different pore systems, porosities were calculated separately for the macropores and aggregates. Pore bodies were presumed to be spherically shaped, and pore throats tobe cylindrical capillaries. Using the pore volumes of each domain,
the domain porosities could be obtained, with the total porosity ofthe dual-porosity medium being simply the sum of the macropore
and micropore domain porosities. For all simulations we kept themean pore water velocity constant at a value of 1.50 × 10−5 m/s.
The molecular diffusion coeﬃcient, D0 , was set at 1.6 × 10−9 m2 /s,
the viscosity at 0.001 Pa s, and the ﬂuid density at 1 g/cm3 .
For all simulations we injected into the pore network a pulseof one pore volume (equivalent to 66 min) of a tracer having a relative concentration of 1.0. BTCs of the average concentrations asobtained with the pore network model were analysed using the
CXTFIT program (Toride et al., 1995) within the STANMOD software of Šimunek et al. (20 0 0), leading to estimates of the solute
dispersivity, λ = Dm/ v, the fraction of mobile water, φ m , and themass transfer coeﬃcient, α , in Eqs. (1) and (2). This in turn also
produced values of the mobile water content, θ m , and the immobile water content, θ im . The intrinsic permeability was determined
for the domain as a whole, as well as for the macropore domain
separately (i.e., excluding the presence of the aggregates). In addition, simulations were performed to obtain the intrinsic permeability of the single aggregates. For the Case II simulations, basedon the simulated pore velocities, the characteristic time scales ofdiffusion and advection were calculated to explore the effects ofthe porewater velocity ratio between the macropore and aggregate
domains on the calculations.

0.1
Pore throat agg.
Pore body agg.
Pore throat macro
Pore body macro

0.3

0.08
0.06

0.2
0.04
0.1

0.02

0
0

0.05

0.1

0.15

0.2

0
0.25

Radius [mm]Fig. 2. Pore body (the right vertical axis) and pore throat (the left vertical axis)
size distributions used for the macropore and aggregate domains of the reference
model.

pore bodies, while the mean pore coordination number was 3.7. A
total of 10 0 0 aggregates were randomly placed within the macropore network. Statistics of the generated pore network are given
in Table 3. The total number of pores within the created computational domain was 3,357,025. The macropore domain porosity
was 0.198, and the mean porosity of the single aggregates 0.204.
Fig. 2 shows the generated pore body and pore throat size distributions of the reference model, while Fig. 3 shows the generated
dual-porosity pore network.
3.2. Case I: porosity effects
Based on the information in Table 1, six different pore body size
distributions were generated as displayed in Fig. 4. Table 4 summarizes the porosities of the micropore and macropore domains
for the Case I simulations. The porosity of the single aggregatesranged from very small (0.082) to very large (0.700). The aggregate
domain porosity ranged from 0.013 to 0.114, and the total domain

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

87Fig. 3. Pore radii of the constructed pore network containing 10 0 0 aggregates and composed of a total of 3,357,025 pores.

Normalized frequency

0.08

Case I-1
Case I-2
Case I-3

0.06

Case I-4
Case I-5
Case I-6

0.04

0.02

0
0.04

0.05

0.06

0.07

0.08

0.09

0.1

0.11

Radius [mm]Fig. 4. Aggregate pore body size distributions for the Case I pore network calculations. The corresponding porosity values are provided in Table 4.

porosity (micropore and macropore domains combined) from 0.214
to 0.315 (the porosity of the macropore domain was kept constant
at 0.198).
Fig. 5 shows the simulated solute concentration distributions
for Cases I-1, I-3, and I-6 after injecting 1.0 pore volume of solute into the medium. The distributions are given at T = 1.1 porevolume, with pore volume deﬁned as T = qt /(φ t L), where L is thelength of the medium (100 mm). Concentration distributions atother times are provided in the Supplementary Material, as wella video clip showing the effect of aggregate porosity. A pulse injection was used to show both solute diffusion into the aggregate
domain and subsequent back diffusion from the aggregate domain
into the macropore domain. The plots in Fig. 5 show increased solute spreading as the porosity of single aggregates increased. CaseI-6 with its very high single aggregate porosity (φ sa = 0.700) was

included as a mere limiting case. The other extreme would be if
no aggregate porosity were to be present (φ sa = 0), in which caseno intra-aggregate diffusion would take place and the standard
advection-dispersion equation presumably would apply.
Fig. 6 shows calculated BTCs (indicated by symbols) obtainedwith the pore network model at 30, 60 and 90 mm distancesfrom the inlet for all Case I scenarios. Also included are the ﬁtted curves obtained with the macroscopic mobile-immobile (MIM)
dual-porosity model given by Eqs. (1) and (2). We optimized forthis purpose the parameters (λm , φ m and α ) simultaneously to thethree curves (i.e., at x = 30, 60 and 90 mm) obtained for each domain with a particular porosity of single aggregates. R2 values for
all Case I calculations, as well for the various Case II examples tobe discussed later, were always higher than 0.9995, thus reﬂecting
excellent descriptions of the pore-scale modelling results using themacroscopic formulation. Fig. 7 shows the ﬁtted BTCs at 90 mm forall six scenarios. Similar BTCs at 30 and 60 mm are provided in the
Supplementary Materials.
The results in Fig. 6 indicate increased spreading of the BTCs
obtained at larger distances. Also, while the BTCs of Case I-1 were
nearly symmetrical, Case I-6 with its highest porosity of single aggregates produced far more asymmetry and tailing in the distributions. Higher porosities of single aggregates (and hence more relatively immobile water) clearly enhanced the extent of preferential ﬂow as exhibited by early breakthrough and late-time tailingof the curves in Fig. 7. This agrees with previous study done on
highly heterogenous rocks by Mehmani et al. (2015) who showedthat an increasing fraction of microporosity causes more tailing of
the observed solute plume. Our ﬁndings are also consistent withseveral other previous studies (e.g., van Genuchten and Wierenga,
1976; Brusseau, 1993; Shukla, 2013).
Table 5 provides the transport parameters obtained for the six
Case I calculations. Results indicate that a higher aggregate porosity caused a decrease in the mobile water fraction, φ m , and an

Table 4Micro- and macro-porosities for the Case I calculations with variable aggregate porosities.

88

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

Fig. 5. Progress of advective transport and dispersion into the domain for Cases I-1, I-3 and I-6 after 1.1 pore volumes. A total of 1.0 pore volumes of tracer was injected.

increase in the mass transfer coeﬃcient, α . The mobile watercontent, θ m , was found to be approximately 0.21 for all Case I simulations, which is slightly higher than the macroporosity of 0.198.
For all of the simulations, the mobile water content was slightly
larger than the macro domain porosity. This is because the calculated fraction of mobile water (φ m ) at the macroscopic level is determined by the ﬂow lines affected by the presence of both macropores and aggregates. Fig. 8 shows that a clear relationship exists
between the mobile water content, φ m , and the fraction of macropores. As compared with the porosity of macropores, the fractionof mobile water content using the ﬁtting resulted in 5–10% larger
values.
Fig. 9 shows the relationship we found between the mass transfer coeﬃcient, α , and the porosity of single aggregates, φ sa . The
mass transfer coeﬃcient increased with aggregate porosity as more
solute mass could be transported into aggregates having higher
porosity values and higher permeabilities as discussed next.
To further analyse and interpret the pore scale modelling results, we calculated the permeability (K) of the macro domain andthe aggregate domain, as well as the permeability ratio between
the aggregates and the macropore region (Kagg /Kmacro ). The perme-ability of the macropore domain was found to be 2.24 × 10−4 mm2 ,with the results for Kagg and Kratio given in Table 6. When theporosity of the single aggregates increased, pores come to be closer
to each other, which resulted in higher aggregate permeability
values and larger permeability ratios. With a larger permeability,
more ﬂuid may ﬂow into and through the aggregates, leading tohigher α and lower φ m values (Fig. 10).
3.3. Case II: aggregate permeability effects
Utilizing the input parameters from Table 1, pore body size distributions were generated for the various examples. We used domains with four different pore-throat size distributions for the aggregates (Fig. 11). Flow and solute transport simulations were carried out for each pore size distribution. Note that the reference
model is Case II-2, which corresponds to Case I-3 of the ﬁrst setof simulations. Table 7 provides the porosity values for the CaseII simulations. The data show that an increase in the pore throat
radii increased the single aggregate porosity as well as the total
porosity.
The progress of one pore volume of solute moving through anddispersing into the domains of Cases II-1, II-2, and II-4 at 1.1 pore

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

Case I-1 [ϕ=0.214]

1.0

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0
0

0.5

1.0

1.5

2.0

2.5

3.0

Case I-3 [ϕ=0.234]

1.0

Concentration [-]Case I-2 [ϕ=0.227]

1.0

0.8

89

0

0.5

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

1.5

2.0

2.5

3.0

Case I-4 [ϕ=0.260]

1.0

0.8

1.0

0
0

0.5

1.0 1.5 2.0 2.5
Case I-5 [ϕ=0.284]

3.0

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0.5

1.0 1.5 2.0 2.5
Case I-6 [ϕ=0.315]

3.0

0

0.5

1.0

3.0

0
0

0.5

1.0

1.5

2.0

2.5

3.0

1.5

2.0

2.5

Pore Volume [-]
30 mm

60 mm

90 mm

Fig. 6. Aggregate porosity effects: BTCs obtained with pore network modelling (symbols) and the ﬁtted MIM transport BTCs (solid lines) for the Case I simulations.

Table 5
Optimized MIM transport parameters for the Case I simulations.

1.0

Porosity effects

Concentration [-]

0.8
0.6
0.4
0.2
0
0

0.5

1.0

1.5

2.0

2.5

3.0

Pore Volume [-]Fig. 7. BTCs obtained with the macroscopic MIM model as ﬁtted to the pore network modelling results at 90 mm for Case I simulations.

volumes are shown in Fig. 12. Similar visualizations at multiple
time levels are provided in the Supplementary Material, as wella video clip showing the effect of aggregate permeability. Flow ofsolute-free water after entry of the concentration pulse into thedomain shows that solutes moved much faster into and throughthe macropore domain for Case II-4, while the residence time ofsolutes in the aggregates was higher for Case II-1. These results indicate that the exchange rate of the solute tracer between the aggregate and macropore domains for Case II-1 was lower than forCases II-2 and II-4. A higher solute exchange rate was observed
especially for Case II-4 due to the fact that solutes in the aggregates were almost instantly leached back into the macropore domain (also see Supplementary Material).
Fig. 13 compares the BTCs for the individual simulations, while
Fig. 14 shows a comparison of the different simulations for Case
II at a longitudinal distance of 90 mm. BTCs obtained at 30 and
60 mm are provided in the Supplementary Materials. The plots

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

1

0.8

-3

R = 0.992

0.7
0.6
0.5
0.5

0.6

0.7

0.8

0.9

Fig. 8. Plot of the mobile water content (φ m ) against the fraction of macroporosity
for Case I simulations. A value of 1 for φ m would mean the absence of any aggregates.

4.5

0.9

4

0.85

3.5

0.8

3

0.75

2.5

0.7

2
1.5
0

0.02

0.04

0.06

0.08

Fig. 10. Values of the mobile water fraction (φ m ) and the mass transfer coeﬃcient(α ) plotted against the permeability ratio (Kagg /Kmacro ) for Case I.

Case II-1
Case II-2= 5.52 10-5 sa + 9.59 10-6
R = 0.996

22

Normalized frequency

0.1

3

0.1

Permeability ratio [-]

10-5

4

[min-1]

0.95

0.65

1

Fraction of macro porosity [-]

5

10-5
5

m

[-]

= 1.05x + 9.97 10

m
2

m

[-]

0.9

1

[min-1]90Case II-3
Case II-4

0.08
0.06
0.04
0.02

1
0
0

0
0

0.1

0.2

0.3

0.4
sa

0.5

0.6

0.7

[-]

0.01

0.02

0.03

0.04

0.05

Radius [mm]Fig. 11. Aggregate pore throat size distributions for Case II pore network calculations. The corresponding porosity values are provided in Table 7.

Fig. 9. Plot of the mass transfer coeﬃcient (α ) against the porosity of single aggregates (φ sa ) for Case I.
Table 6
Calculated permeability values for the aggregates, and the permeability ratio(Kagg /Kmacro ) for Case I.

show that reduced pore throat sizes lead to earlier breakthrough ofsolute, increased attenuation of the peak concentration, and morepositive skewness of the BTCs in the form of tailing (e.g., Case II-1).
Increasing pore throat sizes, on the other hand, resulted in delayed
breakthrough curves with less tailing and more symmetrical BTCs.
This is consistent with the visualizations in Fig. 12 and the Supplementary Materials, which show a more uniform solute front for
Case II-4 as compared to Case II-1. The latter case exhibits more
preferential ﬂow and more pronounced tailing.
Utilizing the BTCs, transport parameters were again obtainedfor each simulation using the CXTFIT program. The results are
given in Table 8. The obtained mobile fraction, φ m , was found tobe higher than the volume fraction of macropores for Cases II-1,

II-2, and II-3, while φ m was lower than the volume fraction ofmacropores for Cases II-4. These results are due to the hydrodynamic effects caused by the presence of aggregates in the domain.
Smaller pore throats create relatively impermeable aggregates,
which causes the overall ﬂuid ﬂow path through the macropores
to become more tortious, leading to higher immobile fractions
at the macroscopic scale. On the other hand, the presence of
aggregates with larger pore throats causes a less tortuous ﬂow
path as aggregates provide a noticeable contribution to the total
ﬂow, leading to lower immobile water fractions at the macroscopic
scale. The mass transfer rate, α , was lowest for Case II-1 and highest for Case II-4. This is in agreement with our earlier observations
about the visualizations.
Table 9 shows calculated pore water velocities of the macropore
and micropore (aggregate) domains. The data indicate that larger
pore throat sizes lead to higher pore water velocities. Fig. 15 shows
plots φ m and α against the velocity ratio, vagg/ vma . The mass transfer rate, α , decreased with a lower velocity ratio. In contrast with
Case I, an increase in α did not lead to longer tailing but produced
more symmetric BTCs. Increasing the aggregate pore water velocities caused the dominant transport process in the micropore domains to change from diffusion dominated to advection dominated.
Interactions between the aggregate and macropore domains consequently increased, leading to higher values of the mass transfer
coeﬃcient, α . We note here that the two-region mobile-immobile
model (MIM) considers the immobile domain to be completely

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

91

Table 7Micro- and macro-porosities for the Case II calculations with variable aggregate permeability.

Fig. 12. Progress of advective transport and dispersion into the domains of Cases II-1, II-2, and II-4 after 1.1 pore volumes. A total of 1.0 pore volumes of tracer was injected.

diffusion dominated, which likely contributed to the ﬂuctuationsin φ m with the velocity ratio.
When solute-free water is injected and the transport regime in
the aggregate domain is diffusion dominated, then the characteristic time for solutes to diffuse back into the macropore domain isgiven by tdif = l2 /Dm , where l is the effective length of the aggregates and Dm the molecular diffusion coeﬃcient. For an advection
dominated transport regime, on the other hand, the characteristictransport time is described by tadv = l/v̄, where l is now the domain length. These time scales are characteristic for the extent of
BTC tailing when the amount of mass stored in the aggregates isnoticeable, and when the velocity ratio is low. The advection anddiffusion time scales for the aggregate and macropore domains arelisted in Table 10. For these calculations we used the mean ﬂuid
velocity in the pore bodies and throats, while for the characteristic
length of the aggregates we used the mean aggregate radius.
We used the transport time scales to identify the dominant transport mechanism in the aggregate domain. Whentadv,agg <tdif,agg , transport in the aggregate domain is primarily advection dominated, but when tadv,agg >tdif,agg the aggregated domain
becomes more dominated by diffusion. The data in Table 10 indicate that Case II-1 was very much diffusion dominated, while theexperiments towards Case II-4 became increasingly more advection
dominated. Tailing will be observed when the characteristic time

92

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

Case II-1 [Κ=1.62×10-8 mm2]

Concentration [-]1.0

Case II-2 [Κ=1.50×10-6 mm2]

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0
0

0.5

1.0

1.5

2.0

2.5

3.0

0

Case II-3 [Κ=1.21×10-5 mm2]

1.0

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

1.0

1.5

2.0

2.5

3.0

Case II-4 [Κ=1.27×10-4 mm2]

1.0

0.8

0.5

0
0

0.5

1.0

1.5

2.0

2.5

3.0

0

0.5

1.0

1.5

2.0

2.5

3.0

Pore Volume [-]
30 mm

60 mm

90 mm

Fig. 13. Aggregate permeability effects: BTCs obtained with pore network modelling (symbols) and the ﬁtted MIM transport BTCs (solid lines) for the Case II simulations.

1.0
Case II-1 [Κ=1.62×10-8mm2]
Case II-2 [Κ=1.50×10-6 mm2]
Case II-3 [Κ=1.21×10-5 mm2]
Case II-4 [Κ=1.27×10-4 mm2]

Concentration [-]

0.8
0.6
0.4
0.2
0
0

0.5

1.0

1.5

2.0

2.5

3.0

Pore Volume [-]Fig. 14. BTCs obtained with the macroscopic MIM model as ﬁtted to the pore network modelling results at 90 mm for Case II.

for diffusive transport in the aggregates is similar or becomes lessin comparison to the advection time scale in the macropore domain.
For Case II-1, when the transport domain is leached, theamount of solute tracer remains relatively high in the micropore
domain since the time of solute tracer to be transported back into
the macropore domain is much lower than the front of the solute-

free water passing through the domain. This effect becomes lessvisible for Case II-2, and is almost absent for Case II-4. Extensive
tailing is observed when the transport regime in the aggregates isdiffusion dominated and the transport regime in the macro domain
is advection dominated. For ﬂow conditions when the velocity ratio between the two domains is close to 1, ﬂuid passes throughthe aggregates as well as through the macropore domain, leading
to relatively symmetrical BTCs. Fluid ﬂow at lower velocity ratios
is relatively much slower in the aggregates, leading to long tails
in the BTC. These ﬁndings agree with other studies, such as thoseby Brusseau (1993) and Bijeljic et al. (2013), who found that thedegree of non-equilibrium and dominance of the velocity variations between the inter-aggregate and intra-aggregate pores inﬂuence the severity of tailing.
Values of the calculated permeabilities are listed in Table 11.
With increasing pore throat sizes, the permeability increased, withthe permeability ratio increasing from a value of 7.24 × 10−5 forCase II-1 to 0.56 for Case II-4. This means that for Case II-4
the ﬂuid penetrated into aggregates almost as easy as it passed
through the macropores. For Case II-1 the ratio is very low, meaning that ﬂuid ﬂow now takes place much less through the aggregate pores as compared to the macropores. This is consistent withthe fact that considerable tailing was observed in the BTCs of Case
II-1, while Case II-4 produced more symmetrical BTCs.

Table 8Optimized MIM transport parameters for the Case II simulations.

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

93

Table 9Mean pore water velocities in the macropore and aggregated domains for Case II.

10-3
1.4

0.92

0.88

1

0.86

0.8

0.84

0.6

0.82

0.4

0.8

0.2

0.78
0

0.05

0.1

0.15

0.2

0.25

[min-1]

1.2

m

[-]

m

0.9

0
0.3

Velocity ratio [-]Fig. 15. Values of the mobile water fraction (φ m ) and the mass transfer coeﬃcient(α ) plotted against the velocity ratio (vagg/ vma ) for Case II.
Table 10
Characteristic advection and diffusion times (in minutes) for the Case II simulations.

Table 11
Calculated permeability values for the aggregates, and the permeability ratio(Kagg /Kmacro ) for Case II.

4. Conclusions
We investigated the effects of multi-scale pore sizes in a dualporosity porous medium on ﬂuid ﬂow and solute transport processes. Several dual-porosity pore networks were constructed byplacing a large number of aggregates into the macropore domain.
After construction of the aggregate domains, ﬂuid ﬂow and solutetransport were simulated using the PoreFlow (Raoof et al., 2013)
pore network model to obtain three-dimensional ﬂow ﬁelds andsolute distributions, as well as the BTCs at several locations along
the ﬂow path. The resulting BTCs were ﬁtted analysed in terms ofthe conventional macroscopic dual-porosity (MIM) transport model
to estimate several key macroscopic transport parameters. Use ofpore scale modelling provided information to link the transport
properties to the underlying pore scale processes, which would bevery diﬃcult and time-consuming to do using only experimental
methods.
A pore network model consisting of 10 0 0 aggregates, each having an average of 966 intra-aggregate pores, was generated to produce a domain with total of 3,357,025 pores. By varying the poresizes, different aggregates with porosity values between 0.082 and
0.700 were obtained. We showed that increasing aggregate porosities did lead to higher values of the mass transfer coeﬃcient andmore tailing in the BTCs. Higher ﬂuid velocities within the aggregate domain also increased the mass transfer coeﬃcient, but
produced less tailing in the calculated BTCs. On the other hand,lower pore water velocities in the micropore domain produced
lower mass transfer coeﬃcients but more extensive tailing in the
BTCs. The velocity ratio between the aggregated and macro domain
could explain the magnitude of tailing observed in the BTCs.
We showed that dual-porosity pore network models are attractive tools for investigating the basic processes governing solute interactions between the macropore domain (containing interaggregate pores) and the micropore or aggregate domain (containing intra-aggregate pores). The developed pore network model
provides a ﬂexible means of analysing the effects of pore scaletransport properties such as aggregate porosity and permeability
on macroscopic ﬂow and transport parameters of dual-porosity
porous media.
While in this study we focused on the effects of aggregate
porosity and permeability on macroscopic transport, further studies are needed to better understand all of the pore scale properties of aggregates. In our work we varied the aggregate porosity
by varying the size of the pore bodies. The aggregate porosity canvary also as a result of changes in the number of pore bodies, thepore coordination number, and the size of the throats. The inﬂuence of the range, variance, skewness, and curtosis of the pore size
distribution on the mass transfer coeﬃcients must also be investigated further. We note also that aggregate domain porosity canbe inﬂuenced by the size and number of aggregates, which havebeen kept constant in this study. We furthermore did not consider
any direct connections between the aggregates. These connections
will have an inﬂuence on the residence time of the solute in theaggregated domain. They very likely will inﬂuence the magnitudeof transferred mass, and how solutes are exchanged between the
micropore and macropore domains. The use of suitable imaging
techniques remains indispensable in all of this to most accurately
obtain the exact micropore structures of the aggregates, and theporous medium as a whole.
Acknowledgements
This work is part of the Veni Talent Scheme awarded to A.
Raoof with project number 016.151.047, which is (partly) ﬁnancedby the Netherlands Organisation for Scientiﬁc Research (NWO).
Supplementary materials
Supplementary material associated with this article can befound, in the online version, at doi:10.1016/j.advwatres.2017.04.013.

94

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95

References
Abbaspour, K.C., Kohler, A., Šimunek, J., Fritsch, M., Schulin, R., 2001. Applicationof a two-dimensional model to simulate ﬂow and transport in a macroporous
agricultural soil with tile drains. Eur. J. Soil Sci. 52, 433–447. http://dx.doi.org/
10.1046/j.1365-2389.20 01.0 0389.x.
Allaire, S.E., Roulier, S., Cessna, A.J., 2009. Quantifying preferential ﬂow in soils: a review of different techniques. J. Hydrol. 378, 179–204. http://dx.doi.org/10.1016/
j.jhydrol.2009.08.013.
Barker, J.A., 1985. Block-geometry functions characterizing transport in densely ﬁssured media. J. Hydrol. 77, 263–279.
Bolt, G.H., 1979. Movement of solutes in soil: principles of adsorption/exchange
chromatography. In: Bolt, G.H. (Ed.), Soil Chemistry, B. Physico-Chemical Models.
Developments in Soil Science 5B. Elsevier, Amsterdam, pp. 295–348.
Bijeljic, B., Mostaghimi, P., Blunt, M.J., 2013. Insights into non-Fickian solute transport in carbonates. Water Resour. Res. 49 (5), 2714–2728. http://dx.doi.org/10.
1002/wrcr.20238.
Brusseau, M.L., 1993. The inﬂuence of solute size, pore water velocity, and intraparticle porosity on solute dispersion and transport in soil. Water Resour. Res.
29 (4), 1071–1080.
Bultreys, T., Van Hoorebeke, L., Cnudde, V., 2015. Multi-scale, micro-computed
tomography-based pore network models to simulate drainage in heterogeneousrocks. Adv. Water Resour. 78, 36–49. http://dx.doi.org/10.1016/j.advwatres.2015.
02.003.
Bultreys, T., De Boever, W., Cnudde, V., 2016. Imaging and image-based ﬂuid transport modeling at the pore scale in geological materials: a practical introduction
to the current state-of-the-art. Earth-Sci. Rev. 155, 93–128. http://dx.doi.org/10.
1016/j.earscirev.2016.02.001.
Cnudde, V., Boone, M.N., 2013. High-resolution X-ray computed tomography in geosciences: a review of the current technology and applications. Earth-Sci. Rev.
123, 1–17. http://dx.doi.org/10.1016/j.earscirev.2013.04.003.
Czachor, H., Charytanowicz, M., Gonet, S., Niewczas, J., Jozefaciuk, G., Lichner, L.,2015. Impact of long-term mineral and organic fertilizer application on the water stability, wettability and porosity of aggregates obtained from two loamy
soils. Eur. J. Soil Sci. 66 (3), 577–588. http://dx.doi.org/10.1111/ejss.12242.
Clothier, B.E., Kirkham, M.B., McLean, J.E., 1992. In situ measurements of the effective transport volume of solute moving through soil. Soil Sci. Soc. Am J. 56,
733–736.
Coats, K.H., Smith, B.D., 1964. Dead end pore volume and dispersion in porous media. Soc. Petrol. Eng. J. 4, 73–84.
Dal Ferro, N., Berti, A., Francioso, O., Ferrari, E., Matthews, G.P., Morari, F., 2012. Investigating the effects of wettability and pore size distribution on aggregate stability: the role of soil organic matter and the humic fraction. Eur. J. Soil Sci. 63(2), 152–164. http://dx.doi.org/10.1111/j.1365-2389.2012.01427.x.
Dal Ferro, N., Charrier, P., Morari, F., 2013. Dual-scale micro-CT assessment of soil
structure in a long-term fertilization experiment. Geoderma 204, 84–93. http://dx.doi.org/10.1016/j.geoderma.2013.04.012.
Fathi, H., Raoof, A., Mansouri, S.H., van Genuchten, M.T., 2017a. Effects of porosity
and water saturation on the effective diffusivity of a cathode catalyst layer. J.
Electrochem. Soc. 164 (4), F298–F305. http://dx.doi.org/10.1149/2.0871704jes.
Fathi, H., Raoof, A., Mansouri, S.H., 2017b. Insights into the role of wettability incathode catalyst layer of proton exchange membrane fuel cell; pore scale immiscible ﬂow and transport processes. J. Power Sources 346, 57–67. http://dx.
doi.org/10.1016/j.jpowsour.2017.03.012.
Flury, M., Flühler, H., Jury, W.A., Leuenberger, J., 1994. Susceptibility of soils to preferential ﬂow of water: a ﬁeld study. Water Resour. Res. 30, 1945–1954.
Gerke, H.H., 2006. Preferential ﬂow descriptions for structured soils. J. Plant Nutr.
Soil Sci. 169, 382–400. http://dx.doi.org/10.1002/jpln.200521955.
Goltz, M.N., Roberts, P.V., 1987. Using the method of moments to analyze three-dimensional diffusion-limited solute transport from temporal and spatial perspectives. Water Resour. Res. 23 (8), 1575–1585.
Hantush, M.M., Mariño, M.A., 1988. Interlayer diffusive transfer and transport ofcontaminants in stratiﬁed formation. I. Theory. J. Hydrogeol. Eng., ASCE. 3 (4),
232–240.
Hendrickx, J.M.H., Flury, M., 2001. Uniform and preferential ﬂow mechanisms in thevadose zone. In: Conceptual Models of Flow and Transport in the Fractured Vadose Zone. National Research Council, National Academy Press, Washington, DC,
pp. 149–187.
Hemes, S., Desbois, G., Urai, J.L., Schröppel, B., Schwarz, J.O., 2015. Multi-scalecharacterization of porosity in Boom Clay (HADES-level, Mol, Belgium) using
a combination of X-ray μ-CT, 2D BIB-SEM and FIB-SEM tomography. Microporous Mesoporous Mater. 208, 1–20. http://dx.doi.org/10.1016/j.micromeso.2015.
01.022.
Hillel, D., 2003. Introduction to Environmental Soil Physics. Academic Press.
Jafari, S., Yamamoto, R., Rahnama, M., 2011. Lattice-Boltzmann method combinedwith smoothed-proﬁle method for particulate suspensions. Phys. Rev. E. 83 (2).
http://dx.doi.org/10.1103/PhysRevE.83.026702.
Jarvis, N.J., 2007. A review of nonequilibrium water ﬂow and solute transport in soil
macropores: principles, controlling factors and consequences for water quality.
Eur. J. Soil Sci. 58, 523–546. http://dx.doi.org/10.1111/j.1365-2389.20 07.0 0915.x.
Jaynes, D.B., Shao, M., 1999. Evaluation of a simple technique for estimating two–domain transport parameters. Soil Sci. 164, 82–91.
Jiang, Z., Dijke, M.I.J., Sorbie, K.S., Couples, G.D., 2013. Representation of multiscale
heterogeneity via multiscale pore networks. Water Resour. Res. 49 (9), 5437–
5449. http://dx.doi.org/10.1002/wrcr.20304.

Khan, F., Enzmann, F., Kersten, M., Wiegmann, A., Steiner, K., 2012. 3D simulation of
the permeability tensor in a soil aggregate on basis of nanotomographic imaging and LBE solver. J. Soils Sediments 12 (1), 86–96. http://dx.doi.org/10.1007/
s11368-011-0435-3.
Köhne, J.M., Köhne, S.K., Šimůnek, J., 2009. A review of model applications for structured soils: a) Water ﬂow and tracer transport. J. Contam. Hydrol. 104, 4–35.
http://dx.doi.org/10.1016/j.jconhyd.2008.10.002.
Mahmoodlu, M.G., Hartog, N., Hassanizadeh, S.M., Raoof, A., 2013. Oxidation ofvolatile organic vapours in air by solid potassium permanganate. Chemosphere91 (11), 1534–1538. http://dx.doi.org/10.1016/j.chemosphere.2012.12.035.
Mahmoodlu, M.G., Hassanizadeh, S.M., Hartog, N., Raoof, A., 2014. Oxidation oftrichloroethylene, toluene, and ethanol vapors by a partially saturated permeable reactive barrier. J. Contam. Hydrol. 164, 193–208. http://dx.doi.org/10.1016/
j.jconhyd.2014.05.013.
Mangalassery, S., Sjögersten, S., Sparkes, D.L., Sturrock, C.J., Mooney, S.J., 2013. The
effect of soil aggregate size on pore structure and its consequence on emissionof greenhouse gases. Soil Till. Res. 132, 39–46. http://dx.doi.org/10.1016/j.still.
2013.05.003.
Maraqa, M.A., 2001. Prediction of mass-transfer coeﬃcient for solute transportin porous media. J. Contam. Hydrol. 53, 163–171. http://dx.doi.org/10.1016/
S0169-7722(01)00107-3.
Martínez, F.S.J., Ortega, F.M., Monreal, F.C., Kravchenko, A.N., Wang, W., 2015. Soilaggregate geometry: measurements and morphology. Geoderma 237, 36–48.
http://dx.doi.org/10.1016/j.geoderma.2014.08.003.
Mehmani, A., Prodanović, M., 2014. The effect of microporosity on transport properties in porous media. Adv. Water Resour. 63, 104–119. http://dx.doi.org/10.1016/
j.advwatres.2013.10.009.
Mehmani, A., Mehmani, Y., Prodanović, M., Balhoff, M., 2015. A forward analysis on
the applicability of tracer breakthrough proﬁles in revealing the pore structureof tight gas sandstone and carbonate rocks. Water Resour. Res. 51 (6), 4751–
4767. http://dx.doi.org/10.1002/2015WR016948.
NRC, 2001. Conceptual Models of Flow and Transport in the Fractured Vadose Zone.
National Research Council, National Academy Press, Washington, DC.
Parker, J.C., Valocchi, A.J., 1986. Constraints on the validity of equilibrium andﬁrst-order kinetic transport models in structured soils. Water Resour. Res. 22(3), 399–407.
Peth, S., Horn, R., Beckmann, F., Donath, T., Fischer, J., Smucker, A.J.M., 2008.
Three-dimensional quantiﬁcation of intra-aggregate pore-space features usingsynchrotron-radiation-based microtomography. Soil Sci. Soc. Am. J. 72 (4), 897–
907. http://dx.doi.org/10.2136/sssaj2007.0130.
Pontedeiro, E.M., van Genuchten, M.Th., Cotta, R.M., Šimůnek, J., 2010. The effects ofpreferential ﬂow and soil texture on risk assessments of a NORM waste disposalsite. J. Hazard. Mater 174, 648–655. http://dx.doi.org/10.1016/j.jhazmat.2009.09.
100.
Prodanović, M., Mehmani, A., Sheppard, A.P., 2015. Imaged-based multiscale network modelling of microporosity in carbonates. Geol. Soc. Lond. Special Publications 406 (1), 95–113. http://dx.doi.org/10.1144/SP406.9.
Rabbi, S.F., Wilson, B.R., Lockwood, P.V., Daniel, H., Young, I.M., 2015. Aggregate hierarchy and carbon mineralization in two Oxisols of New South Wales, Australia.
Soil Till. Res. 146, 193–203. http://dx.doi.org/10.1016/j.still.2014.10.008.
Raoof, A., Hassanizadeh, S.M., 2010. A new method for generating pore-networkmodels of porous media. Transp. Porous Media 81 (3), 391–407. http://dx.doi.
org/10.10 07/s11242-0 09-9412-3.
Raoof, A., Hassanizadeh, S.M., Leijnse, A., 2010. Upscaling transport of adsorbing solutes in porous media: pore-network modeling. Vadose Zone J. 9 (3), 624–636.
http://dx.doi.org/10.2136/vzj2010.0026.
Raoof, A., Hassanizadeh, S.M., 2012. A new formulation for pore-network modeling of two-phase ﬂow. Water Resour. Res. 48 (1). http://dx.doi.org/10.1029/
2010WR010180.
Raoof, A., Hassanizadeh, S.M., 2013. Saturation-dependent solute dispersivity inporous media: pore-scale processes. Water Resour. Res. 49 (4), 1943–1951.
http://dx.doi.org/10.1002/wrcr.20152.
Raoof, A., Nick, H.M., Hassanizadeh, S.M., Spiers, C.J., 2013. PoreFlow: a complex
pore-network model for simulation of reactive transport in variably saturated
porous media. Comp. Geosci. 61, 160–174. http://dx.doi.org/10.1016/j.cageo.2013.
08.005.
Rasmuson, A., Neretnieks, I., 1980. Exact solution of a model for diffusion in particles and longitudinal dispersion in packed beds. Am. Inst. Chem. Eng. J. 26,
686–690.
Sahimi, M., 2011. Flow and Transport in Porous Media and Fractured Rock: FromClassical Methods to Modern Approaches. Wiley-VCH. http://dx.doi.org/10.1002/
9783527636693.
Shukla, M.K., 2013. Soil Physics: An Introduction. CRC Press.
Šimůnek, J., Jarvis, N.J., van Genuchten, M.Th., Gärdenäs, A., 2003. Review and comparison of models for describing non-equilibrium and preferential ﬂow andtransport in the vadose zone. J. Hydrol. 272 (1), 14–35. http://dx.doi.org/10.1016/
S0 022-1694(02)0 0252-4.
Šimunek, J., van Genuchten, M.Th., Šejna, M., Toride, N., Leij, F.J., 20 0 0. The STANMOD computer software for evaluating solute transport in porous media using analytical solutions of convection-dispersion equation. In: Version 2.0, IGWMC-TPS-71, Int. Ground Water Modeling Center (IGWMC). Colorado School ofMines, Golden, Colorado, p. 32.
Šimunek, J., van Genuchten, M.Th, 2008. Modeling nonequilibrium ﬂow and transport processes using HYDRUS. Vadose Zone J. 7 (2), 782–797. http://dx.doi.org/
10.2136/vzj20 07.0 074.

E.T. de Vries et al. / Advances in Water Resources 105 (2017) 82–95
Sudicky, E.A., Frind, E.O., 1982. Contaminant transport in fractured porous media:
analytical solutions for a system of parallel fractures. Water Resour. Res. 18 (6),1634–1642.
Tang, D.H., Frind, E.O., Sudicky, E.A., 1981. Contaminant transport in fractured porous
media: analytical solution for a single fracture. Water Resour. Res. 17, 555–564.
Tisdall, J.M., Oades, J., 1982. Organic matter and water-stable aggregates in soils. J.
Soil Sci. 33 (2), 141–163.
Toride, N., Leij, F.J., van Genuchten, M.Th, 1995. The CXTFIT code for estimatingtransport parameters from laboratory or ﬁeld tracer experiments. In: Version
2.10. U.S. Salinity Laboratory. USDA-ARS, Riverside, CA, p. 121.
Vasilyev, L., Raoof, A., Nordbotten, J.M., 2012. Effect of mean network coordination
number on dispersivity characteristics. Transp. Porous Media 95 (2), 447–463.
http://dx.doi.org/10.1007/s11242-012-0054-5.
van Genuchten, M.Th., Dalton, F.N., 1986. Models for simulating salt movement inaggregated ﬁeld soils. Geoderma 38 (1), 165–183.
van Genuchten, M.Th, 1985. General approach for modeling solute transport instructured soils. Memoires Int. Assoc. Hydrogeol. 17 (2), 513–526.
van Genuchten, M.Th., Wierenga, P.J., 1976. Mass transfer studies in sorbing porous
media, I. Analytical solutions. Soil Sci. Soc. Am. J. 40 (4), 473–480.

95

van Genuchten, M.Th., Tang, D.H., Guennelon, R., 1984. Some exact solutions for solute transport through soils containing large cylindrical macropores. Water Resour. Res. 20 (3), 335–346.
Vogel, T., Sanda, M., Dusek, J., Dohnal, M., Votrubova, J., 2010. Using oxygen-18 to
study the role of perefential ﬂow in the formation of hillslope runoff. Vadose
Zone J. 9, 252–259. http://dx.doi.org/10.2136/vzj20 09.0 066.
Wildenschild, D., Sheppard, A.P., 2013. X-ray imaging and analysis techniques forquantifying pore-scale structure and processes in subsurface porous medium
systems. Adv. Water Resour. 50, 217–246. http://dx.doi.org/10.1016/j.advwatres.
2012.07.018.
Zhang, Q., Hassanizadeh, S.M., Karadimitriou, N.K., Raoof, A., Liu, B., Kleingeld, P.J.,
Imhof, A., 2013. Retention and remobilization of colloids during steady-state andtransient two-phase ﬂow. Water Resour. Res. 49 (12), 8005–8016. http://dx.doi.
org/10.1002/2013WR014345.
Zhou, H., Peng, X., Perfect, E., Xiao, T., Peng, G., 2013. Effects of organic and inorganic
fertilization on soil aggregation in an Ultisol as characterized by synchrotron
based X-ray micro-computed tomography. Geoderma 195, 23–30. http://dx.doi.
org/10.1016/j.geoderma.2014.08.003.
Zhou, B.B., Li, Y., Wang, Q., Jiang, Y.L., Li, S., 2014. Simulation of chloride transport inan aggregated soil using three conceptual models. Arab. J. Geosci. 7 (7), 2539–
2546. http://dx.doi.org/10.1007/s12517- 013- 0970- x.