Microsoft Word - 43poletti.docx


 CHEMICAL ENGINEERING TRANSACTIONS  
 

VOL. 47, 2016 

A publication of 

 
The Italian Association 

of Chemical Engineering 
Online at www.aidic.it/cet 

Guest Editors: Angelo Chianese, Luca Di Palma, Elisabetta Petrucci, Marco Stoller
Copyright © 2016, AIDIC Servizi S.r.l., 
ISBN 978-88-95608-38-9; ISSN 2283-9216 

Microscale Colloidal Transport Simulations for Groundwater 
Remediation 

Eleonora Crevacorea, Gianluca Boccardob, Daniele L. Marchisiob, Rajandrea 
Sethi*a 
a
DIATI, Politecnico di Torino, C.so Duca degli Abruzzi 24, 10129 Torino 

b
DISAT, Politecnico di Torino, C.so Duca degli Abruzzi 24, 10129 Torino 

rajandrea.sethi@polito.it 

The injection of zero valent iron particles is a promising approach for the remediation of contaminated aquifer 
systems. In this framework microscale numerical simulations are a powerful tool to investigate transport and 
deposition of colloidal suspensions. The present analysis compare different porous media structures and 
highlights the role played by the geometry of the medium on the physics of the problem. Extrapolating 
macroscale parameters of interest in practical applications, results are compared with classical literature 
correlations, that are inadequate to describe the effect of a complex 3D structure. 

1. Introduction 
The investigation of transport and deposition of colloidal particles in porous media is a key factor in several 
fields, from environmental engineering to medicine. Just to make a few examples, the colloid filtration theory 
(CFT) finds application in catalytic processes (Bensaid et al., 2010), investigation of clogging of depth filters 
and wells, enhanced oil recovery (Shi et al., 2013) and thermo-radiotherapy, (Gordon et al., 2014). 
In the field of groundwater remediation the injection of colloidal nano and micro zero-valent iron particles (nZVI 
and mZVI) is a very promising technique. Indeed iron particles are able to degrade a broad range of 
contaminants, including the most recalcitrant ones, (Di Molfetta and Sethi, 2006, Comba et al., 2011). Zero-
valent iron particles can be injected close to the source of contamination and their delivery to the 
contaminated plume can enhanced by groundwater flow, (Kocur et al., 2014, Krol et al., 2013, O’Carroll et al., 
2013, Vecchia et al., 2009, Gastone et al., 2014, Tosco et al., 2014, Luna et al., 2015). In order to design a 
successful intervention, an intimate understanding of the phenomena governing the problem is needed. In this 
perspective, microscale simulations make for a powerful tool to simulate colloidal transport and deposition 
under a broad range of operative conditions that allow to understand, develop and derive constitutive laws and 
design macro-scale models, (Boccardo et al., 2014, Messina et al., 2015, Icardi et al., 2014). 
In colloid filtration theory there is an extended literature concerning the transport and deposition of colloidal 
suspensions; the focus of these works is on single spherical collectors, (Happel, 1959, Nelson and Ginn, 
2005, Nelson and Ginn, 2011, Tufenkji and Elimelech, 2004, Yao et al., 1971, Johnson et al., 2007, Song and 
Elimelech, 1995, Molnar et al., 2015). The appeal of this classical approach lies in the use of a detailed 
description of the physics of the problem, since it is possible to consider the mutual interaction of different 
transport mechanism, like Brownian diffusion, interception, and sedimentation. On the other hand because of 
the simplified representation of the porous medium many simplifications are introduced; the starting point to 
overcome these simplifications is the analysis of full 3D geometries made of spherical collectors. 
In this work we present a preliminary study of two distinct periodic porous media: this kind of geometries is 
suitable to reproduce homogeneous and isotropic media and represents an intermediate step between 
simplified geometries and realistic media. Numerical simulation of particle injection were performed, and from 
the comparison of the results in different porous media geometries it is possible to highlight that, not simply 
the presence of a multiplicity of collectors, but mainly their spatial configuration greatly affect transport and 
deposition mechanism of colloidal suspensions. As a consequence, a particular attention is devoted to the 

                               
 
 

 

 
   

                                                  
DOI: 10.3303/CET1647046

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Please cite this article as: Crevacore E., Boccardo G., Marchisio D., Sethi R., 2016, Microscale colloidal transport simulations for groundwater 
remediation, Chemical Engineering Transactions, 47, 271-276  DOI: 10.3303/CET1647046

