Author Archives: Martin

What is a “diffusion path” anyway?

In a previous post I discussed how parts of the bentonite1 research community unjustifiably explain variation in (effective) diffusion rates as changes in “diffusion paths”. But what do authors really mean when using the term “diffusion path”?

In the geochemical/reactive transport literature, “diffusive pathways” are usually introduced when discussing the (presumably) related concept of tortuosity. For example, Steefel and Maher (2009) present a figure very similar to the one below, with the caption “Tortuous diffusion paths in porous material.”

Tortuous diffusion paths in porous material (Steefel and Maher (2009))

The text explains further that these are paths “the solute […] follow[s] in tortuous media”2. Several questions immediately arise. For example:

  • What is it, exactly, that follows these paths?
  • Why do the paths have a direction?
  • Why are these particular paths singled out? What stops e.g. the “Long Path” from taking this obvious shortcut:
Alternative path in Steefel and Maher's illustration

Let’s start with a hopefully obvious statement: Individual molecules or ions do not follow “paths” in a diffusive process, but conduct random motion. So paths as those in the figure above are certainly not trajectories of single particles.3

An answer to what is illustrated may be found in Van Brakel and Heertjes (1974) a main reference in the bentonite research community when discussing tortuosity etc.4 In this work, the system is analyzed in steady-state, and the following description is given for “diffusion paths”

Assume the pore space of the porous medium to be completely filled with what we will call diffusion paths. The main direction of the diffusion paths is the same as that of the concentration gradient. In the way the diffusion paths wind through the pore space they can be compared with the streamlines for laminar flow in porous media.

Here it may sound as if the authors reject Fickian diffusion (which is always parallel to the concentration gradient), but it is rather that they use the term “concentration gradient” for denoting the externally applied constant concentration difference, required to maintain a steady-state flux. To me, this is quite confusing, because the comparison of diffusion paths with streamlines implies that they are directed along a (negative) concentration gradient. Two “gradients” must thus be kept in mind simultaneously: the external concentration difference, and the actual gradient on the pore-scale.

But with this distinction made, it is clear what Van Brakel and Heertjes mean by “diffusion paths”, and that their aim is to reduce a more complex 3D problem to an effective 1D description. It also becomes relatively clear that this is the way the term is used in much of the bentonite literature. It explains e.g. why the “paths” in the above picture (and others) have a direction: they must be thought of as the steady-state flow on the pore scale, with an implied constant concentration difference applied across the sample, making the macroscopic flow effectively 1D.

This implied reduction to steady-state transport in 1D is often found in the literature, e.g. in Shackelford and Moore (2013)

This increased tortuosity reduces the macroscopic concentration gradient (i.e., increases the distance over which the concentration difference is applied) and, therefore, reduces the diffusive mass flux relative to that which would exist in the absence of the porous medium.

or in Altman et al. (2015)

\(\tau\) is a geometrical factor (\(\le 1\)) representing the reduction in the effective concentration gradient (d(Me)/dx) due to the fact that diffusion paths through a porous medium will generally be greater, i.e. more tortuous, than the straight-line distance between the system boundaries, i.e. dx.

I am quite puzzled by this description for several reasons. Firstly, I find it unsatisfying that these definitions require the system to be in steady-state. Information on the influence of geometry is of course contained in the diffusion coefficient itself, independent of any external concentration differences. To associate “tortuosity” with such concentration differences, rather than with the mobility of the diffusing substance, seems inadequate to me. Moreover, the procedure of reducing a “macroscopic” concentration gradient due to path length seems to only work for an isolated path. At least, the procedure must become more involved for a system of connected paths, something I’ve not seen commented on by authors adopting this concept.

Secondly, note that a “diffusion path” — with the steady-state definition — simply indicates net mass transfer of diffusing substance. The absence of a “diffusion path” in a region does not mean that dissolved constituents don’t migrate there, but only that flux contributions in different directions add up to zero net mass transfer. I did a silly random walk simulation to illustrate this point (the concentration of walkers is kept at a constant finite value on the left side of the domain, while it is kept at zero on the right side)

Silly random walk

Note that with the definitions here discussed, we must accept that the vertical section in this illustration is not a diffusion path. This situation is quite distinct from laminar advective transport, where — if I’m thinking correctly — the absence of a streamline implies the absence of motion.

Thirdly, if you consider a porous system to be a network of thin cylinders, I guess the steady-state flux will basically resemble the system itself (interconnected 1D-spaghetti configured in 3D). I suspect that this is the view many authors have in mind when speaking of “diffusion paths”. But, if so, why not simply speak of “paths”? Note also that the pore volume of smectitic systems mainly consists of 2D water films configured in 3D (it is lasagna, not spaghetti).

Lastly, what about non-steady-state transport? Concepts like the ones discussed here are also used when describing closed-cell diffusion tests , but are seldom (never?) defined in any other way. How could e.g. “tortuosity” reduce a macroscopically applied (1D) gradient in this case? And what is even meant by “diffusion paths”, if these are defined in steady-state? Since non-steady state is the general case, it would be more satisfying if quantities obtained under such circumstances were applied to the steady-state, rather than the other way around.

To get a feel for how pore geometry influences diffusion in non-steady-state, I conducted some more random walk simulations. In the animation below is compared random walks in an unrestricted 2D plane (blue) with random walks on a square net (red; strip width: 1 length unit, square size: 20 length units)

Random walks in 2D and on net

To quantify the diffusivity, we plot the average of the square of the displacements, \(\langle r^2 \rangle\), as a function of time5 in the two systems

