Fundamental theoretical study of nanoscale ion transport through

charged porous rocks: Part I diffusion

P. Pivonka, D. Smith & B. Gardiner

Department of Civil and Environmental Engineering, The University of Melbourne, Parkville, Melbourne, Australia

ABSTRACT: Understanding of diffusive mass transport through porous rocks is essential for many applications
in geomechanics and environmental engineering (e.g., the assessment of contaminant transport through porous
rocks or the development of ore deposits). Such systems are generally characterized by aqueous solutions
containing multiple ions filling the pore space of the porous material. However, the surfaces of the porous rock
may be charged leading to anion/cation exclusion depending on the charge of the ion transported through the
solution.
In this paper, the fundamental equations describing coupled multi-ion diffusion and migration through charged
porous rocks are given for ideal solution systems (i.e., the Poisson-Nernst-Planck (PNP) equations. This system
of governing equations is applied at the nanoscale. The PNP equations are generalized for the description of mass
transfer in porous rocks using the volume averaging method applied to periodic media. This volume averaging
procedure accounts for the geometry, the porosity, and the tortuosity of the material. The influence of various
effects such as multiple ions, particle size, and size of surface charge on particles of the porous material on ion
transport behaviour under steady state conditions is examined in detail using numerical simulation tools.

1 INTRODUCTION we use the periodic media approach (Auriault, 1991),

(Whitaker, 1999). Investigation of a unit cell (repre-
In many fields of engineering steady state diffusive sentative volume element) allows for conclusions to
mass transport is described by Ficks first law (Cussler, be made about the macroscale behaviour. Two exam-
1997). However, this neglects electrical driving forces ple problems will be investigated: (i) diffussive mass
(which may be caused by the presence of an internal transport in uncharged porous media (pure solution
electrical field as a result of couplings between indi- system) and (ii) diffusive transport in charged porous
vidual ions and/or the presence of an external electrical media. For the charged porous medium the influence
field). These effects can be taken into account using of the microstructure, i.e., the particle size, shape,
the Nernst-Planck (N-P) equations (Newman, 1991). and orientation on diffusive mass transfer is investi-
In order to obtain a unique description of diffusive gated. Furthermore the influence of the particle charge
transport the N-P equations must be supplemented by density is studied.
an additional equation. The choice of this equation
is made on electrochemical rather than on a math-
ematical basis. The most rigorous way is to define
the variation of the electric potential according to the 2 TRANSPORT EQUATIONS FOR DIFFUSION
spacial distribution of electrical charges (Helfferich, AND MIGRATION
1962), (Newman, 1991). This relation leads to the
well known Poisson (P) equation of electrostatics. Use 2.1 Governing equations at the microscale
of this set of equations allows for the description of Diffusion is the movement of a species under a gradient
diffuse double-layers encountered in the vicinity of of concentration, or more correctly speaking, chemi-
charged particles such as charged clay particles. The cal potential, from regions of high chemical potential
presence of such diffuse double-layers in porous mate- towards regions of lower chemical potential. On the
rials may lead to anion/cation exclusion depending on other hand, migration is the movement of ionic species
the charge of ions in solution (Smith et al., 2002). along the gradient of an electric field. This two trans-
In this paper we apply the PNP system of equations port mechanisms can be combined by the notion of
to describe diffusive mass transfer. For this purpose the electrochemical potential. In the following we

417
Copyright 2005 Taylor & Francis Group plc, London, UK
derive the transport equations for ideal solutions (i.e., E is the electric field. Eqns. (4)(6) only hold for
we neglect the influence of activity). special situations depending on the problem of inves-
The constitutive equation describing multi-ion dif- tigation (for example electroneutrality is observed in
fusion in the presence of electrochemical forces is the all solutions except in thin double-layers near elec-
Nernst-Planck (NP) equation (Cussler, 1997), that is, trodes and charged particles in porous materials). The
most rigorous way is to relate the variation of the
electric potential to the spatial distribution of the elec-
tric charges. This relationship is given by the Poisson
equation (Helfferich, 1962):
where ji , ci , and ci is the molar flux density, the
concentration of the i-th ion, and the concentration gra-
dient of the i-th ion, respectively. Di is the free-solution
diffusion coefficient of the i-th ion. The quantities zi where the permittivity is the product of the permit-
and are the valence of the i-th ion and the voltage tivity of free space, 0 , and the relative permittivity
gradient. F, R and T are the Faraday constant, gas con- of water w . q is the electrical charge density. For
stant, and the absolute temperature, respectively. The the special case of a binary electrolyte (i.e., N = 2)
NP equation is a generalization of Ficks first law, tak- the influence of these assumptions on the numerical
ing into account electrical effects on charged solutes. results were reported in (Kato, 1995). It was found
In Eqn. (1) the voltage may arise from the application that the electro-neutrality condition (4) and the con-
of an external electrical field and/or from an inter- stant electric field assumption (5) are asymptotic cases
nal electrical field caused by interactions of ions. This of the Poisson Eqn. (7) (see (MacGillivray, 1968)
internal electrical field (also called diffusion potential) (Kato, 1995) for details). In the literature the coupled
is generated by different self-diffusion coefficients of Poisson-Nernst-Planck equations are often denoted as
individual ions, which leads to the formation of time PNP equations. Because of the complexity of the PNP
averaged ion dipoles during diffusion. However, the PDEs numerical methods are employed for their solu-
ions in solution are influenced by the electric potential tion. A comprehensive review on numerical models
irrespective of its origin, so the two voltages can be for solving the PNP system of equations is given in
summed and represented by a single quantity . (Samson et al., 1999).
The flux Eqn. (1) are supplemented by mass con- For the case of a binary electrolyte, i.e., a solution
servation equations. Mass conservation for each ion of a single salt composed of one kind of cation and one
species is expressed by, kind of anion, in the absence of net current flow (nil
current condition) the gradient of the electric potential
can be expressed as (Newman, 1991):

where t denotes the time. Inserting Eqn. (1) into (2)

leads the governing partial differential equation for where the superscripts indicate cations and anions,
diffusive mass transport: respectively, and c denotes the salt concentration. +
and are the stoichiometric coefficients of cations
and anions. Integration of Eqn. (8) between two arbi-
trary points leads to the diffusion potential (liquid
junction potential or Nernst potential) for a binary
electrolyte as:
In order to complete this system of equations, one has
to define a relationship that couples the transport of
all ionic species to one another. Commonly applied
constraint equations are:
The sign of the diffusion potential depends on the
integration domain, i.e., c1 and c2 , which in the case
of one dimensional diffusion are the upstream and
downstream concentrations of the salt, respectively.
Insertion of Eqn. (8) into the NP Eqn. (1) delivers the
steady state diffusion equation for binary electrolytes:

418
Copyright 2005 Taylor & Francis Group plc, London, UK
where ji are the molar flux densities of cations and where Dapp,i and app denote the homogenized (appar-
anions and DNP denotes the salt diffusion coefficient. ent) diffusion coefficient and permittivity of the REV.
Eqn. (10) indicates that a salt dissolved in a solu- X is the fixed volumetric surface charge density related
tion behaves like a single ionic species because of to particle surface charges and is the sign of the
the electro-neutrality requirement. However, at the surface charge (i.e., negative for negatively charged
molecular level, different self-diffusion coefficients of particles and positive for positively charged particles).
cations and anions result in separation of the species. For uncharged porous materials X equals zero. For
This leads to creation of a minute dipole density which this case the apparent diffusion coefficient can be
then prevents further separation. The dipole density estimated as (Bear and Bachmat, 1991):
creates a potential (the diffusion potential, Eqn. (9))
which acts to speed up the ion with the smaller self-
diffusion coefficient and slow down the ion with the
larger self-diffusion coefficient. with the porosity and the (purely geometrical) quan-
tity representing the tortuosity of the pore space of
the material.
2.2 Macroscopic diffusive transport: volume
averaging approach
3 NUMERICAL ANALYSES
Eqn. (1), (2), and (7) describe electro-diffusive trans-
port of ions at the microscale. However, most engi- In the following, two steady state diffusive trans-
neering applications deal with the investigation of port problems are considered: diffusion of a binary
macroscopic transport properties such as measure- electrolyte (NaCl) in a pure solution system (Sub-
ment of macroscopic diffusion coefficients, conduc- section 3.1) and in a porous medium contain-
tivities, and membrane potentials. There are various ing charged particles (Subsection 3.2). Parameters
alternatives to describe diffusive transport quantities used for the coupled PNP system of equations
of porous materials at the macroscale. Given the intrin- are given as: 0 = 8.85e12 C2 /(J m) , w = 78,
sic complexity of the pore structure of most geo- and D+ = 1.33e09 m2 /s, D = 2.03e 09 m2 /s, F =
biomaterials, modeling of the complete pore struc- 96500 C/mol, T = 293 K, and R = 8.31 J/(K mol).
ture constitutes a formidable task (Bear and Bachmat,
1991). A commonly applied procedure is the vol-
ume averaging method (Bear and Bachmat, 1991) 3.1 Diffusive transport in pure solution system
(Whitaker, 1999), where microscopic quantities are The first example deals with (one dimensional) diffu-
averaged over a representative volume element (RVE). sive flow of a binary electrolyte (NaCl) in a 2D domain.
The notion of the RVE forms the basis of the volume The length L of the unit cell is 150 nm. The electro-
averaging method (homogenization). The RVE must static boundary conditions (for Eqn. (7)) where chosen
be large enough in order to capture material properties such that zero charge density was applied at the left
(such as the porosity, aggregate fraction, and thickness and right boundary, whereas voltage boundary condi-
of double layers) correctly. In the following, a volume tions were specified at the upper and lower boundary
averaged quantity is defined as: (see Figure 1). Zero reference voltage was applied

Integration of the governing equations over the volume

of the RVE leads the macroscopic transport equations
(Samson et al., 1999):

Figure 1. Diffusive transport in pure solution system:

geometric dimensions and boundary conditions.

419
Copyright 2005 Taylor & Francis Group plc, London, UK
 Figure 4. Diffusive transport in charged porous media: peri-
odic microstructure and unit cell with spherical inclusion
(figures not to scale).

Independent of the ratio c1 /c2 a linear concentra-

tion distribution for anions and cations is obtained (see
Figure 2). On the other hand, the nonlinearity of the
voltage strongly depends on the applied concentration
Figure 2. Diffusive transport in pure solution system: con- ratio (Figure 3). For large values of c1 /c2 steep gradi-
centration distributions of anions and cations (at cross-section ents of voltage near the lower boundary are obtained
x = 0.75e 7 m).
(Figure 3). These strong nonlinearities near the lower
boundary lead to inaccurate numerical results using
coarse FE meshes. A numerical investigation using
various FE meshes which have been refined at the
lower boundary (i.e., where high nonlinearities are
expected) showed good agreement between numerical
and analytical results. Finally, we note that the salt dif-
fusion coefficient (see Eqn. (10)) is independent of the
ratio c1 /c2 and constant, i.e., Dsalt = 1.61 109 m2 /s.

3.2 Diffusive transport in charged porous media

The second example deals with the investigation of
diffusive transport through a charged porous medium
in a 2D domain. In the course of the periodic medium
approach the PDEs are applied to one particular unit
cell (see Figure 4). The size of the unit cell is 150 nm
Figure 3. Diffusive transport in pure solution system: volt- square, with a spherical inclusion of R = 15 nm. The
age distributions for various ratios of c1 /c2 (at cross-section boundary conditions for the concentration and voltage
x = 0.75e 7 m). are also shown in the figure (c1 > c2 , c2 = 0 mol/m3 ).
Negative surface charges s (in C/m2 ) ranging between
0 and 0.01 C/m2 are applied at the boundary of the
at the upper boundary condition, whereas at the inclusion. Zero reference voltage was applied at the
lower boundary was computed numerically (using a upper boundary, whereas 2 was computed from the
bisectioning algorithm) from the nil current condition nil current condition (Eqn. (6)).
(Eqn. (6)). The first investigation deals with the influence of
A triangular finite element mesh consisting of 772 the applied concentration gradient |(c2 c1 )/L| and
elements was used for the numerical simulations. surface charge density s on the ion flux density.
As shown in the previous section it turns out that for Figure 5 shows a strong increase of the ion flux
the special case of a binary electrolytic solution this density with increasing concentration gradient irre-
voltage can be expressed analytically (see Eqn.(9)). spective of the surface charge. The ion flux decreases
However, from the numerical studies it was found with increasing surface charge. However, for large
that dependent of the concentration ratio c1 /c2 and concentration gradients the ion fluxes of the charged
the chosen finite element mesh rather large deviations porous unit cell approach the ion flux obtained from
between numerical and analytical determined voltage the uncharged case.
may be obtained. The reason for the sensitivity of This behaviour can be explained best by looking
the numerical results on c1 /c2 can be explained best at the concentrations of cations and anions along a
by looking at the concentration distribution and the cross-section of the unit cell for a constant concentra-
voltage distribution (see Figures 2 and 3). tion gradient (|(c2 c1 )/L| = 1.0/L). Figure 6 clearly

420
Copyright 2005 Taylor & Francis Group plc, London, UK
 Figure 7. Diffusive transport in charged porous media: volt-
Figure 5. Diffusive transport in charged porous media: age along cross-section x = 0.75e 7 m for charged and
dependence of cation (anion) flux density on concentration uncharged particle.
gradient and surface charge.

x 10-3
12
=0
= -0.001
10
flux [mol/(m2 s)]

4
0.6 0.7 0.8 0.9 1
porosity

Figure 6. Diffusive transport in charged porous media: Figure 8. Diffusive transport in charged porous media:
concentration of cations (anions) along cross-section dependence of ion flux density on porosity for charged and
x = 0.75e 7 m for charged and uncharged particle. uncharged particle.

shows that for a surface charge s = 0.001 the cation keeping the surface charge constant leads to a com-
concentration at the upstream side of the particle is pression of the (diffuse) double layer and, hence, to an
approximately 1.10 mol/m3 compared to a concentra- increase in flux density (see Figure 5).
tion of 0.70 mol/m3 for the uncharged case. This is Finally, we investigate the influence of particle size
an almost 1.6 fold increase in concentration on the (i.e., porosity) on the ion flux density for a constant
particle surface. On the other hand, the anion concen- concentration gradient. Figure shows the dependence
tration (at the upstream side) is reduced 1.7 fold, from of ion flux on the porosity. Increasing the particle size
0.70 mol/m3 to 0.42 mol/m3 . This exclusion of anions (i.e., decreasing the porosity) leads to a rather strong
from the pore space is commonly referred to as anion reduction of ion flux. For the investigated case the
exclusion (see (Westermann-Clark and Christoforou, influence of surface charge on the ion flux is not very
1986) for details). pronounced.
The change of ion concentration around the charged
particle is accompanied by a strong increase of the
voltage near the particle surface (see Figure 7). The 4 CONCLUSIONS
length defined as the (normal) distance from the par-
ticle surface at which the voltage reduces to 1/e of the From the numerical results obtained from the cou-
surface voltage is generally denoted as double-layer pled PNP equations the following conclusions can be
thickness. Increasing the concentration gradient while drawn: It was found that the gradient of the electric

421
Copyright 2005 Taylor & Francis Group plc, London, UK
potential strongly depends on the applied gradient of Cussler, E. 1997. Diffusion mass transfer in fluid systems, 2nd
concentration. Large gradients in voltage may lead to Edition, Cambridge University Press, New York, USA.
numerical inaccuracies using coarse FE meshes. For Helfferich, F. 1962. Ion exchange, Series in Advanced Chem-
electro-diffusive ion transport through charged porous istry, McGraw-Hill, New York, USA.
Kato, M. 1995. Numerical analysis of the Nernst-Planck-
materials it was found that the fluxes of anions and Poisson system, Journal of Theoretical Biology 177,
cations strongly depend on the gradient of concentra- 299304.
tion and on the surface charge density on particles. MacGillivray, A. 1968. Nernst-Planck equations and the
For the case of large concentration gradients the electroneutrality and Donnan equilibrium assumptions,
influence of surface charge becomes less pronounced Journal of Chemical Physics 48, 29032907.
due to compression of diffuse double layers. Nega- Newman, J. 1991. Electrochemical Systems, 2nd Edition,
tively charged particles lead to anion exclusion and Prentice-Hall, Englewood Cliffs, NJ, USA.
cation inclusion. For the investigated range of sur- Samson, E., Marchand, J., Robert, J.-L., Bournazel, J.-P.
face charge a similar dependence of ion flux on the 1999. Modelling ion diffusion mechanisms in porous
media, International Journal for Numerical Methods in
particle size (porosity) as obtained for the uncharged Engineering 46, 20432060.
case was observed. Finally it is noted that the pro- Samson, E., Marchand, J., Beaudoin, J. (1999). Describing
posed model provides a powerful tool for investigation ion diffusion mechanisms in cement-based materials using
of coupled diffusion and migration of ions. In order the homogenization technique, Cement and Concrete
to further evaluate the model, comparison between Research 29, 13411345.
(macroscopic) experimental data such as effective dif- Smith, D., Pivonka, P., Jungnickel, C., Fityus, S. 2004. Theo-
fusion coefficients and membrane potentials and the retical analysis of anion exclusion and diffusive transport
numerical results will be made after upscaling. through a platy clay soil, Transport in Porous Media, 57,
251277.
Westermann-Clark, G. and Christoforou, C. (1986), The
exclusion-diffusion potential in charged porous mem-
REFERENCES branes, Journal of Electroanalytical Chemistry 198,
213231.
Auriault, J.-L.1991. Heterogeneous medium. Is an equiv- Whitaker, S. 1999. The method of volume averaging, Kluwer
alent macroscopic description possible?, International Academic Publishers, Dordrecht, The Netherlands.
Journal of Engineering Science 29 (7), 785795.
Bear, J. and Bachmat, Y. 1991. Introduction to modelling of
transport phenomena in porous media, Vol. 4, Kluwer
Academic Publishers, Dordrecht, The Netherlands.

422
Copyright 2005 Taylor & Francis Group plc, London, UK