271



proper evaluation of the deposition efficiency, as it is a macroscale parameter of interest in the design of 
remediation technologies.  

2. Simulation approach 
In the present study, the porous media structure is modelled with a 3D arrangement of spherical collectors of 
equal size, held fixed in ordered cubic packing structures. In particular two configurations inspired to 
crystallography were chosen: the face-centered cubic (FCC) and the body-centered cubic (BCC) structures, 
Figure(1). These geometries are characterized by different values of the specific surface (i.e.: the ratio of the 
collectors surface over the total volume of the medium); moreover they have, at the center of the cubic 
element, three plane of symmetry and this geometrical peculiarity can be exploited to reduce the 
computational cost of numerical simulations, since just one fourth of the cubic element is sufficient to describe 
the whole geometry. 
The starting point of our analysis is the determination of the fluid velocity field within the porous medium; in 
particular we solved the Navier-Stokes and the continuity equations, that under steady state conditions read: ∙ = − + ∆  , (1) 
∙ = 0 , (2) 

where ρ is the water density (kg·m-3), u is the effective velocity (m·s-1), p is the pressure (kg·m-1·s-2)) and μ is 
the dynamic viscosity (kg·m-1·s-1)). To solve these equations, proper boundary conditions are needed. The 
dimensionless parameter associated to Eq(1) is the Reynolds number; it is the ratio between inertial and 
viscous forces. When Re << 1, viscous forces are predominant and the flow is laminar (creeping flow); under 
this condition the Darcy law holds true and there is a linear relationship between the superficial velocity 
discharge and the pressure gradient. On the grain surface the no-slip condition is imposed, that is the fluid 
velocity is equal to zero. In the Eulerian framework, the transport of colloidal particles can be described with 
the advection diffusion equation: , + ∙ , = ∙ ,  , (3) 
where c is the solute concentration per volume of liquid (kg·m-3) and Dm  is the molecular diffusion coefficient 
(m2·s-1). The underlying assumption is that the solute concentration is sufficiently small so that its transport 
can be described with Fick’s law and the solute itself is treated as a colloid, that is a diluted suspension of 
spherical particles. It is generally possible to evaluate the molecular diffusion with the Stokes-Einstein 
equation: = 3 	, (4) 
where dp is the colloidal particle diameter (m). In this work no external forces, like gravity and van der Waals, 
were considered. To solve Eq(3) both initial and boundary conditions are needed. To simulate a continuous 
injection, a constant concentration is imposed at the inlet face during the whole simulation. Two different type 
of phenomena can be simulated, depending on the conditions imposed on the grains surface: to simulate the 
transport of an conservative tracer, the concentration gradient, whereas to simulate an instantaneous 
heterogeneous reaction (i.e. perfect sink condition), concentration equal to zero is imposed. The 
dimensionless number associated to Eq(3) is the Péclet number, describing the magnitude of advection with 
respect to diffusion. 
In the framework of the colloid filtration theory, the grains composing the porous medium are seen as 
collectors, and the deposition and removal of particles is described as a two-step process: particle transport 
from the bulk of the fluid to the collector surface, and attachment. These steps are quantitatively described by 
the collector efficiency η  and the attachment efficiency α; the total removal efficiency is then defined as η = αη . In presence of favorable physical-chemical conditions,α is equal to one and total efficiency is simply 
equal to η0. When Brownian diffusion is the leading mechanism, (Levich and Technica, 1962), defined η as the 
ratio between the flux normal to the grains surface (diffusive flux) and that normal to the inlet face (diffusive 
and advective):  
 

272



 

Figure 1: Face-centered (top) and body-centered (bottom) packing structures, with flow-filed contour plots. 

	=	 ∇ ( , ) ∙ 	 			 ( )	 ( , ) ∙ + ∇ ( , ) ∙ . (7) 