<r2> vs t for 2D and isotropic net

We see that the diffusivity — which is directly proportional to the slope of these curves — is very close to twice as large in the unrestricted case as compared with diffusion on the net. From such a result it may be tempting to conclude that this reduction by a factor of two is due to longer “diffusion paths” on the net (and relate it to \(\sqrt{2}\), which conveniently is the ratio between the side and the diagonal of a square). But note that the diffusional process is isotropic also on the net, as demonstrated by identical slopes of the angle-resolved \(\langle r^2 \rangle\)-curves. Thus, interpreted in terms of “diffusion paths” on the net — however these should be defined — the conclusion is that the “paths” have the same length in any direction.

But the situation is easily analyzed from the simple model underlying the figures displayed above: in the 2D-plane, the random walk process has a maximized variance, because movements in the \(x-\) and \(y-\)directions are uncorrelated. The net geometry, on the other hand, correlates these variables: if a walker has free passage in the \(y\)-direction, it is restricted in the \(x\)-direction, and vice versa. Thus, the diffusivity is not diminished due to longer “paths”, but because the geometrical restrictions reduce the variance of the underlying process. This effect will depend on the relative reduction of dimensions: with line-like pores in a 3D configuration, the reduction factor becomes 1/3 (I guess this is what what is alluded to for a “homogeneous isotropic pore network” in the often-cited work Dykhuzien and Casey (1989) ), but for the case relevant for bentonite — diffusion in 2D-planes configured in 3D — the factor is 2/3 (which I haven’t seen stated anywhere). I furthermore don’t understand why such a factor should be termed “tortuosity”, because there is nothing intrinsically “tortuous” about it (in a sense one could even argue that individual trajectories in the unrestricted 2D-plane are more “tortuous” than the ones on the net).

By making the net non-isotropic, e.g. by replacing the squares by rectangles like this

Rectangluar net

the correlation between the \(x\)- and \(y\)-variables alters (it is now twice as likely for a walker to have no restriction in the \(y\)-direction as in the \(x\)-direction), which is directly reflected in the diffusivities

<r2> vs t for square and rectangular net

The diffusivity in the \(y\)-direction is now about twice as large as the diffusivity in the \(x\)-direction. Also the diffusivity in the diagonal directions is significantly larger than in the \(x\)-direction. Following a naive definition of “tortuosity”, this result may seem surprising (is the “solute” in the \(x\)-direction following a longer path than than in the diagonal directions?). Still, with a correct averaging procedure I guess the diffusivity can be related to “paths” on the net (However, I still don’t understand how to differ “diffusion paths” from geometrical paths).

With these simulations I simply want to argue for that it seems considerably more subtle and complex to relate pore geometry to diffusivity, than how it typically is presented in the bentonite literature. To be frank, I consider most talk about “diffusion paths” in the bentonite literature, as well as most definitions of various geometric factors, to be just that: talk. There is an established “tradition” to mention certain concepts (geometric factors, tortuosity, constrictivity, paths…), but in the end the introduced factors are usually only functioning as fudge factors, leading to unjustified claims about the nature of bentonite. Similarly, discussions on actual values of such factors are in principle always only qualitative.

As an example, Choi and Oscarson (1996) interpret different values of diffusivity measured in Na- and Ca-bentonite directly as a difference in “tortuosity”:

We attribute this to the larger quasicrystal, or particle, size of Ca- compared to Na-bentonite. Hence, Ca-bentonite has a greater proportion of relatively large pores; this was confirmed by Hg intrusion porosimetry. This means the diffusion pathways in Ca-bentonite are less tortuous than those in Na-bentonite.

But they could have been considerably more quantitative than this. In the paper, tortuosity is defined as \(\tau = L^2/L_e^2\), where “\(L\) is the straight-line macroscopic distance between two points defining the transport path, and \(L_e\) is the actual, microscopic or effective distance between the same two points.” Tortuosity is furthermore evaluated from HTO diffusion to \(\tau_{\ce{Na}} = 0.062\) in Na-bentonite, and \(\tau_{\ce{Ca}} = 0.117\) in Ca-bentonite. Combining these expressions gives \(L_{e,\ce{Na}} = 1.37\cdot L_{e,\ce{Ca}}\).

What is implicitly said in this work is thus that the “actual, microscopic or effective distance between two points” is 1.37 times longer in Na-bentonite as compared with Ca-bentonite. I mean that it would be suitable for authors making these kind of (implicit) claims to provide a quantitative idea of how the pore space is modified in order to achieve this particular alteration of distances. To me, it is not even obvious why “larger quasicrystals” implies shorter “diffusion paths” — note that the effect of the “net” geometries above are scale independent.

But rather than making a quantitative discussion, Choi and Oscarson (1996) give the following caveat

In reality, \(\tau\) may account for more than just the pore geometry of the clay. Another factor that may be included in \(\tau\) is, for instance, the variation in the viscosity of the solution within the pore space (Kemper et al., 1964).

I find this an incredible statement. It is similar to saying that, “in reality”, Earth’s gravity constant (\(g\)) may include effects of air resistance.

Footnotes

[1] “Bentonite” is used in the following as an abbreviation of “Bentonite and claystone”.

[2] That is at least my interpretation. A fuller quotation reads: “[Tortuosity] is defined as the ratio of the path length the solute would follow in water alone, \(L\), relative to the tortuous path length it would follow in porous media, \(L_e\)” (while the following equation actually contains the square of this ratio).