Then, the following correlation valid for a single spherical collector immersed in an infinite fluid medium, were 
proposed: η	=	4.04	 / . (8) 
Starting from this, many other authors presented new relations that try to include the influence of other 
collectors present is the system. A very popular model is Happel’s, (Happel, 1959), that reads as: η	=	4.04	 / 	 / 	, (9) 
where 	 = (1 − ) /  and =	 	  . 
Dealing with media composed by multiple collectors, neither the single collector efficiency as derived by 
Levich nor the one proposed by Happel can be directly applied to evaluate the total efficiency of the medium. 
In this work, no other deposition mechanism, like interception and sedimentation will be considered. 

3. Test cases and numerical details 
All the simulations presented in this work have been performed with the finite volume open-source CFD code 
OpenFOAM. 
We built a computational domain for each of the presented geometrical configurations. In order to do that we 
chose a reasonable grains diameter dg=200 µm and porosity ε=0.4 for both cases; the latter in order to avoid 
solid compenetration. As anticipated, we took advantage of structure symmetries, and considered a domain 
equal to one fourth of the cubic cell. Moreover, to guarantee a proper spatial evolution of the phenomena 
investigated, the computational domain is formed by a series of four elementary modules. It is important to 
notice that the computational domain is periodic, with a period equal to the side of the cubic structure, as 
detailed in Table(1). We assume that the medium is completely saturated with water at ambient temperature. 
Pressure drop values were chosen in order to obtain effective velocities ranging from 1·10-5 to 1·10-2 m·s-1; the 
corresponding Reynolds number are presented in Table(2).  
The resulting flow fields were used in both type of transport simulations, both with and without heterogeneous 
reaction. For each geometry we chose two different type of monodisperse colloidal suspensions, with colloidal 
particle diameters of 100 nm and 400 nm respectively, resulting in molecular diffusion values of 4.392·10-12 
and 1.098·10-12m2·s-1 respectively. In order to reach a steady state concentration, the length of each 
simulation is equal to five characteristic advective times τ (i.e.: the ratio between the domain length and the 
superficial velocity). 
For the following analysis, it is important to notice that the main flow direction is along the x-axis, and so it is 
the meaningful particle concentration gradient. 
The concentrations obtained from the performed simulations are microscale values. Since the final purpose is 
to extract information to be used in macroscale models, we extrapolated macroscale values of concentration; 
in particular we are interested in influent and effluent concentrations. Given the domains configuration, these 
values are defined as: 

273



( , ) =	1 , 	 , (12) 
where C is the macroscale concentration, c is the microscale concentration and A is the integration area at the 
inlet-outlet of each module; in particular for x=0 we obtained the inlet concentration, C0, and for x = 4L, we 
have the outlet concentration, Cout. 

Table 1:  Geometrical characteristics of computational domains. 

Structure Chosen porosity Grains diameter (µm) Period length (µm) Total domain length (µm) 
FCC 0.4  200 303.4  1213.6 
BCC 0.4 200 240.8  963.2 

 

Table 2:  Operating conditions: effective velocities and dimensionless groups. 

 FCC  BCC  
Effective 
velocity 
(m·s-1) 

Re Pe 
D= 4.392·10-12 

Pe 
D= 1.098·10-12 

 (s) Re Pe 
D= 4.392·10-12 

Pe 
D= 1.098·10-12 

 (s) 

1·10-2 8.21·10-1 2.78·106 1.11·107 0.120 8.17·10-1 2.19·106 8.77·106 0.096 
1·10-3 8.21·10-2 2.78·105 1.11·106 1.207 8.17·10-2 2.19·105 8.77·105 0.963 
1·10-4 8.21·10-3 2.78·104 1.11·105 12.071 8.17·10-3 2.19·104 8.77·104 9.631 
1·10-5 8.21·10-4 2.78·103 1.11·104 120.711 8.17·10-4 2.19·103 8.77·103 96.315