[3] I am not fully convinced that all authors keep this in mind at all times. How should e.g. the following passage from Charlet et al. (2017) be interpreted: “An important geometric parameter is the tortuosity factor, \(\tau\), that quantifies the travelled distance of a dissolved constituent through the pore network compared to actual distance between two points.”? Or this one from Van loon et al. (2018) : “The tortuosity is a measure for path lengthening and takes into account that molecules have to diffuse around grains and thus take a longer way.”?

[4] Which is quite amazing, considering that this paper deals with diffusion in a gas phase in macroporous systems.

[5] Since all walkers start in the same point (the origin) the data show finite-size effects for small times. The presented data is therefore taken after an initiation time, labeled \(t^\star\).

Evidence for anions in montmorillonite interlayers (swelling pressure, part II)

It is easy to find models assuming montmorillonite interlayers devoid of “anions” . Here I will present empirical evidence that such an assumption is incorrect. Before doing so, just a quick remark on the term “anions” in this context. If anions reside in interlayers, they certainly do so accompanied by excess cations, in order to maintain overall charge neutrality. Thus, when speaking of “anions” in the interlayer we really mean “salt” (= anion(s) + cation(s)). In the following I will use the term “salt” because it better reflects the overall charge neutral character of the process (we are not interested in pushing a handful of negative charge into an interlayer).

The nature of bentonite swelling

The evidence for salt having access to interlayers follows directly from the observed swelling pressure response to changes in external salinity. It is therefore important to first understand the thermodynamic basis for swelling pressure, which I wrote about in an earlier post (the same nomenclature is adopted here). In essence, swelling is a consequence of balancing the water chemical potential1 in the clay with that in the external solution2, and swelling pressure quantifies the difference in chemical potential between the external solution and the non-pressurized bentonite sample, as illustrated here

chemical potentials in non-pressurizied bentoniote and in external solution

Since the chemical potential in the external solution depends on the salt content, we generally expect a response in swelling pressure when altering external salinity.

Labeling the salt concentration \(c^{ext}\), we write the chemical potential of the external solution in terms of an osmotic pressure3

\begin{equation}
\mu_w^{ext} = \mu_0 – P_{osm.}^{ext}(c^{ext})\cdot v_w
\tag{1}
\end{equation}

where \(v_w\) is the partial molar volume of water. \(P_{osm}^{ext}\) is not the pressure in the external solution, but the pressure that would be required to keep the solution in equilibrium with pure water. The actual pressure in this compartment is the same as for the reference state: \(P_0\). It may seem confusing to use a “pressure” to specify the chemical potential, but we will see that it has its benefits. Experimentally we have full control of \(P_{osm}^{ext}\) by choosing an appropriate \(c^{ext}\).

Response in an indifferent clay

With salt in the external solution, the big question is what happens to the chemical potential of the clay. We will start by assuming (incorrectly) that external salt cannot access the interlayers. This means that the chemical potential of the (non-pressurized) bentonite does not change when the external salinity changes. We refer to this hypothetical bentonite as indifferent. In analogy with the external solution, we write the chemical potential of the indifferent non-pressurized bentonite as4 \begin{equation} \mu_w^{int}(P_0) = \mu_0 -P_s^0\cdot v_w \;\; \;\; \;\; \text{(indifferent clay)} \end{equation}

were \(P_s^0\) is the swelling pressure in case of pure water as external solution. By assumption, \(\mu_w^{int}(P_0)\) does not depend on the external salinity (it is independent of \(P_{osm}^{ext}\)). The chemical potential in the indifferent clay at an elevated pressure \(P\) is

\begin{equation}
\mu_w^{int}(P) = \mu_0 – P_s^0\cdot v_w +(P-P_0)\cdot v_w \;\;
\;\; \;\; \text{(indifferent clay)}
\tag{2}
\end{equation}

The swelling pressure (defined as the difference in pressure between bentonite and external solution, when the two are in equilibrium: \(P_s \equiv P_{eq} – P_0\)) in an indifferent clay is given by equating eqs. 1 and 2, giving the neat formula

\begin{equation}
P_s(c^{ext}) = P_s^0 – P_{osm}^{ext}(c^{ext}) \;\; \;\; \;\;
\text{(indifferent clay)}
\end{equation}

Note the following:

  • Although an indifferent clay is not affected by salt, it certainly has a swelling pressure response, demonstrating that swelling pressure depends as much on the external solution as it does on the clay.
  • Since swelling pressure in this case decreases linearly with the osmotic pressure of the external solution, it is predicted to vanish when the osmotic pressure equals \(P_s^0\).
  • External osmotic pressures larger than \(P_s^0\) implies “drying” of the clay (water transport from the clay into the external compartment)

If the above derivation feels a bit messy, with all the different types of pressure quantities to keep track of, here is a hopefully helpful animation

Animation swelling pressure response without anions in interlayer

Real swelling pressure response

Equipped with the swelling pressure response of an indifferent clay, let’s compare with the real response: The swelling pressure response in real bentonite deviates strongly from the indifferent clay response. This is seen e.g. here for Na-montmorillonite equilibrated in sequence with NaCl solutions of increasing concentration5 (data from Karnland et al., 2005 )

Swelling pressure response to salinities mid range densities

Swelling pressure indeed drops with increased concentration, but the drop is not linear in \(P^{ext}_{osm}\), and is weaker as compared with the indifferent clay response (shown by dashed lines). All samples in the diagram above exert swelling pressure when \(P^{ext}_{osm} \gg P_s^0\), i.e. far beyond the point where the swelling pressure in an indifferent clay is lost.

The deviation of the observed response from that of an indifferent clay directly demonstrates that the clay is affected by salt, i.e. that the chemical potential of the non-pressurized clay depends on the external salt concentration. The only reasonable way for salt to influence the chemical potential in the bentonite is of course to reside in the interlayer pores. Consequently, the observed swelling pressure response proves that salt from the external solution enters the interlayer pores.

Here is an illustration of how the chemical potentials relate to the swelling pressure in real bentonite

Swelling pressure repsonse real bentonite

Although the observed swelling pressure response in itself is sufficient to dismiss the idea that salt does not have access to interlayers, the study by Karnland et al., (2005) provides a much broader verification of the thermodynamic nature of swelling pressure. In particular, the chemical potential was measured (by means of vapor pressure) separately in the same samples as used for swelling pressure tests, after they had been isolated and unloaded. The terms in the relation \(P_s = \left(\mu_w^{ext} – \mu_w^{int}(P_0) \right)/v_w\) were thus checked independently, as indicated here

Measurements performed in Karnland et al. (2005)

A striking confirmation of salt residing in interlayers is given by the observation that the chemical potential in the non-pressurized samples is lower than that in the corresponding external solution, as well as that in non-pressurized samples of similar density, but equilibrated with pure water.

Another interesting observation is that the sample with the highest density behaves qualitatively similar to the others: although the external osmotic pressure never exceeded \(P_s^0\) (\(\approx\)56 MPa), the response strongly deviates from that of an indifferent clay

swelling pressure response to salinity high density

Because the pore space of samples this dense (\(2.02\;\mathrm{g/cm^3}\)) mainly consists of mono- and bihydrated interlayers, this similarity in response shows that salt has access also to such pores.

Implications

The issue of whether “anions” have access to montmorillonite interlayers has — for some reason — been a “hot” topic within the bentonite research community for a long time, and a majority of contemporary models rest on some version of the assumption that “anions” does not have access to the full pore volume. But, as far as I can see, this whole idea is based on misconceptions. I guess that saying so may sound quite grandiose, but note that swelling pressure is not at all considered in most chemical models of bentonite. And if it is, it is usually treated incorrectly. As an example, here is what Bradbury and Baeyens (2003) writes in a very influential publication

One of the main premises in the approach proposed here is that highly compacted bentonite can function as an efficient semi-permeable membrane (Horseman et al., 1996). This implies that the re-saturation of compacted bentonite involves predominantly the movement of water molecules and not solute molecules. Thus, to a first approximation, the composition of the external saturating aqueous phase should be a second-order effect which has little influence on the initial compacted bentonite porewater composition.

If the composition of the re-saturating water were to play an important role in determining the porewater composition, then it should also have a significant influence on swelling (Bolt, 1979). Dixon (2000) recently reviewed the role of salinity on the development of swelling pressure in bentonite buffer and backfill materials. He concluded that provided the initial dry densities were greater than 900 \(\mathrm{kg\;m^{-3}}\), the swelling pressures developed are unaffected for groundwater salinities \(< 75 \;\mathrm{g\;l^{-1}}\). Even brines appear to have little or no influence for initial dry densities \(>1500 \;\mathrm{kg\;m^{-3}}\).

But, as we just have learned, a system with a weak swelling pressure response necessarily has a significant contribution to its water chemical potential due to externally provided salt. In contrast, the approximation discussed in the first paragraph of the quotation — which is basically that of an indifferent clay — maximizes the swelling pressure response. Thus, the discussed “main premise” does not hold, and the provided empirical “support” is actually an argument for the opposite (i.e. that salt has access to the clay).

Footnotes

[1] In the following I will write only “chemical potential” — it is always the chemical potential of water that is referred to.

[2] This is just a complicated way of saying that swelling is (an effect of) osmosis.

[3] Some may say that \(P_{osm}^{ext}\) is simply the “suction” of the solution, but I am not a fan of using that concept in this context. I will comment on “suction” in a later blog post. Update 210816: “suction” is discussed here.

[4] The density dependence of the chemical potential in the bentonite is not explicitly stated here, in order to keep the formulas readable, but we assume throughout that the bentonite has some specific water-to-solid mass ratio \(w\).

[5] The NaCl concentrations are 0.0 M, 0.1 M, 0.3 M, 1.0 M, and 3.0 M.

The problem with geometric factors

Many papers on diffusion in compacted bentonite and claystone1 assume a direct connection between diffusivities measured in clay (\(D_\mathrm{clay}\)), and the geometrical configuration of the pore space. It seems close to mandatory to include
an
expression of this type \begin{equation} D_\mathrm{clay} = \frac{1}{G}\cdot D_\mathrm{w}\tag{1} \end{equation}

where \(D_\mathrm{w}\) is the corresponding diffusion coefficient in bulk water (usually taken as the limiting value at infinite dilution), and \(G\) is a factor postulated to have pure geometric meaning. The geometric multiplier is written in a variety of ways, sometimes in terms of other factors which individually are supposed to convey specific geometric information, e.g. \(1/G = \delta/\tau^2\), where \(\delta\) is the “constrictivity” and \(\tau^2\) is the “tortuosity factor”2.