4. Results 
A first meaningful evaluation that highlights the influence of the geometrical configuration on transport 
phenomena comes from the comparison of the breakthrough curves, representing the temporal evolution of 
the effluent concentration Cout, normalized with respect to the influent one, C0. Concerning tracer test 
simulations, as depicted in Figure(2), the breakthrough curves show a very different behaviours depending on 
the geometry. In particular, while in the face-centered case the curve presents a sigmoidal behaviour and the 
dispersion flux can be modelled with Fick’s law, this is not true in the body-centered case. In the BCC 
configuration at the beginning of injection there is a swift increase of concentration, while after a time equal to τ, concentration rise with time slows considerably. 
As stated before, to design macroscale models the knowledge of the deposition efficiency is of pivotal 
importance. In particular, we employed Eq(7) in order to evaluate the total efficiency for each geometrical 
structure. Moreover we looked for the spatial variation of , evaluating it over each of the four modules 
composing our periodic geometries. From the comparison of the total deposition efficiency of face-centered 
and body-centered geometry, the FCC show higher values of . This difference can be ascribed to a different 
effective porosity characterizing the geometry. That is, since the influence of the packing structure on the flow 
field is not entirely negligible, to properly evaluate transport phenomena, the dimension of the portion of fluid in 
which the advective transport does not take place should be carefully considered. In particular, while in the 
face-centered geometry (in the proposed range of operating conditions) the effective porosity is equal to the 
geometric one, in the body-centered case it is lower. In order to compare the values obtained with literature 
results, the efficiency obtained were normalized with respect to the number of collectors present in each 
configuration, Figure(3). The values for both the geometries collapse over each other and while they are 
sufficiently close to theoretical values, it is possible to appreciate that the trend is different and less regular. 
At last, as depicted in Figure(4), we provide the evolution of  respect to x; in particular it is evaluated in four 
points, that correspond to the outlet of each periodic module composing the medium. As expected, this 
analysis 

274



 

Figure 2: Tracer simulation – comparison of breakthrough curves of FCC (continuos line) and BCC (dashed 
line) configuration, Pèclet number of 2.78·104. 

shows that the efficiency decreases with increasing distance from the injection point. The differences between 
efficiency at the first module and that at the last one are slightly accentuated in the body-centered 
configuration, but they mainly depend on the operating conditions, that is on the extension and shape of the 
concentration front; indeed greater differences are found at low Pèclet numbers. The visualisation of the 
contour plot of concentration, Figure(4), gives a clear idea of the phenomena. 

5. Conclusions 
The present work highlights the role of the distribution of the grains composing the porous medium on the 
transport of dilutes and the deposition of particles. From the analysis of transport it emerges that the 
distribution of the grains composing the medium greatly affects the breakthrough of colloids. This effect can be 
seen also on the parameter η, where multiple-collector structures with the same porosity but different grain 
distribution result in a different deposition efficiency; when a normalization of these values over the specific 
surface area is performed, however, the results collapse on a single curve. 
 

 

Figure 3: On the left, comparison between total efficiency for the FCC (continuos line) and the BCC (dashed 
line) geometry; on the right, the single collector efficiency of Happel’s model (dotted line), evaluated at ε=0.4, 
is compared to the normalized efficiency for both geometries. 

 
Figure 4: Spatial variation of the efficiency coefficient for the FCC case (top) and the BCC case (bottom) 
 
 
 
 

275



Reference 

 
BENSAID, S., MARCHISIO, D. L. & FINO, D. 2010. Numerical simulation of soot filtration and combustion 

within diesel particulate filters. Chemical Engineering Science, 65, 357-363. 
BOCCARDO, G., MARCHISIO, D. L. & SETHI, R. 2014. Microscale simulation of particle deposition in porous 

media. Journal of colloid and interface science, 417, 227-237. 
COMBA, S., DI MOLFETTA, A. & SETHI, R. 2011. A comparison between field applications of nano-, micro-, 

and millimetric zero-valent iron for the remediation of contaminated aquifers. Water, Air, & Soil Pollution, 
215, 595-607. 

DI MOLFETTA, A. & SETHI, R. 2006. Clamshell excavation of a permeable reactive barrier. Environmental 
Geology, 50, 361-369. 

GASTONE, F., TOSCO, T. & SETHI, R. 2014. Green stabilization of microscale iron particles using guar gum: 
bulk rheology, sedimentation rate and enzymatic degradation. Journal of colloid and interface science, 
421, 33-43. 

GORDON, A. C., LEWANDOWSKI, R. J., SALEM, R., DAY, D. E., OMARY, R. A. & LARSON, A. C. 2014. 
Localized Hyperthermia with Iron Oxide–Doped Yttrium Microparticles: Steps toward Image-Guided 
Thermoradiotherapy in Liver Cancer. Journal of Vascular and Interventional Radiology, 25, 397-404. 