However, the above described procedure is futile:

  1. A geometric factor can not be deduced a priori, because the pore space geometry of any realistic system is way too complicated. So-called tortuosities have only been calculated in some trivial cases, relating e.g. to highly symmetric networks of identical cylinders, which has basically no resemblance to the chaotic smectite-lasagna expected to be found in bentonite. Indeed, many papers on these issues admit that the geometric factor in practice functions as a “fudge factor” (to be “optimized”).
  2. Even if a geometric factor could be deduced from only considering the pore space configuration, it is not expected to relate to \(D_\mathrm{clay}\) in accordance with eq. 1. If we imagine removing all influence of pore geometry on the diffusivity in the clay, what remains is diffusion in an environment heavily influenced by the presence of exchangeable cations. Thus, a correct geometric factor should not multiply \(D_\mathrm{w}\), but the diffusion coefficient measured on the short time and length scale in the clay. Such measurements exist — at least for water diffusion — and show values different from \(D_\mathrm{w}\).

Because it is practically impossible to independently deduce a geometric factor, it is de facto defined by the ratio \(D_\mathrm{clay}/D_\mathrm{w}\) in many studies. I dare to say that in any work within the bentonite research community, a geometric factor has no other meaning than \(D_\mathrm{clay}/D_\mathrm{w}\) (or some function of this parameter, depending on definitions).

This means — at best — that the use of a geometric factor does nothing for a model description, but only replaces the parameter \(D_\mathrm{clay}\) with the parameter \(D_\mathrm{clay}/D_\mathrm{w}\); it simply corresponds to the awkward choice of specifying diffusion coefficients in units of \(D_\mathrm{w}\).

Often, however, the procedure has worse consequences, because eq. 1 (as used in practice) represents an unjustified claim; it is not correct to use the value of \(D_\mathrm{clay}/D_\mathrm{w}\) to make specific statements about the pore space geometry, while ignoring other possible mechanisms contributing to the value of \(D_\mathrm{clay}\). If diffusion instead is described without invoking a geometric factor, \(D_\mathrm{clay}\) itself functions as as a model parameter, and no implicit claims are made about why it has the value it has.

There are many examples in the literature where this confusion shows up. In e.g. Appelo and Wersin (2007), it is concluded that tortuosity for “anions” is larger than for water, while tortuosity for “cations” is smaller than for water. This conclusion is based on fitting of geometric factors without independent information on the pore space geometry:

[…] the pore may become filled entirely with a diffuse double layer when it narrows sufficiently. This constricts the passage of anions, and since the anions must circumnavigate the obstacle, they have greater tortuosity than tritium. This explains that a model with a tortuosity factor for iodide that is 1.4 times higher than for tritium better matches the data.

[…] in interlayers and pore constrictions the cations pass in relatively larger amounts than in free porewater, and consequently, they have smaller tortuosity than tritium.

Similar descriptions are given in Altmann et al. (2012) (in a section describing “advances” in process understanding with the FUNMIG project)

At low ionic strengths, cations are mostly located in the diffuse layer. This layer is of relatively low volume and connected throughout the porosity — cations need less time to explore this small volume and therefore to pass from one side of the media to the other.[…]

Anions are excluded from the diffuse layer and, therefore,also have a smaller volume to explore than neutral tracers. They however must pass through small pores in order to explore the total porosity, and these pores act as electrostatic throats […], forming an energy barrier to anion movement, reducing the probability of passage (resulting in reduced \(D_e\)) and increasing the time needed to explore all of the porosity.

And the same type of argument — that different diffusive behavior of cations and anions is caused by differences in “pathways” — is found in Charlet et al. (2017)

The different pathways of anions and cations in these regions lead to different effective diffusion rates (Van Loon et al., 2007; Mazurek et al., 2009)

But effective diffusion rates — which can basically be made to vanish e.g. for chloride, as has been demonstrated in bentonite — do not vary because of changes in “pathways”, but because of changes in boundary conditions, due to the presence of interfaces to external solutions.

To investigate the behavior of \(D_\mathrm{clay}/D_\mathrm{w}\) one must of course consider real — not effective — diffusion coefficients (unfortunately, actual diffusion coefficients are usually referred to as “apparent” in the bentonite research community). Real diffusion coefficients reflect the mobility of the diffusing species, and can be measured e.g. in non-steady-state closed-cell tests, which are not influenced by interfaces to external solutions3. A very good appreciation for how such diffusion coefficients behave is given by the vast amount of data measured by Kozaki et al. in Na-montmorillonite (Kunipia-F). Some of this data (sources: 1, 2, 3, 4, 5) is plotted below, in terms of \(D_\mathrm{clay}/D_\mathrm{w}\) (using a linear \(y\)-axis, for a change).

Kozaki's diffusion data.

This plot clearly shows that if \(D_\mathrm{clay}/D_\mathrm{w}\) is interpreted as having pure geometrical meaning, the conclusion must be that chloride and sodium have basically identical “tortuosity”, while the “tortuosity” of the cations calcium and cesium is considerably larger.4 Thus, the closed-cell data disproves any claim that “tortuosity” for anions is larger than for cations. Moreover, from this data it is quite a stretch to state that “tortuosity” for anions is larger than for water, and it is plain wrong to state that “tortuosity” for cations is smaller than for for water. (For the record, I’m strongly convinced that both calcium and cesium occupy the same pore volume as sodium in this material.)

I should emphasize that I am not primarily arguing for a better geometrical interpretation of \(D_\mathrm{clay}/D_\mathrm{w}\), but rather that this whole procedure (using eq. 1) should be avoided. This is not to say that pore geometry doesn’t influence clay diffusivity — most reasonably it does — but analysis based on eq. 1 will be flawed, because of the unjustified model assumptions.