HAPPEL, J. 1959. Viscous flow relative to arrays of cylinders. AIChE Journal, 5, 174-177. 
ICARDI, M., BOCCARDO, G., MARCHISIO, D. L., TOSCO, T. & SETHI, R. 2014. Pore-scale simulation of 

fluid flow and solute dispersion in three-dimensional porous media. Physical Review E, 90, 013032. 
JOHNSON, W., LI, X. & YAL, G. 2007. Colloid retention in porous media: Mechanistic confirmation of wedging 

and retention in zones of flow stagnation. Environmental science & technology, 41, 1279-1287. 
KOCUR, C. M., CHOWDHURY, A. I., SAKULCHAICHAROEN, N., BOPARAI, H. K., WEBER, K. P., SHARMA, 

P., KROL, M. M., AUSTRINS, L., PEACE, C. & SLEEP, B. E. 2014. Characterization of nZVI mobility in a 
field scale test. Environmental science & technology, 48, 2862-2869. 

KROL, M. M., OLENIUK, A. J., KOCUR, C. M., SLEEP, B. E., BENNETT, P., XIONG, Z. & O’CARROLL, D. M. 
2013. A field-validated model for in situ transport of polymer-stabilized nZVI and implications for 
subsurface injection. Environmental science & technology, 47, 7332-7340. 

LEVICH, V. G. & TECHNICA, S. 1962. Physicochemical hydrodynamics, Prentice-hall Englewood Cliffs, NJ. 
LUNA, M., GASTONE, F., TOSCO, T., SETHI, R., VELIMIROVIC, M., GEMOETS, J., MUYSHONDT, R., 

SAPION, H., KLAAS, N. & BASTIAENS, L. 2015. Pressure-controlled injection of guar gum stabilized 
microscale zerovalent iron for groundwater remediation. Journal of Contaminant Hydrology. 

MESSINA, F., MARCHISIO, D. L. & SETHI, R. 2015. An extended and total flux normalized correlation 
equation for predicting single-collector efficiency. Journal of colloid and interface science, 446, 185-193. 

MOLNAR, I. L., JOHNSON, W. P., GERHARD, J. I., WILLSON, C. S. & O'CARROLL, D. M. 2015. Predicting 
colloid transport through saturated porous media: A critical review. Water Resources Research, 51, 6804-
6845. 

NELSON, K. E. & GINN, T. R. 2005. Colloid filtration theory and the Happel sphere-in-cell model revisited with 
direct numerical simulation of colloids. Langmuir, 21, 2173-2184. 

NELSON, K. E. & GINN, T. R. 2011. New collector efficiency equation for colloid filtration in both natural and 
engineered flow conditions. Water Resources Research, 47. 

O’CARROLL, D., SLEEP, B., KROL, M., BOPARAI, H. & KOCUR, C. 2013. Nanoscale zero valent iron and 
bimetallic particles for contaminated site remediation. Advances in Water Resources, 51, 104-122. 

SHI, X., PRODANOVIĆ, M., HOLDER, J., GRAY, K. & DICARLO, D. 2013. Coupled solid and fluid mechanics 
modeling of formation damage near wellbore. Journal of Petroleum Science and Engineering, 112, 88-96. 

SONG, L. & ELIMELECH, M. 1995. Particle deposition onto a permeable surface in laminar flow. Journal of 
colloid and interface science, 173, 165-180. 

TOSCO, T., GASTONE, F. & SETHI, R. 2014. Guar gum solutions for improved delivery of iron particles in 
porous media (Part 2): Iron transport tests and modeling in radial geometry. Journal of Contaminant 
Hydrology, 166, 34-51. 

TUFENKJI, N. & ELIMELECH, M. 2004. Correlation equation for predicting single-collector efficiency in 
physicochemical filtration in saturated porous media. Environmental science & technology, 38, 529-536. 

VECCHIA, E. D., LUNA, M. & SETHI, R. 2009. Transport in porous media of highly concentrated iron micro-
and nanoparticles in the presence of xanthan gum. Environmental science & technology, 43, 8942-8947. 

YAO, K.-M., HABIBIAN, M. T. & O'MELIA, C. R. 1971. Water and waste water filtration. Concepts and 
applications. Environmental science & technology, 5, 1105-1112. 

 
 

276