When it comes to geometrical understanding of diffusion in bentonite, I think it would be more rewarding to focus on the quite large variation of reported diffusivities in seemingly identical systems. As an example, the figure below shows the diffusivity of cesium in very similar, pure Na-montmorillonite, measured in four different closed-cell studies: apart from Kozaki et al.’s data (blue), also Sato et al. (1992) (green), and Tachi and Yotsuji (2014) (orange).

Cs diffusion data from Kozaki et al.  Sato et al., and Tachi and Yotsuji

In aggregation, this figure demonstrates almost an order of magnitude variability for the diffusion coefficient of cesium, even in very pure systems. If this kind of variability cannot be explained, it is of course not very relevant to develop models describing (smaller) differences in diffusivity between different types of ions in different types of bentonites.

If you claim to understand the variation seen in this plot, rather than simply finding an appropriate factor to blame (particle size?), here is a challenge: Give a prescription for how to prepare a Kunipia-F/P sample in order to minimize/maximize the resulting diffusivity.5

Footnotes

[1] In the following, the term “bentonite” is used as a shorthand for “compacted bentonite and claystone”.

[2] The nomenclature here is a mess, but I leave it to people who really want to use these concepts to actually define them in detail.

[3] They can also be deduced from the transient behavior in “steady-state” tests, if the transport capacity of the confining filters is not limiting the process.

[4] Note that “tortuosity”, as used e.g. by Appelo and Wersin (2007), is related to the inverse of \(D_\mathrm{clay}/D_\mathrm{w}\).

[5] The diffusion direction should still be along the axis of sample compaction, of course.

Swelling pressure, part I

I am puzzled by how bentonite swelling pressure is presented in present day academic works.

Here, I would like to revisit the pure thermodynamic description of swelling pressure, which I think may help in resolving several misconceptions about swelling pressure.

Of course, thermodynamics cannot answer what the microscopic mechanism of swelling is, but puts focus on other — often relevant — aspects of the phenomenon. We thus take as input that, at the same pressure and temperature, the water chemical potential2 is lowered in compacted bentonite as compared with pure water, and we ignore the (microscopic) reason for why this is the case. We write the chemical potential in non-pressurized3 bentonite as \begin{equation} \mu_w(w,P_0) = \mu_0 + \Delta \mu(w,P_0) \end{equation}

where \(\mu_0\) is a reference potential of pure bulk water at pressure \(P_0\) (isothermal conditions are assumed, and temperature will be left out of this discussion), and \(w\) is the water-to-solid mass ratio. Note that \(\Delta \mu(w,P_0)\) is a negative quantity.

The chemical potential in a pressurized system is given by integrating \(d\mu_w = v_wdP\), where \(v_w\) is the partial molar volume of water, giving4 \begin{equation} \mu_w(w,P) = \mu_0 + \Delta \mu(w,P_0) + v_w\cdot (P-P_0) \end{equation}

In order to define swelling pressure, we require that the bentonite is confined to a certain volume while still having access to externally supplied water, i.e. that it is separated from an external water source by a semi-permeable component. This may sound abstract, but is in fact how any type of swelling pressure test is set up: water is supplied to the sample via e.g. sintered metal filters.

With this boundary condition, a relation between swelling pressure and the chemical potential is easily obtained by invoking the condition that, at equilibrium, the chemical potential is the same everywhere. Assuming an external reservoir of pure water at pressure \(P_0\), its chemical potential is \(\mu_0\), and the equilibrium condition reads \begin{equation} \mu_w(w,P_{eq}) = \mu_0 + \Delta \mu(w,P_0) + v_w\cdot (P_{eq}-P_0) = \mu_0 \end{equation}

where \(P_{eq}\) is the pressure in the bentonite at thermodynamic equilibrium.

Defining the swelling pressure as \(P_s = P_{eq}-P_0\) we get the desired relation5 \begin{equation} P_s = -\frac{\Delta \mu(w,P_0)}{v_w} \tag{4} \end{equation}

Alternatively this relation can be expressed in terms of activity (related to the chemical potential as \(\mu = \mu_0 +RT\ln a\)) \begin{equation} P_s = -\frac{RT}{v_w}\ln a (w,P_0) \tag{5} \end{equation}

or, if the activity is expressed in terms of the vapor pressure, \(P_v\), in equilibrium with the sample, \begin{equation} P_s = -\frac{RT}{v_w}\ln \frac{P_v}{P_{v0}} \tag{6} \end{equation}

where \(P_{v0}\) is the corresponding vapor pressure of pure bulk water.

The above relation has been presented in the literature for a long time. But, as far as I am aware, direct interpretation of experimental data using eq. 4 is more scarce. Spostio (72) compares swelling pressures in Na-montmorillonite (reported by Warkentin et al 57) with water activities measured in the materials (reported by Klute and Richards 62) and concludes a “quite satisfactory” agreement of eq. 4 (the highest pressures were on the order of 1 MPa). He moreover comments

Future measurements of \(P_S\) and \(\Delta \mu_w\) for pure clays and soils as a function of water content would do much to help assess the merit of equation (11) [eq. 4 here].

Such “future” measurements were indeed presented by Bucher et al (1989), for “natural” bentonites in a density range including very high pressures (\(\sim 40\) Mpa). For “MX-80” the data looks like this

Here the value of \(v_w\) was set equal to the molar volume of bulk water when applying eq. 6. It is interesting to note that this value, which is necessarily correct in the limit of low density, appears to be valid for densities as large as \(2\;\mathrm{g/cm^3}\).

The clearest demonstration of the validity of eq. 4 is in my opinion the study by Karnland et al. (2005), where swelling pressure and vapor pressure were measured on the same samples. The result for Na-montmorillonite is shown below (again, the value of bulk water molar volume was used for \(v_w\)).

The above plots make it clear that the description underlying eq. 4 (or eq. 5, or eq. 6) is valid for bentonite, at any density. An important consequence of this insight — and something I think is often not emphasized enough — is that swelling pressure depends as much on the external solution as it does on the bentonite.

Measuring the response in swelling pressure to changes in the external solution is therefore a powerful method for exploring the physico-chemical behavior of bentonite. I will return to this point in later blog posts, in particular when discussing the “controversial” issue whether “anions” have access to montmorillonite interlayers.

The animation below summarizes the thermodynamic view of the development of swelling pressure: the external reservoir fixes the value of the water chemical potential, and in order for the bentonite sample to attain this level, its pressure increases.

Footnotes

[1] You can even find a statement saying that clay swelling has been proved to be “due to long-range interaction between particle surfaces and the water” (I don’t agree).

[2] In the following I will simply write “chemical potential”. Here the water chemical potential is the only one involved.

[3] Here “non-pressurized” means being at the reference pressure \(P_0\). In practice \(P_0\) is usually atmospheric absolute pressure.

[4] Here it is assumed that \(v_w\) is independent of pressure. Also, using \(w\) as thermodynamic variable implies that the water chemical potential is measured in units of energy per mass, which requires this volume factor to be the partial specific volume of water. Here we assume that the chemical potential is measured in units energy per mol, but use \(w\) for quantifying the amount of water in the clay, since it is the more commonly used variable in the bentonite world. The amount of moles of water is of course in strict one-to-one correspondence with the water mass.

[5] What is said here is that swelling pressure generally is identified as an osmotic pressure. I will expand on this in a future blog post.

You cannot add fluxes from different domains in multiporosity models

Most models of compacted bentonite and claystone1 assume that the material contains different domains with distinctly different properties. In particular, it is common to assume that the system contains both an aqueous domain (a bulk water solution) and one or several domains associated with the solid surfaces.2 Models of this type is here called multiporosity models.

A “standard” procedure in multiporosity models is to assume that the diffusive flux in any domain can be fully expressed using quantities that relates only to that same domain, while a “total” flux is given by summing the contributions from each domain. In e.g. the case of having one bulk water domain and one surface domain, the diffusive flux is typically written something like

$$ \label{eq1} j = j_{bulk} + j_{surface} = -D_e\cdot\nabla c -D_s \cdot\nabla c_s \tag{1} $$

were \(D_e\) and \(c\) denote, respectively, the diffusion coefficient and concentration in the bulk water domain, and the corresponding parameters in the surface domain are denoted \(D_s\) and \(c_s\). It is important to keep in mind that what is referred to as a “total” flux in these descriptions is the experimentally accessible, macroscopic flux (here simply labelled \(j\)).

This approach is basically adopted by all proponents of multiporosity models, and has been so for a long time, see e.g.:

There is, however, a big problem with this procedure:

Equation 1 does not make sense

The statement of eq. 1 — that the macroscopic, measurable flux has two independent contributions — is only true if the different domains are completely isolated from each other. This can e.g. be seen by making \(D_s\) vanishingly small in eq. 1; the flux is still non-zero due to the contribution from the completely unaffected \(j_{bulk}\)-term. In contrast, for non-isolated domains the expected behavior is that a macroscopic diffusive flux will vanish if the diffusivity of any domain is made to approach zero.

Consequently, a large part of the compacted bentonite research field seems to treat diffusion by assuming the different domains to be isolated “pipes”. But this is of course not how bentonite works! If you postulate multiple diffusive domains in the clay, you must certainly allow for diffusing species to be able to be transferred between them. And here is where it gets mysterious: the default choice is to also assume that the domains are non-isolated. In fact, the domains are typically assumed to be in equilibrium. In e.g. the case of having one bulk water and one surface domain, the equilibrium is usually expressed something like

$$
c_s = K\cdot c
\tag{2}
$$

where \(K\) is some sort of partition coefficient.

Thus, the “standard” treatment is to combine eq. 1 — even though it is only valid for isolated domains — with an equilibrium assumption (eq. 2) to arrive at expressions which typically look something like this

$$
j = \left (D_e + K\cdot D_s\right)\cdot \nabla c
\tag{3}
$$

Note how this equation has the appearance of an addition of diffusion coefficients, as the “total” flux now is related to a single concentration gradient (which in practice always is a bulk solution concentration).

Since eq. 3 results from combining two incompatible assumptions (isolation and equilibrium), it should come as no surprise that it is not a valid description for multiporosity models.

I cannot get my head around why this erroneous description has prevailed within the compacted clay research field for so long. The problem of deriving “effective” properties of heterogeneous systems appears in many research fields, and applies to many other quantities than diffusivity, e.g. heat conductivity, electrical conductance, elastic moduli, etc. It should therefore be clear also in the compacted clay world that obtaining the effective diffusivity for a system as a whole is considerably more involved than simply adding diffusivities from each contributing domain.

Some illustrations

To discuss the issue further, let’s make some illustrations. A molecule or ion in a porous system conducts random “jiggling” motion, exploring all regions to which it has access, as illustrated in this figure

This animation shows a simulated random walk at the microscopic scale, in a “made-up” generic porous system, consisting of solids (brown), a bulk water domain (blue), and a surface domain (pink). The green dots show the trace of the random walker; light green when the walker is located in the bulk water domain (with higher mobility), and darker green when it is located in the surface domain (with lower mobility). Note how the walker frequently switches between the domains — which is exactly what is required for the two domains to be in equilibrium.

As diffusivity is fundamentally related to the statistical properties of the random walk,3 it should be crystal clear that the macroscopic diffusion coefficient (reflecting mobility on length scales beyond what is shown in the animation) depends in a complicated way on the mobility of the involved domains, just as it depends on the detailed geometrical configuration. Oddly, the complex geometrical dependency is usually acknowledged within the bentonite research field (by means of e.g. formation and/or tortuosity factors), which makes it even more incomprehensible that the issue of multiple domain diffusivities is treated as it is.

In contrast, the figure below illustrates the microscopic behavior implied by eq. 1.

Here is shown one random walker in the bulk domain (green trace) and one in the surface domain (red trace). Since no exchange occurs between the domains, they are not in equilibrium (the concept does not apply). The lack of exchange between the domains is of course completely non-physical, but it is only under this assumption that the system is described by two independent fluxes (eq. 1). It may also be noted that the two types of walkers experience completely different geometrical configurations (i.e. formation factors), none of which corresponds to the geometrical configuration in the case of non-isolated domains.

A more reasonable formulation

If you insist on having a multiporosity model, the correct way to go about is to assume equilibrium within a so called representative elementary volume (REV), and to formulate the model on the macroscopic — averaged — level. On this level of description, it is legitimate to speak of both a bulk water and a surface concentration at any given point in the model (they represent REV averages), but there is of course only a single flux, as well as a single value of the chemical potential at any given point. The diffusive flux can be written quite generally as

\begin{equation} j = -\phi
D_\mathrm{macr.}\cdot\frac{\bar{c}}{RT}\nabla\mu
\end{equation}

where \(\bar{c}\) is the macroscopic, averaged, concentration, \(\phi\) is the porosity, and \(\mu\) is the (electro-) chemical potential. \(D_\mathrm{macr.}\) is the macroscopic diffusion coefficient, containing statistical information on how the diffusing substance traverses the porous system on length scales larger than the REV (in particular, it contains averaged information on the mobility in all involved domains, as well as their detailed geometrical configuration).

The value of \(\mu\) can be obtained from considering the microscopic scale. Since the system by assumption is in equilibrium within a REV, \(\mu\) has the same value everywhere within this volume. On the other hand, the (electro-) chemical potential may also be expressed e.g. as

\begin{equation}
\mu(\xi ) = \mu_0 + RT\cdot\ln a(\xi ) + zF\psi(\xi ) = \mathrm{const}
\end{equation}

where \(\xi\) is a microscopic variable denoting the position within the REV, \(\mu_0\) is a reference potential, \(a(\xi )\) is the chemical activity and \(\psi(\xi )\) the electrostatic potential. In the bulk water domain, the electrostatic potential is zero by definition (ignoring possible couplings between different types of charged species), and \(\mu\) can here be evaluated as (where the index “b” is used to denote bulk)

\begin{equation}
\mu_b = \mu_0 + RT\cdot\ln a_b
\tag{4}
\end{equation}

But this is by definition the value of the (electro-)chemical potential everywhere within the REV.

Thus, going back to the macroscopic scale and using eq. 4 for \(\mu\), we have for the flux

\begin{equation}
j = -\phi D_\mathrm{macr.}\cdot\frac{\bar{c}}{a_b}\nabla a_b
\end{equation}

Ignoring also possible gradients in activity coefficients, the flux can be simplified as

\begin{equation}
j = -\phi D_\mathrm{macr.}\cdot\frac{\bar{c}}{c}\nabla c
\end{equation}

Where, as before, \(c\) denotes a bulk solution concentration. For the case of having one bulk domain, and one surface domain, \(\bar{c}\) can be written

\begin{equation}
\bar{c} = \frac{\phi_b}{\phi}c + \frac{\phi_s}{\phi}c_s
\end{equation}

where \(\phi_b\) and \(\phi_s\) are the porosities for the bulk and surface domains, respectively (\(\phi = \phi_b + \phi_s\)). Using also eq. 2, the flux can finally be written

\begin{equation}
j = -D_\mathrm{macr.}\left (\phi_b+\phi_s\cdot K\right )\nabla c
\tag{5}
\end{equation}

This expression can be compared to the result from the “standard” procedure, eq. 3. Although eq. 3 is derived using incompatible, non-physical assumptions, it is similar to eq. 5 , but involves two (\(D_e\) and \(D_s\)) rather than a single diffusion parameter (\(D_\mathrm{macr.}\)). If these two parameters are used merely as fitting parameters (which they basically always are), eq. 3 can always be made to resemble eq. 5. In this sense, eq. 3 is not wrong. Rather, it is not even wrong.


Footnotes

[1] In the following I will only say bentonite, but I mean any type of clay system with significant swelling/ion exchange properties.

[2] I strongly oppose the way bulk water phases are included in such models, but this post is not a critique of that model choice. Here, I instead point out that virtually all published descriptions of multiporosity models handle the flux incorrectly.

[3] Specifically, the average of the square of the displacement of the walker is proportional to time, with the diffusion coefficient as proportionality constant (apart from a dimensional factor).