Category Archives: Multi-porosity models

Extracting anion equilibrium concentrations from through-diffusion tests

Recently, we discussed reported equilibrium chloride concentrations in sodium dominated bentonite, and identified a need to assess the individual studies. As most data is obtained from through-diffusion experiments, we here take a general look at how anion equilibrium is a part of the through-diffusion set-up, and how we can use reported model parameters to extract the experimentally accessible equilibrium concentrations.

We define the experimentally accessible concentration of a chemical species in a bentonite sample as

\begin{equation} \bar{c} = \frac{n}{m_\mathrm{w}} \end{equation}

where $n$ is the total amount of the species,¹ and $m_{w}$ is the total water mass in the clay.² It should be clear that $\bar{c}$, which we will refer to as the clay concentration, is accessible without relying on any particular model concept.

An equilibrium concentration is defined as the corresponding clay concentration (i.e. $\bar{c}$) of a species when the clay is in equilibrium with an external solution with species concentration $c^\mathrm{ext}$. A convenient way to express this equilibrium is in terms of the ratio $\bar{c}/c^\mathrm{ext}$.

The through-diffusion set-up

A through-diffusion set-up consists of a (bentonite) sample sandwiched between a source and a target reservoir, as illustrated schematically here (for some arbitrary time):

The sample length is labeled $L$, and we assume the sample to be initially empty of the diffusing species. A test is started by adding a suitable amount of the diffusing species to the source reservoir. Diffusion through the bentonite is thereafter monitored by recording the concentration evolution in the target reservoir,³ giving an estimation of the flux out of the sample ($j^\mathrm{out}$). The clay concentration for anions is typically lower than the corresponding concentration in the source reservoir.

Although a through-diffusion test is not in full equilibrium (by definition), local equilibrium prevails between clay and external solution⁴ at the interface to the source reservoir ($x=0$). Thus, even if the source concentration varies, we expect the ratio $\bar{c}(0)/c^\mathrm{source}$ to stay constant during the course of the test.⁵

The effective porosity diffusion model

Our primary goal is to extract the concentration ratio $\bar{c}(0)/c^\mathrm{source}$ from reported through-diffusion parameters. These parameters are in many anion studies specific to the “effective porosity” model, rather than being accessible directly from the experiments. We therefore need to examine this particular model.

The effective porosity model divides the pore space into a bulk water domain and a domain that is assumed inaccessible to anions. The porosity of the bulk water domain is often referred to as the “effective” or the “anion-accessible” porosity, and here we label it $\epsilon_\mathrm{eff}$.

Anions are assumed to diffuse in the bulk water domain according to Fick’s first law

\begin{equation} \label{eq:Fick1_eff} j = -\epsilon_\mathrm{eff} \cdot D_p \cdot \nabla c^\mathrm{bulk} \tag{1} \end{equation}

where $D_p$ is the pore diffusivity in the bulk water phase. This relation is alternatively expressed as $j = -D_e \cdot \nabla c^\mathrm{bulk}$, which defines the effective diffusivity $D_e = \epsilon_\mathrm{eff} \cdot D_p$.

Diffusion is assumed to be the only mechanism altering the concentration, leading to Fick’s second law

\begin{equation} \label{eq:Fick2_eff} \frac{\partial c^\mathrm{bulk}}{\partial t} = D_p\cdot \nabla^2 c^\mathrm{bulk} \tag{2} \end{equation}

Connection with experimentally accessible quantities

The bulk water concentration in the effective porosity model relates to the experimentally accessible concentration as

\begin{equation} \label{eq:cbar_epsilon} \bar{c} = \frac{\epsilon_\mathrm{eff}}{\phi} c^\mathrm{bulk} \tag{3} \end{equation}

where $\phi$ is the physical porosity of the sample. Since a bulk water concentration varies continuously across interfaces to external solutions, we have $c^\mathrm{bulk}(0) = c^\mathrm{source}$ at the source reservoir, giving

\begin{equation} \label{eq:cbar_epsilon0} \frac{\bar{c}(0)}{c^\mathrm{source}} = \frac{\epsilon_\mathrm{eff}} {\phi} \tag{4} \end{equation}

This equation shows that the effective porosity parameter quantifies the anion equilibrium concentration that we want to extract. That is not to say that the model is valid (more on that later), but that we can use eq. 4 to translate reported model parameters to an experimentally accessible quantity.

In principle, we could finish the analysis here, and use eq. eq. 4 as our main result. But most researchers do not evaluate the effective porosity in the direct way suggested by this equation (they may not even measure $\bar{c}$). Instead, they evaluate $\epsilon_\mathrm{eff}$ from a fitting procedure that also includes the diffusivity as a parameter. It is therefore fruitful to also include the transport aspects of the through-diffusion test in our analysis.

From closed-cell diffusion tests, we know that the clay concentration evolves according to Fick’s second law, both for many cations and anions. We will therefore take as an experimental fact that $\bar{c}$ evolves according to

\begin{equation} \label{eq:Fick2_exp} \frac{\partial \bar{c}}{\partial t} = D_\mathrm{macr.} \nabla^2 \bar{c} \tag{5} \end{equation}

This equation defines the diffusion coefficient $D_\mathrm{macr.}$, which should be understood as an empirical quantity.

Combining eqs. 3 and 2 shows that $D_p$ governs the evolution of $\bar{c}$ in the effective porosity model (if $\epsilon_\mathrm{eff}/\phi$ can be considered a constant). A successful fit of the effective porosity model to experimental data thus provides an estimate of $D_\mathrm{macr.}$ (cf. eq. 5), and we may write

\begin{equation} D_p = D_\mathrm{macr.} \tag{6} \end{equation}

With the additional assumption of constant reservoir concentrations, eq. 2 has a relatively simple analytical solution, and the corresponding outflux reads

\begin{equation} \label{eq:flux_analytic} j^\mathrm{out}(t) = j^\mathrm{ss} \left ( 1 + 2\sum_{n=1}^\infty \left (-1 \right)^n e^{-\frac{\pi^2n^2 D_\mathrm{p} t}{L^2}} \right ) \tag{7} \end{equation}

where $j^\mathrm{ss}$ is the steady-state flux. In steady-state, $c^\mathrm{bulk}$ is distributed linearly across the sample, and we can express the gradient in eq. 1 using the reservoir concentrations, giving

\begin{equation} j^\mathrm{ss} = \epsilon_\mathrm{eff} \cdot D_\mathrm{p} \cdot \frac{c^\mathrm{source}}{L} \tag{8} \end{equation}

where we have assumed zero target concentration.

Treating $j^\mathrm{ss}$ as an empirical parameter (it is certainly accessible experimentally), and using eq. 6, we get another expression for $\epsilon_\mathrm{eff}$ in terms of experimentally accessible quantities

\begin{equation} \epsilon_\mathrm{eff} = \frac{j^\mathrm{ss}\cdot L}{c^\mathrm{source} \cdot D_\mathrm{macr.} } \tag{9} \end{equation}

This relation (together with eqs. 4 and 6) demonstrates that if we fit eq. 7 using $D_p$ and $j^\mathrm{ss}$ as fitting parameters, the equilibrium relation we seek is given by

\begin{equation} \label{eq:exp_estimate} \frac{\bar{c}(0)}{c^\mathrm{source}} = \frac{j^\mathrm{ss}\cdot L} {\phi \cdot c^\mathrm{source} \cdot D_\mathrm{macr.} } \tag{10} \end{equation}

This procedure may look almost magical, since any explicit reference to the effective porosity model has now disappeared; eq. 10 can be viewed as a relation involving only experimentally accessible quantities.

But the validity of eq. 10 reflects the empirical fact that the (steady-state) flux can be expressed using the gradient in $\bar{c}$ and the physical porosity. The effective porosity model can be successfully fitted to anion through-diffusion data simply because it complies with this fact. Consequently, a successful fit does not validate the effective porosity concept, and essentially any description for which the flux can be expressed as $j = -\phi\cdot D_p \cdot \nabla\bar{c}$ will be able to fit to the data.

We may thus consider a generic model for which eq. 5 is valid and for which a steady-state flux is related to the external concentration difference as

\begin{equation} \label{eq:jss_general} j_\mathrm{ss} = – \beta\cdot D_p \cdot \frac{c^\mathrm{target} – c^\mathrm{source}}{L} \tag{11} \end{equation}

where $\beta$ is an arbitrary constant. Fitting such a model, using $\beta$ and $D_p$ as parameters, will give an estimate of $\bar{c}(0)/c^\mathrm{source}$ ($=\beta / \phi$).

Note that the system does not have to reach steady-state — eq. 11 only states how the model relates a steady-state flux to the reservoir concentrations. Moreover, the model being fitted is generally numerical (analytical solutions are rare), and may account for e.g. possible variation of concentrations in the reservoirs, or transport in the filters connecting the clay and the external solutions.

The effective porosity model emerges from this general description by interpreting $\beta$ as quantifying the volume of a bulk water phase within the bentonite sample. But $\beta$ can just as well be interpreted e.g. as an ion equilibrium coefficient ($\phi\cdot \Xi = \beta$), showing that this description also complies with the homogeneous mixture model.

Additional comments on the effective porosity model

The effective porosity model can usually be successfully fitted to anion through-diffusion data (that’s why it exists). The reason is not because the data behaves in a manner that is difficult to capture without assuming that anions are exclusively located in a bulk water domain, but simply because this model complies with eqs. 5 and 11. We have seen that also the homogeneous mixture model — which makes the very different choice of having no bulk water at all within the bentonite — will fit the data equally well: the two fitting exercises are equivalent, connected via the parameter identification $\epsilon_\mathrm{eff} \leftrightarrow \phi\cdot\Xi$.

Given the weak validation of the effective porosity model, I find it concerning that most anion through-diffusion studies are nevertheless reported in a way that not only assumes the anion-accessible porosity concept to be valid, but that treats $\epsilon_\mathrm{eff}$ basically as an experimentally measured quantity.

Perhaps even more remarkable is that authors frequently treat the effective porosity model as was it some version of the traditional diffusion-sorption model. This is often done by introducing a so-called rock capacity factor $\alpha$ — which can take on the values $\alpha = \phi + \rho\cdot K_d$ for cations, and $\alpha = \epsilon_\mathrm{eff}$ for anions — and write $D_e = \alpha D_a$, where $D_a$ is the “apparent” diffusion coefficient. The reasoning seems to go something like this: since the parameter in the governing equation in one model can be written as $D_e/\epsilon_\mathrm{eff}$, and as $D_e/(\phi + \rho\cdot K_d)$ in the other, one can view $\epsilon_\mathrm{eff}$ as being due to negative sorption ($K_d < 0$).

But such a mixing of completely different mechanisms (volume restriction vs. sorption) is just a parameter hack that throws most process understanding out the window! In particular, it hides the fact that the effective porosity and diffusion-sorption models are incompatible: their respective bulk water domains have different volumes. Furthermore, this lumping together of models has led to that anion diffusion coefficients routinely are reported as “apparent”, although they are not; the underlying model contains a pore diffusivity (eq. 2). As I have stated before, the term “apparent” is supposed to convey the meaning that what appears as pure diffusion is actually the combined result of diffusion, sorption, and immobilization. Sadly, in the bentonite literature, “apparent diffusivity” often means “actual diffusivity”.

Footnotes

[1] For anions, the total amount is relatively easy to measure by e.g. aqueous extraction. Cations, on the other hand, will stick to the clay, and need to be exchanged with some other type of cation (not initially present). In any case, the total amount of a species ($n$) can in principle be obtained experimentally, in an unambiguous manner.

[2] Another reasonable choice would be to divide by the total sample volume.

[3] If the test is designed as to have a significant change of the source concentration, it is a good idea to also measure the concentration evolution in this reservoir.

[4] Here we assume that the transfer resistance of the filter is negligible.

[5] Provided that the rest of the aqueous chemistry remains constant, which is not always the case. For instance, cation exchange may occur during the course of the test, if the set-up involves more than one type of cation, and there may be ongoing mineral dissolution.

The danger of log-log plots — measuring and modeling “apparent” diffusivity

Leave a reply

In a previous blog post, we discussed how the diffusivity of simple cations¹ has a small, or even negligible, dependence on background concentration (or, equivalently, on $K_d$), and how this observation motivates modeling compacted bentonite as a homogeneous system, containing only interlayer pores.

Despite the indisputable fact that “$D_a$”² for simple ions does not depend much on $K_d$, the results have seldom been modeled using a homogeneous bentonite model. Instead there are numerous attempts in the bentonite literature to both measure and model a variation of “$D_a$” with $K_d$, usually with a conclusion (or implication) that “$D_a$” depends significantly on $K_d$. In this post we re-examine some of these studies.

The claimed $K_d$-dependency is often “supported” by the so-called surface diffusion model. I have previously shown that this model is incorrect.³ Here we don’t concern ourselves with the inconsistencies, but just accept the resulting expression as the model to which authors claim to fit data. This model expression is

\begin{equation} D_a = \frac{D_p + \frac{\rho K_d}{\phi} D_s}{1+\frac{\rho K_d}{\phi}} \tag{1} \end{equation}

where $D_p$ and $D_s$ are individual domain diffusivities for bulk water and surface regions, respectively, $\rho$ is dry density, $\phi$ porosity, and $K_d$, of course, is assumed to quantify the distribution of ions between bulk water and surfaces as $s = K_d\cdot c^\mathrm{bulk}$, where $s$ is the amount of ions on the surface (per unit dry mass), and $c^\mathrm{bulk}$ is the corresponding bulk water concentration.

Muurinen et al. (1985)

Muurinen et al. (1985) measured diffusivity in high density bentonite samples at various background concentrations, using a type of closed-cell set-up. They also measured corresponding values of $K_d$ in batch “sorption” tests. The results for cesium, in samples with density in the range $1870 \;\mathrm{kg/m^3}$ — $2030 \;\mathrm{kg/m^3}$, are presented in the article in a figure similar to this:

The markers show experimental data, and the solid curve shows the model (eq. 1) with $D_p = 1.2 \cdot 10^{-10}\;\mathrm{m^2/s}$⁴ and $D_s = 4.3\cdot 10^{-13}\;\mathrm{m^2/s}$.

The published plot may give the impression of a systematic variation of $D_a$ for cesium, and that this variation is captured by the model. But the data is plotted with a logarithmic y-axis, which certainly is not motivated. Let’s see how the plot looks with a linear y-axis (we keep the logarithmic x-axis, to clearly see the model variation).

Now the impression is quite different: this way of plotting reveals that the experimental data only cover a part where the model does not vary significantly. With the adopted range on the x-axis (as used in the article) we actually don’t see the full variation of the model curve. Extending the x-axis gives the full picture:

With the full model variation exposed, it is evident that the model fits the data only in a most superficial way. The model “fits” only because it has insignificant $K_d$-dependency in the covered range, in similarity with the measurements.

The defining feature of the model is that the diffusivity is supposed to transition from one specific value at high $K_d$, to a significantly different value at low $K_d$. As no such transition is indicated in the data, the above “fit” does not validate the model.

Muurinen et al. (1985) also measured diffusion of strontium in two samples of density $1740 \;\mathrm{kg/m^3}$. The figures below show the data and corresponding model curve.

The left diagram is similar to how the data is presented in the article, while the right diagram utilizes a linear y-axis and shows the full model variation. The line shows the surface diffusion model with parameters $D_p = 1.2 \cdot 10^{-10}\;\mathrm{m^2/s}$ and $D_s = 8.8 \cdot 10^{-12}\;\mathrm{m^2/s}$. In this case it is clear even from the published plot that the experimental data shows no significant variation.

The only reasonable conclusion to make from the above data is that cesium and strontium diffusivity does not significantly depend on $K_d$ (which implies a homogeneous system). This is actually also done in the article:

The apparent diffusivities of strontium and cesium do not change much when the salt concentration used for the saturation of the samples is changed and the sorption factors change. The surface diffusion model agrees fairly well with the observed diffusion-sorption behaviour.

I agree with the first sentence but not with the second. In my mind, the two sentences contradict each other. From the above plots, however, it is trivial to see that the surface diffusion model does not agree (in any reasonable sense) with observations.

Eriksen et al. (1999)

Although Muurinen et al. (1985) concluded insignificant $K_d$-dependency on the diffusion coefficients for strontium and cesium, researchers have continued throughout the years to fit the surface diffusion model to experimental data on these and other ions.

Eriksen et al. (1999) present old and new diffusion data for strontium and cesium (and sodium), fitted and plotted in the same way as in Muurinen et al. (1985). Here are the evaluated diffusivities for cesium plotted against evaluated $K_d$, as presented in the article, and re-plotted in different ways with a linear y-scale:

The curve shows the surface diffusion model (eq. 1), with parameters $D_p = 8 \cdot 10^{-10}\;\mathrm{m^2/s}$ and $D_s = 6 \cdot 10^{-13}\;\mathrm{m^2/s}$. The points labeled “Eriksen 99” are original data obtained from through-diffusion tests on “MX-80” bentonite at dry density 1800 $\mathrm{kg/m^3}$.⁵ The source for the data points labeled “Muurinen 94” is the PhD thesis of A. Muurinen.⁶

The upper left plot shows the data as presented in the article; again, a logarithmic y-axis is used. In this case, a zoomed-in view with a linear y-axis (upper right diagram) may still give the impression that the data has a systematic variation that is captured by the model. But viewing the whole range reveals that the model is fitted to data where variation is negligible (bottom diagrams), just as in Muurinen et al. (1985).

Data and model for strontium presented in Eriksen et al. (1999) look like this:

The model (line) has parameters $D_p = 3 \cdot 10^{-10}\;\mathrm{m^2/s}$ and $D_s = 1 \cdot 10^{-11}\;\mathrm{m^2/s}$, and the source for the data points labeled “Eriksen 84” is found here.

In this case, not even the diagram presented in the article (left) seems to support the promoted model. This is also confirmed when utilizing a linear y-axis, and showing the full model variation (right diagram).

Eriksen et al. (1999) conclude that strontium diffusivities are basically independent of $K_d$, but claim, in contrast to Muurinen et al. (1985), that cesium diffusivity depends significantly on $K_d$:

[I]n the $K_d$ interval 0.01 to 1 the apparent $\mathrm{Cs}^+$ diffusivity decreases by approximately one order of magnitude whereas for $\mathrm{Na}^+$ and $\mathrm{Sr}^{2+}$ the apparent diffusivity is virtually constant.

They also claim that the surface diffusion model fits the data:

$D_\mathrm{a}$ curves for $\mathrm{Cs}^+$ and $\mathrm{Sr}^{2+}$, calculated using a Eq. (6) [eq. 1 here], are plotted in Fig. 4. As can be seen, good fits to experimental data were obtained […]

Note that the variation in the model for cesium is motivated by three data points with relatively high diffusivity and basically the same $K_d \sim 0.05\;\mathrm{m^3/kg}$. It seems like the model has been fitted to these points, while the point at $K_d \sim 0.02\;\mathrm{m^3/kg}$ has been mainly neglected. The resulting model has a huge bulk water diffusivity ($D_p$), which is about 7 times larger than in the corresponding fit in Muurinen et al. (1985), and only 2.5 times smaller than the diffusivity for cesium in pure water.

Note that, if you claim that the surface diffusion model fits in this case, you implicitly claim that the observed variation — which still is negligible on the scale of the full model variation — is caused by the influence of this enormous (for a 1800 $\mathrm{kg/m^3}$ sample) bulk pore water diffusivity; with a more “reasonable” value for $D_p$, the model no longer fits. There are consequently valid reasons to doubt that the claimed $K_d$ dependence is real. We will return to this fit in the next section.

Gimmi & Kosakowski (2011)

We have now seen several examples of authors erroneously claiming (or implying) that a surface diffusion model is valid, when the actual data for “$D_a$” has no significant $K_d$-dependency. For reasons I cannot get my head around, this flawed treatment is still in play.

Rather than identifying the obvious problem with the previously presented fits, Gimmi and Kosakowski (2011) instead extended the idea of expressing the diffusivity as a function of $K_d$ by using scaled, dimensionless quantities

\begin{equation} D_\mathrm{arw} = \frac{D_\mathrm{a}\tau_w}{D_0} \tag{2} \end{equation}

\begin{equation} \kappa = \frac{\rho K_d}{\phi} \end{equation}

where $D_0$ is the corresponding diffusivity in pure water and $\tau_w$ is the “tortuosity factor” for water in the system of interest. This factor is simply the ratio between the water diffusivity in the system of interest and the water diffusivity in pure water (I have written about the problem with factors like this here).

The idea — it seems — is that using $D_\mathrm{arw}$ and $\kappa$ as variables should make it possible to directly compare the mobility of a given species in systems differing in density, clay content, etc.

Even though it makes some sense that the diffusivity of a specific species scales with the diffusivity of water in the same system, the above procedure inevitably introduces more variation in the data — both because an additional measured quantity (water diffusivity) is involved when evaluating the scaled diffusivity, but also because water diffusivity may depend differently on density as compared with the diffusivity of the species under study.

Also Gimmi and Kosakowski (2011) use the flawed surface diffusion model for analysis, and their expression for $D_\mathrm{arw}$ is

\begin{equation} D_\mathrm{arw} = \frac{1+\mu_s\kappa}{1+\kappa} \tag{3} \end{equation}

where $\mu_s = D_s\tau_w/D_0$ is a “relative surface mobility”. This equation is obtained from eq. 1, by dividing by $D_p$ and assuming $D_p = D_0/\tau_w$.

Gimmi and Kosakowski (2011) fit eq. 3 to a large set of collected data, measured in various types of material, including bentonites, clay rocks, and clayey soils. This is their result for cesium⁷ (the model curve is eq. 3 with $\mu_s = 0.031$⁸)

Viewed as a whole, this data is more scattered as compared with the previous studies. This is reasonably an effect of the larger diversity of the samples, but also an effect of multiplying the “raw” diffusion coefficient with the factor $\tau_w$ (eq. 2).

Just as in the previous studies we have looked at, the published plot (similar to the left diagram) may give the impression of a systematic variation of the diffusivity with $K_d$ (it contains partly the same data). But just as before, a linear y-axis (right diagram) reveals that the model is fitted only to data where variation is negligible.

Note that the three data points that contributed to the majority of the variation in the fitted model in Eriksen et al. (1999) here appear as outliers.⁹ The variation with $K_d$ for cesium claimed in that study is thus invalidated by this larger data set.

As we have noted already, the only reasonable conclusion to draw from this data is that there is no systematic $K_d$-dependency on diffusivity of cesium or strontium, and that it does not — in any reasonable sense — fit the surface diffusion model. Yet, also Gimmi and Kosakowski (2011) imply that the surface diffusion is valid:

The data presented here show a general agreement with a simple surface diffusion model, especially when considering the large errors associated with the $D_\mathrm{erw}$ and $D_\mathrm{arw}$.

This paper, however, contains an even worse “fit” to strontium data, as compared to the earlier studies (the left diagram is similar to the how it is presented in the article, the right diagram uses a linear y-axis; the line is eq. 3 with $\mu_s = 0.24$⁸):

This data does not suggest a variation in accordance with the adopted model even when plotted in a log-log diagram. With a linear y-axis, the dependence rather seems to be the opposite: $D_\mathrm{arw}$ appears to increase with $\kappa$. However, I suspect that this is a not a “real” dependence, but rather an effect of trying to construct a “relative” diffusivity; note that while $\kappa$ spans four orders of magnitude, $D_\mathrm{arw}$ scatters only by a factor of 5 or 6. Nevertheless, how this data can be claimed to show “general agreement” with the surface diffusion model is a mystery to me.

The view is similar for sodium (the left diagram is similar to the how it is presented in the article, the right diagram uses a linear y-axis; the line is eq. 3 with $\mu_s=0.52$⁸):

Even if the model in this case only displays minor variation, it can hardly be claimed to fit the data: again, the data suggests a diffusivity that increases with $\kappa$. But a significant amount of these data points have $D_\mathrm{arw} > 1$, which is not likely to be true, as it indicates that the relative mobility for sodium is larger than for water. Consequently, the major contribution of the variation seen in this data is most probably noise.

Gimmi and Kosakowski (2011) also examined diffusivity for calcium, and the data looks like this (the left diagram is similar to the how it is presented in the article, the right diagram uses a linear y-axis; the line is eq. 3 with $\mu_s=0.1$⁸):

Here it looks like the data, to some extent, behaves in accordance with the model also when plotted with linear y-axis covering the full model variation. However, there are significantly less data reported for calcium (as compared with cesium, strontium, and sodium) and the model variation is supported only by a few data points¹⁰. I therefore put my bet on that if calcium diffusivity is studied in more detail, the dependence suggested by the above plot will turn out to be spurious.¹¹

Some thoughts

I am more than convinced that the only reasonable starting point for modeling saturated bentonite is a homogeneous description. I had nevertheless expected to at least have to come up with an argument against the multi-porous view put forward in the considered publications (and in many others). I am therefore quite surprised to find that this argument is already provided by the data in the very same publications (and even by the statements, sometimes): there is nothing in the data here reviewed that seriously suggests that cation diffusion is influenced by a heterogeneous pore structure.

Still, the unsupported idea that cations in compacted bentonite are supposed to diffuse in two (or more) different types of water domains has evidently propagated through the scientific literature for decades, and a multi-porous view is mainstream in modern bentonite research. It is difficult to not feel disheartened when faced with this situation. What would it take for researchers to begin scrutinize their assumptions? Is nobody interested in the topics we are supposed to study?

Footnotes

[1] Unfortunately, a quantity which by many is incorrectly interpreted as an “apparent” diffusivity.

[2] I use quotation marks to indicate that $D_a$ is a parameter in the traditional diffusion-sorption model, a model not valid for compacted bentonite. Still, this parameter is often reported as if it was a directly measured quantity.

[3] I have also derived a correct version of the surface diffusion model, which does not involve apparent diffusivity.

[4] The article states $\epsilon D_p = 3.5\cdot 10^{-11}\; \mathrm{m^2/s}$, where $\epsilon$ is the porosity. $D_p = 1.2\cdot 10^{-10} \; \mathrm{m^2/s}$ corresponds to $\epsilon = 0.29$.

[5] In this study, both $K_d$ and $D_a$ were evaluated by fitting the traditional diffusion-sorption model to concentration measurements.

[6] I have had no access to this document, and I have not verified e.g. sample density (this data set is different from that presented in the previous section). Instead, I have read these values from the diagram in Eriksen et al. (1999).

[7] They actually divide their cesium data into two categories, which show quite different mobility. The data shown here — which includes bentonite samples — is for systems categorized as being “non-illite” or having Cs concentration above “trace”.

[8] According to the article table, the fitted values for $\mu_s$ are 0.52 (Na), 0.39 (Sr), 0.087 (Ca), and 0.015 (Cs). The plotted lines, however, appear to instead use what is listed as “mean $\mu_s$”. Here, I have used these $\mu_s$-values: 0.52 (Na), 0.24 (Sr), 0.1 (Ca), and 0.031 (Cs).

[9] This cluster contains a fourth data point, from Jensen and Radke (1988).

[10] All data for calcium is essentially from only two different sources: Staunton (1990) and Oscarsson (1994).

[11] It would also be more than amazing if it turns out — after it is verified that Cs, Na, and (especially) Sr show no significant $K_d$ dependence — that Ca diffusivity actually varies in accordance with the flawed surface-diffusion model!

Swelling pressure, part V: Suction

Leave a reply

There are several “descriptions” of bentonite swelling. While a few of them actually denies any significant role played by the exchangeable cations, most of these descriptions treat the exchangeable ions as part of an osmotic system. I have earlier discussed how the terms “osmotic” or “osmosis” may cause some confusion in different contexts, and discussed the confusion surrounding the treatment of electrostatic forces.

In this blog post I discuss the description of bentonite swelling often adopted in the fields of soil mechanics and geotechnical engineering. In particular, we focus on the concept of suction, which is central in these research fields, while being basically absent in others.

As far as I understand, suction is just the water chemical potential “disguised” as a pressure variable; although I have trouble finding clear-cut definitions, it seems clear that suction is directly inherited from the “water potential” concept, which has been central in soil science for a long time. Applied to bentonite, the geotechnical description is thus not principally different from the osmotic approach that I have presented previously. But the way the suction concept is (and isn’t) applied may cause unnecessary confusion regarding the swelling mechanisms. I think a root for this confusion is that suction involves both osmotic and capillary mechanisms.

Matric suction (capillary suction)

Matric suction is typically associated with capillarity, a fundamental mechanism in many conventional soil materials under so-called unsaturated conditions. A conventional soil with a significant amount of small enough pores shows capillary condensation, i.e. it contains liquid water below the condensation point for ordinary bulk water. Naturally, the equilibrium vapor pressure increases with the amount of water in the soil, as the pores containing liquid water become larger. For conventional soils, it therefore makes sense to speak of the degree of saturation of a sample, and to relate saturation and equilibrium vapor pressure by means of a water retention curve. Underlying this picture is the notion that the solid parts constitute a “soil skeleton” (the matrix), and that the soil can be viewed as a vessel that can be more or less filled with water.

The pressure of the capillary water is lower than that of the surrounding air, and is related to the curvature of the interfaces between the two phases (menisci), as expressed by the Young-Laplace equation. For a spherically symmetric meniscus this equation reads

\begin{equation} \Delta p = p_w – p_a = \frac{2\sigma}{r} \tag{1} \end{equation}

where $p_a$ and $p_w$ denote the pressures of air and capillary water, respectively, $\sigma$ is the surface tension, and $r$ is the radius of curvature of the interface. The sign of $r$ depends on whether the interface bulges inwards (“concave”, $r<0$) or outwards (“convex”, $r>0$). For capillary water, $r$ is negative and $\Delta p$ — which is also called the Laplace pressure — is a negative quantity.

As far as I understand, matric suction is simply defined as the negative Laplace pressure, i.e.

\begin{equation} s_m = p_a – p_w \tag{2} \end{equation}

With this definition, suction has a straightforward physical meaning as quantifying the difference in pressure of the two fluids occupying the pore space, and clearly relates to the everyday use of the word.

Suction — in this capillary sense — gives a simple principal explanation for (apparent) cohesion in e.g. unsaturated sand: individual grains are pushed together by the air-water pressure difference, as schematically illustrated here (the yellow stuff is supposed to be two grains of sand, and the blue stuff water)

Net force for two sand grains exposed to matric suction

It is reasonable to assume that the net force transmitted by the soil skeleton — usually quantified using the concept of effective stress — governs several mechanical properties of the soil sample, e.g. shear strength. The above description also makes it reasonable to assume that effective stress depends on suction.

Thus, in unsaturated conventional soil are quantities like degree of saturation, pore size distribution, (matric) suction, effective stress, and shear strength very much associated. Another way of saying this is that there is an optimal combination of water content and particle size distribution for constructing the perfect sand castle.

The chemical potential of the capillary water is related to matric suction. Choosing pure bulk water under pressure $p_a$¹ as reference, the chemical potential of the liquid phase in the soil is obtained by integrating the Gibbs-Duhem equation from $p_a$ to $p_w$

\begin{equation} \mu = \mu_0 + \int_{p_a}^{p_w}v dP = \mu_0 + v\cdot \left(p_w – p_a \right) = \mu_0 -v\cdot s_m \tag{3} \end{equation}

where $\mu_0$ is the reference chemical potential, $v$ is the molar volume of water, and we have assumed incompressibility.

The above expression shows that matric suction in this case directly quantifies the (relative) water chemical potential. Note, however, that eq. 3 does not define matric suction; $s_m$ is defined as a pressure difference between two phases (eq. 2), and happens to quantify the chemical potential under the present circumstances (pure capillary water).

A chemical potential can generally be expressed in terms of activity ($a$)

\begin{equation} \mu = \mu_0 + RT \ln a \tag{4} \end{equation}

For our case, water activity is to a very good approximation equal to relative humidity, the ratio between the vapor pressures in the state under consideration and in the reference state, i.e. $a = p_v/p_{v,0}$. Combining eqs. 3 and 4, we see that the vapor pressure in this case is related to matric suction as

\begin{equation} \frac{p_v}{p_{v,0}} = e^{-v\cdot s_m/RT} \end{equation}

Using the Young-Laplace equation (eq. 1) for $s_m$ we can also write this as

\begin{equation} \frac{p_v}{p_{v,0}} = e^{\frac{2v\sigma}{RTr}} \end{equation}

This is the so-called Kelvin equation, which relates the equilibrium vapor pressure to the curvature of an air-pure water interface. Note that, since $r<0$ for capillary water, the vapor pressure is lower than the corresponding bulk value ($p_v < p_{v,0}$).

Osmotic suction and total suction

So far, we have discussed suction in a capillary context, and related it to water chemical potential or vapor pressure. Now consider how the picture changes if the pores in our conventional soil contain saline water. Matric suction — i.e. the actual pressure difference between the pore solution and the surrounding air, sticking with eq. 2 as the definition — is in general different from the pure water case, because solutes influence surface tension. Also, water activity (vapor pressure) is different from the pure water case, but there is no longer a direct relation between water activity and matric suction, because water activity is independently altered by the presence of solutes.

The water chemical potential of a saline bulk solution (i.e. with no capillary effects), can be written in terms of the osmotic pressure, $\pi(c)$

\begin{equation} \mu(c) = \mu_0 – v\cdot\pi(c) \tag{5} \end{equation}

where we have assumed a salt concentration $c$, and indicated that the osmotic pressure, and hence the chemical potential, depends on this concentration.

Although eq. 5 is of the same form as eq. 3, matric suction and osmotic pressure are very different quantities. The former is defined under circumstances where an actual pressure difference prevail between the air and water phases. In contrast, there is no pressure difference between the phases in a container containing both a solution and a gas phase. $\pi(c)$ corresponds to the elevated pressure that must be applied for the solution to be in equilibrium with pure water kept at the reference pressure.

Despite the different natures of matric suction and osmotic pressure, the fields of geotechnical engineering and soil mechanics insist on also referring to $\pi(c)$ as a suction variable: the osmotic suction. Similarly, total suction is defined as the sum of matric and osmotic suction

\begin{equation} \Psi = s_m + \pi(c) \end{equation}

These definitions seem to have no other purpose than to be able to write the water chemical potential generally as

\begin{equation} \mu = \mu_0 -v\cdot \Psi \tag{6} \end{equation}

Total suction is thus de facto defined simply as the (relative) value of the water chemical potential, expressed as a pressure (I think this is completely analogous to “total water potential” in soil science).

Eq. 6 shows that $\Psi$ is directly related to water activity, or vapor pressure, and we can write

\begin{equation} \frac{p_v}{p_{v,0}} = e^{-v\cdot \Psi/RT} \tag{7} \end{equation}

This relation is quite often erroneously referred to as the Kelvin equation (or “Kelvin’s law”) in the bentonite literature. But note that the above equation just restates the definition of water activity, because $v\cdot\Psi$ cannot be reduced to anything more concrete than the relative value of the water chemical potential. The Kelvin equation, on the other hand, expresses something more concrete: the equilibrium vapor pressure for a curved air-water interface. Some clay literature refer to the above relation as the “Psychrometric law”, but that name seems not established in other fields.²

A definition is motivated by its usefulness, and total change in water chemical potential is of course central when considering e.g. moisture movement in soil. My non-geotechnical brain, however, is not fond of extending the “suction” variable in the way outlined above. To start with, there is already a variable to use: the water chemical potential. Also, “total suction” no longer has the direct relation to the everyday use of the word suction: there is no “sucking” going on in a saline bulk solution,³ while in a capillary there is. Furthermore, with a saline pore solution there is no direct relation between (total) suction and e.g. effective stress or shear strength.

Although both matric suction and osmotic pressure under certain circumstances can be measured in a direct way, it seems that (total) suction usually is quantified by measuring/controlling the vapor pressure with which the soil sample is in equilibrium. Actually, one of the more comprehensive definitions of various “suctions” that I have been able to find — in Fredlund et al. (2012) — speaks only of various vapor pressures (although based on the capillary and osmotic concepts):⁴

Matric or capillary component of free energy: Matric suction is the equivalent suction derived from the measurement of the partial pressure of the water vapor in equilibrium with the soil-water relative to the partial pressure of the water vapor in equilibrium with a solution identical in composition with the soil-water.

Osmotic (or solute) component of free energy: Osmotic suction is the equivalent suction derived from the measurement of the partial pressure of the water vapor in equilibrium with a solution identical in composition with the soil-water relative to the partial pressure of water vapor in equilibrium with free pure water.

Total suction or free energy of soil-water: Total suction is the equivalent suction derived from the measurement of the partial pressure of the water vapor in equilibrium with the soil-water relative to the partial pressure of water vapor in equilibrium with free pure water.

It seems that such operational definitions of suction has made the term synonymous with “vapor pressure depression” in large parts of the bentonite scientific literature.

Suction in bentonite

In the above discussion we had mainly a conventional soil in mind. When applying the suction concepts to bentonite,⁵ I think there are a few additional pitfalls/sources for confusion. Firstly, note that the definitions discussed previously involve “a solution identical in composition with the soil-water”. But soil-water that contains appreciable amounts of exchangeable ions — as is the case for bentonite — cannot be realized as an external solution.

It seems that this “complication” is treated by assuming that an external solution in equilibrium with a bentonite sample is the soil-water (this is analogous to how many geochemists use the term “porewater” in bentonite contexts). Not surprisingly, this treatment has bizarre consequences. The conclusion for e.g. a salt free bentonite sample — which is in equilibrium with pure water — is that it lacks osmotic suction, and that its lowered vapor pressure (when isolated and unloaded) is completely due to matric suction! I think this is such an odd outcome that it is worth repeating: A system dominated by interlayer pores, containing dissolved cations at very high concentrations, is described as lacking osmotic pressure! It is not uncommon to find descriptions like this one (from Lang et al. (2019))

The total suction of unsaturated soils consists of matric and osmotic suctions (Yong and Warkentin, 1975; Fredlund et al., 2012; Lu and Likos, 2004). In clays, the matric suction is due to surface tension, adsorptive forces and osmotic forces (i.e. the diffuse double layer forces), whereas the osmotic suction is due to the presence of dissolved solutes in the pore water.

We apparently live in a world where “matric suction” consists of “osmotic forces”, while the same “osmotic forces” do not contribute to “osmotic suction”. Except when the clay contains excess ions, in which case we have an arbitrary combination of the two “suctions” (note also that “osmotic suction” and “osmotic swelling” are two quite different things).

Although the above consequence is odd, it is still only a matter of definition: accepting that “matric suction” involves osmotic forces (which I don’t recommend), the description may still be adequate in principle; after all, “total suction” quantifies the reduction of the water chemical potential.

But the focus on “matric suction” also reveals a conceptual view of bentonite structure that I find problematic: it suggests a first order approximation of bentonite as a conventional soil, i.e. as an assemblage of solid grains separated from an aqueous phase (and a gas phase). This “matric” view is fully in line with the idea of “free water” in bentonite, and it is quite clear that this is a prevailing view in the geotechnical, as well as in the geochemical, literature. For instance, with the formulation “the presence of dissolved solutes in the pore water” in the above quotation, the “pore water” the authors have in mind is a charge neutral bulk water solution.

With the “matric” conceptual view, the degree of saturation becomes a central variable in much soil mechanical analyses of bentonite. When dealing with actual unsaturated bentonite samples, I guess this makes sense, but once a sample is saturated this variable has lost much of its meaning.⁶ Consider e.g. the different expected behaviors if drying e.g. a water saturated metal filter or a saturated bentonite sample.

The different nature of drying a metal filter compared with drying a saturated bentonite sample

The equilibrium vapor pressure of both these systems is lower than the corresponding pure bulk water value. For the metal filter, the lowered water activity is of course due to capillarity, i.e. there is an actual pressure reduction in the water phase (matric suction!). When lowering the external vapor pressure below the equilibrium point (i.e. drying), capillary water migrates out of the filter, while the metal structure itself remains intact. In this case, as the system remains defined in a reasonable way, it is motivated to speak of the saturation state of the filter.

For a drying bentonite sample, the behavior is not as well defined, and depends on how the drying is performed and on initial water content. For a quasi-static process, where the external vapor pressure is lowered in small steps at an arbitrary slow rate, it should be clear that the entire sample will respond simply by shrinking. In this case it does not make much sense to speak of the sample as still being saturated, nor to speak of it as having become unsaturated.

For a more “violent” drying process, e.g. placing the bentonite sample in an oven at 105 °C , it is also clear that — rather than resulting in a neatly shrunken, dense piece of clay — the sample now will suffer from macroscopic cracks and other deformations. Neither in this case does it make much sense to try to define the degree of saturation, in relation to the sample initially put in the oven.

Note also that if we, instead of drying, increase the external vapor pressure from the initial equilibrium value, the metal filter will not respond much at all, while the bentonite sample immediately will begin to swell.

I hope that this example has made it clear, not only that the degree of saturation is in general ill-defined for bentonite, but also that a bentonite sample behaves more as an aqueous solution rather than as a conventional soil: if we alter external vapor pressure, an aqueous solution responds by either “swelling” (taking up water) or “shrinking” (giving off water). A main aspect of this conceptual view of bentonite — which we may call the “osmotic” view — is that water does not form a separate phase.⁷ This was pointed out e.g. by Bolt and Miller (1958) (referring to this type of system as an “ideal clay-water system”)

In contrast to the familiar case described is the ideal clay-water system in which the particles are not in direct contact but are separated by layers of water. Removal of water from such a system does not introduce a third phase but merely causes the particles to move closer to one another with the pores remaining water saturated.

From these considerations it follows that a generally consistent treatment is to relate bentonite water activity to water content, rather than to degree of saturation.

Another consequence of adopting a “matric” view of bentonite (i.e. to include osmotic forces in “matric suction”) is that “matric suction” loses its direct connection with effective stress. This can be illustrated by taking the “osmotic” view: just as the mechanical properties of an aqueous solution (e.g. viscosity) do not depend critically on whether or not it is under (osmotic) pressure, we should not expect e.g. bentonite shear strength to be directly related to swelling pressure.⁸

Footnotes

[1] Often, the air is at atmospheric pressure, in which case the reference is the ordinary standard state.

[2] The relatively common misspelling “Psychometric law” is kind of funny.

[3] The cautious reader may remark that saline solutions do “suck”, in terms of osmosis. But note the following: 1) Osmosis requires a semi-permeable membrane, separating the solution from an external water source. We have said nothing about the presence of such a component in the present discussion. The way osmotic suction sometimes is described in the literature makes me suspect that some authors are under the impression that the mere presence of a solute causes a pressure reduction in the liquid. 2) In the presence of a semi-permeable membrane, osmosis has no problem occurring without a pressure difference between between the two compartments. 3) For cases when the solution is acted on by an increased hydrostatic pressure, water is transported from lower to higher pressure. It is difficult to say that there is any “sucking” in such a process (I would argue that the establishment of a pressure difference is an effect, rather than a cause, in the case of osmosis) 4) The idea that a solution has a well-defined partial water pressure is wrong.

[4] I’m still not fully satisfied with this definition: It may be noted that the definitions are somewhat circular (“matric suction is the equivalent suction…”), so they still require that we have in mind that “suction” also is defined in terms of a certain vapor pressure ratio (e.g. eq. 7). Note also that the headings speak of “free energy”. Perhaps I am nitpicking, but (free) energy is an extensive quantity, while suction (pressure) is intensive. Thus, “free energy” here really mean “specific free energy” (or “partial free energy”, i.e. chemical potential). I think the soil science literature in general is quite sloppy with making this distinction.

[5] “Bentonite” is used in the following as an abbreviation for bentonite and claystone, or any clay system with significant cation exchange capacity.

[6] If you press bentonite granules to form a cohesive sample you certainly end up with a system having both water filled interlayer pores and air-filled macropores (or perhaps an even more complex pore structure). This blog post mainly concerns saturated bentonite, by which I mean bentonite material which does not contain any gas phase. We can thus speak of saturated bentonite, although a degree of saturation variable is not well defined.

[7] Rather, montmorillonite and water form a homogeneous mixture.

[8] However, bentonite strength relates indirectly to swelling pressure (under specific conditions) because both quantities depends on a third: density.

Donnan equilibrium and the homogeneous mixture model

Leave a reply

We can directly apply the homogeneous mixture model for bentonite to isolated systems — e.g. closed-cell diffusion tests — as discussed previously. For systems involving external solutions we must also handle the chemical equilibrium at solution/bentonite interfaces.

I have presented a framework for calculating the chemical equilibrium between an external solution and a bentonite component in the homogeneous mixture model here. In this post I will discuss and illustrate some aspects of that work.

Overview

We assume a homogeneous bentonite domain in contact with an external solution, with the clay particles prevented from crossing the domain interface. For real systems, this corresponds to the frequently encountered set-up with bentonite confined in a sample holder by means of e.g. a metal filter. From the assumptions of the homogeneous model — that all ions are mobile and allowed to cross the domain interface — it follows that the type of equilibrium to consider is the famous Donnan equilibrium. I have discussed the Donnan effect and its relevance for bentonite quite extensively here.

Since the adopted model assumes a homogeneous bentonite domain, the only region where Donnan equilibrium comes into play is at the interface between the bentonite and the external solution. This is quite different from how Donnan equilibrium calculations are implemented in many multi-porous models, where the equilibrium is internal to the clay — between assumed “macro” and “micro” compartments of the pore structure. The need for performing Donnan equilibrium calculations is thus minimized in the homogeneous mixture model (as mentioned, isolated systems require no such calculations). Note also that the semi-permeable mechanism in multi-porous models is required to act on the pore-scale. I have never seen any description or explanation how such a mechanism is supposed to work.¹ In the homogeneous mixture model, on the other hand, the semi-permeable interface corresponds directly to a macroscopic and experimentally well-defined component: the confining filter.

The problem to be solved can be illustrated like this

Schematic illustration of an external solution in contact with a homogeneous bentonite domain

The aim is to relate the set of species concentrations in the external solution ($\{c_i^\mathrm{ext}\}$) to those in the clay domain ($\{c_i^\mathrm{int}\}$) when the system is in equilibrium. This is done by applying the standard approach to Donnan equilibrium, as found in textbooks on the subject. If there is anything “radical” about this framework, it is thus not in the way Donnan equilibrium is implemented, but rather in treating bentonite as a single phase: this approach is formally equivalent to assuming the bentonite to be an aqueous solution.

Chemical equilibrium

I prefer to formulate the Donnan equilibrium framework in a way that separates effects due to difference in the local chemical environment from effects due to differences in electrostatic potential between the two compartments. An important reason for focusing on this separation is that the local environment affects the chemistry under all circumstances, while the (relative) value of the electrostatic potential only is relevant when bentonite is contacted with an external solution. We therefore express the chemical equilibrium as

\begin{equation} \frac{c_i^\mathrm{int}}{c_i^\mathrm{ext}} = \frac{\gamma_i^\mathrm{ext}}{\gamma_i^\mathrm{int}}\cdot e^{-\frac{z_iF\psi^\star}{RT}} \tag{1} \end{equation}

This formula is achieved by setting the electro-chemical potential equal for each species in the two compartments. Here $\gamma_i$ denotes the activity coefficient for species $i$, and $\psi^*$ is the electrostatic potential difference between the compartments, which we refer to as the Donnan potential.

I find it convenient to rewrite this expression using some fancy Greek letters

\begin{equation} \label{eq:chem_eq2} \Xi_i = \Gamma_i \cdot f_D^{-z_i} \tag{2} \end{equation}

Here I call $\Xi_i = c_i^\mathrm{int}/c_i^\mathrm{ext}$ the ion equilibrium coefficient for species $i$. This quantity expresses the essence of ion equilibrium in the homogeneous mixture model, and will appear in many places in the analysis. $\Xi_i$ has two factors:

$\Gamma_i = \gamma_i^\mathrm{ext}/\gamma_i^\mathrm{int}$ expresses the chemical aspect of the equilibrium: when $\Gamma_i$ is large ($>1$), the species has a chemical preference for residing in the interlayer pores, and when $\Gamma_i$ is small ($<1$), the species has a preference for the external solution. In general, $\Gamma_i$ for any specific species $i$ is a function of all species concentrations in the system.
$f_D^{-z_i}$, where $f_D = e^{\frac{F\psi^\star}{RT}}$ is a dimensionless transformation of the Donnan potential (this is basically the Nernst equation), which we here call the Donnan factor. $f_D$ expresses the electrostatic aspect of the equilibrium, and is the same for all species. The effect on $\Xi_i$, however, is different for species of different charge number, because of the exponent $-z_i$ in the full expression.

I want to emphasize that eqs. 1 and 2 express the exact same thing: chemical equilibrium between the two compartments.

Illustrations

To get a feel for the quantity $\Xi$, here is a hopefully useful animation

Relation beteween internal and external concentration for varying Xi

It may also be helpful to see the influence of $f_D$ on the equilibrium. Since the Donnan potential is negative, $f_D$ is less than unity and typical values in relevant bentonite systems is $f_D \sim$ 0.01 — 0.4. Due to the exponent $-z_i$ in eq. 2, this influence on the equilibrium looks quite different for species with different valency. For mono- and di-valent cations, the behavior looks like this (here is put $\Gamma = 1$ for both species)

Variation of internal cation concentrations with varying Donnan factor

The typical behavior for cations is that the internal concentration is much larger than the corresponding external concentration (at $f_D = 0.01$ in the above animation, the internal concentration for the di-valent cation is enhanced by a factor $\Xi = 10 000$!). For anions, the internal concentration is instead lower than the external concentration,² as shown here ($\Gamma = 1$ for both species)

Variation of internal anion concentration with the Donnan factor

Equation for $f_D$

For a complete description, we need an equation for calculating $f_D$. This is derived by requiring charge neutrality in the two compartments and looks like

\begin{equation*} \sum_i z_i\cdot\Gamma_i \cdot c_i^\mathrm{ext} \cdot f_D^{-z_i} – c_{IL} = 0 \tag{3} \end{equation*}

where

\begin{equation*} c_{IL} = \frac{CEC}{F \cdot w} \end{equation*}

is the structural charge present in the clay (i.e. negative montmorillonite layer charge) expressed as a monovalent interlayer concentration. Here $CEC$ is the cation exchange capacity of the clay component, $w$ the water-to-solid mass ratio,³ and $F$ is the Faraday constant.

The way eq. 3 is formulated implies that the external concentrations should be used as input to the calculation. This is typically the case as the external concentrations are under experimental control.

In typical geochemical systems it is required to account for aqueous species with valency at least in the range -2 — +2 (e.g. $\mathrm{Ca}^{2+}$, $\mathrm{Na}^{+}$, $\mathrm{Cl}^{-}$, $\mathrm{SO_4}^{2-}$), which implies that the equation for calculating $f_D$ is generally a polynomial equation of degree four or higher.

An important special case is the 1:1 system — e.g. pure Na-montmorillonite contacted with a NaCl solution — which has an equation for $f_D$ of only degree two, and thus has a relatively simple analytical solution

\begin{equation*} f_D = \frac{c_{IL}}{2c^\mathrm{ext} \Gamma_\mathrm{Cl}} \left ( \sqrt{1+ \frac{4(c^\mathrm{ext})^2 \Gamma_\mathrm{Na}\Gamma_\mathrm{Cl}} {c_{IL}^2}} – 1 \right ) \end{equation*}

With the machinery in place for calculating the Donnan potential, here is an animation demonstrating the response in internal sodium and chloride concentrations as the external NaCl concentration is varied. In this calculation $c_{IL} = 2$ M, and $\Gamma_\mathrm{Na} = \Gamma_\mathrm{Cl} = 1$

Relation between internal and external Na and Cl concentrations

Comment on through-diffusion

To me, the last illustration makes it absolutely clear that Donnan equilibrium and the homogeneous mixture model provide the correct principal explanation for e.g. the behavior of tracer ions in through-diffusion tests. If you choose to relate the flux in through-diffusion tests to the external concentration difference — which is basically done in all published studies, via the parameter $D_e$ — you will evaluate large “diffusivities” for cations and small “diffusivities” for anions. These “diffusivities” will, moreover, have the opposite dependence on background concentration: the cation flux diverges in the low background concentration limit,⁴ while the anion flux approaches zero.

But this behavior is seen to be caused by differently induced internal concentration gradients. If fluxes are related to these gradients — which they of course should, if you strive for an actual Fickian description — you find that the diffusivities are no different from what is evaluated in closed-cell tests. Relating the steady-state flux to the external concentration difference in the homogeneous mixture model gives (assuming zero tracer concentration on the outflow side)

\begin{equation*} j_\mathrm{ss} = -\phi\cdot D_c \cdot \nabla c^\mathrm{int} = \phi\cdot D_c \cdot\Xi\cdot \frac{c^\mathrm{source}}{L} \end{equation*}

where $c^\mathrm{source}$ denotes the tracer concentration in the external solution on the inflow side, $\phi$ is the porosity, $D_c$ is the pore diffusivity in the interlayer domain, and $L$ is the length of the bentonite sample. From the above equation can directly be identified

\begin{equation} D_e = \phi\cdot\Xi\cdot D_c \end{equation}

$D_e$ is thus not a diffusion coefficient, but basically a measure of $\Xi$.

Note that this explanation for the behavior of $D_e$ does not invoke any notion of an anion accessible volume, nor any “sorption” concept for cations.⁵

Additional comments

When I first published on Donnan equilibrium in bentonite, I was a bit confused and singled out the term “Donnan equilibrium” to refer to anions only, while calling the corresponding cation equilibrium “ion-exchange equilibrium”. To refer to “both” types of equilibrium we used the term “ion equilibrium”.⁶ Of course, Donnan equilibrium applies to ions of any charge and, being better informed, I should have used a more stringent terminology. In later publications I have tried to make amends by pointing out that the process of cation exchange is part of the establishment of Donnan equilibrium.

Being new to the Donnan equilibrium world, I also invented some of my own nomenclature and symbols: e.g. I named the ratio between internal and external concentration the ion equilibrium coefficient ($\Xi$). Conventionally, if I now have understood correctly, this concentration ratio is referred to as the “Donnan ratio”, and is usually labeled $r$ (although I’ve also seen $K$).

But the term “Donnan ratio” seems to be used slightly differently in different contexts, e.g. defined either as $c^\mathrm{int}/c^\mathrm{ext}$ or as $c^\mathrm{ext}/c^\mathrm{int}$, and is sometimes related more directly to the Donnan potential (if no distinction is made between activities and concentrations, we can write $f_D^{-z_i} = c_i^\mathrm{int}/c_i^\mathrm{ext}$). I therefore will continue to use the term “ion equilibrium coefficient” — with label $\Xi$ — in the context of bentonite systems. This usage has also been picked up in some other clay publications. The ion equilibrium coefficient should be understood as strictly defined as $\Xi = c^\mathrm{int}/c^\mathrm{ext}$ for any species, and never to define, or being defined by, the Donnan potential.

To emphasize the difference between effects due to the presence of a Donnan potential and effects due to different local chemical environments, I will refer to $f_D$ as the Donnan factor. (This term does not seem to be used conventionally for any other quantity, although there are examples where it is used as a synonym for Donnan ratio.)

Finally, as in any other approach, the current framework requires a description for the activity coefficients. For activity coefficients in the external solution, there are quite a number of models already available. For the interlayer, modeling — and measuring! — activities is an open research area (at least I hope that this research area is open).

Footnotes

[1] This is just one of several major “loose ends” in most multi-porous models. I have earlier discussed the lack of treatment of swelling, and the incorrect treatment of fluxes in different domains. Update (220622): The lack of a semi-permeable component in multi-porosity models is further discussed here.

[2] This does not have to be the case in principle, if $\Gamma$ for the anion is large, at the same time as the external concentration is not too low.

[3] Hence, it is implied that we use concentration units based on water mass (molality).

[4] What actually happens is that the transport resistance in the filters begins to dominate.

[5] Speaking of “sorption”, we have noted before that this term nowadays is used to mean any type of uptake between bulk water and some other domain (where the species may or may not be immobile). In this sense, there is “sorption” in the homogeneous mixture model (for both cations and anions), but only at interfaces to external solutions. It thus translates to a boundary condition, rather than being part of the transport dynamics within the clay (which makes life much simpler from a numeric perspective). Update (220622): The homogeneous mixture model is extended to deal with ions that truly sorbs here.

[6] It turns out Donnan himself actually used this terminology (“ionic equilibria”)

Sorption part II: Letting go of the bulk water

Leave a reply

Disclaimer: The following discussion applies fully to ions that only interact with bentonite by means of being part of an electric double layer. Here such ions are called “simple” ions. Species with more specific chemical interactions will be discussed in separate blog posts.

The “surface diffusion” model is not suitable for compacted bentonite

In the previous post on sorption¹ we derived a correct “surface diffusion” model. The equation describing the concentration evolution in such a model is a real Fick’s second law, meaning that it only contains the actual diffusion coefficient (apart from the concentration itself)

\begin{equation}
\frac{\partial c}{\partial t} = D_\mathrm{sd} \cdot\nabla^2 c \tag{1}
\end{equation}

Note that $c$ in this equation still denotes the concentration in the presumed bulk water,² while $D_\mathrm{sd}$ relates to the mobility, on the macroscopic scale, of a diffusing species in a system consisting of both bulk water and surfaces.³

Conceptually, eq. 1 states that there is no sorption in a surface diffusion model, in the sense that species do not get immobilized. Still, the concept of sorption is frequently used in the context of surface diffusion, giving rise to phrases such as “How Mobile Are Sorbed Cations in Clays and Clay Rocks?”. The term “sorption” has evidently shifted from referring to an immobilization process, to only mean the uptake of species from a bulk water domain to some other domain (where the species may or may not be mobile). In turn, the role of the parameter $K_d$ is completely shifted: in the traditional model it quantifies retardation of the diffusive flux, while in a surface diffusion model it quantifies enhancement of the flux (in a certain sense).

A correct⁴ surface diffusion model resolves several of the inconsistencies experienced when applying the traditional diffusion-sorption model to cation diffusion in bentonite. In particular, the parameter referred to as $D_e$ may grow indefinitely without violating physics (because it is no longer a real diffusion coefficient), and the insensitivity of $D_\mathrm{sd}$ to $K_d$ may be understood because $D_\mathrm{sd}$ is the real diffusion coefficient (it is not an “apparent” diffusivity, which is expected to be influenced by a varying amount of immobilization).

Still, a surface diffusion model is not a very satisfying description of bentonite, because it assumes the entire pore volume to be bulk water. To me, it seems absurd to base a bentonite model on bulk water, as the most characteristic phenomenon in this material — swelling — relies on it not being in equilibrium with a bulk water solution (at the same pressure). It is also understood that the “surfaces” in a surface diffusion model correspond to montmorillonite interlayer spaces — here defined as the regions where the exchangeable ions reside⁵ — which are known to dominate the pore volume in any relevant system.

Indeed, assuming that diffusion occurs both in bulk water and on surfaces, it is expected that $D_\mathrm{sd}$ actually should vary significantly with background concentration, because a diffusing ion is then assumed to spend considerably different times in the two domains, depending on the value of $K_d$.⁶

Using the sodium diffusion data from Tachi and Yotsuji (2014) as an example, $\rho\cdot K_d$ varies from $\sim 70$ to $\sim 1$, when the background concentration (NaCl) is varied from 0.01 M to 0.5 M (at constant dry density $\rho=800\;\mathrm{kg/m^3}$). Interpreting this in terms of a surface diffusion model, a tracer is supposed to spend about 1% of the time in the bulk water phase when the background concentration is 0.01 M, and about 41% of the time there when the background concentration is 0.5 M⁷. But the evaluated values of $D_\mathrm{sd}$ (referred to as “$D_a$” in Tachi and Yotsuji (2014)) show a variation less than a factor 2 over the same background concentration range.

Insignificant dependence of $D_\mathrm{sd}$ on background concentration is found generally in the literature data, as seen here (data sources: 1, 2, 3, 4, 5)

Diffusion coefficients as a function of background concentration for Sr, Cs, and Na.

These plots show the deviation from the average of the macroscopically observed diffusion coefficients ($D_\mathrm{macr.}$). These diffusion coefficients are most often reported and interpreted as “$D_a$”, but it should be clear from the above discussion that they equally well can be interpreted as $D_\mathrm{sd}$. The plots thus show the variation of $D_\mathrm{sd}$, in test series where $D_\mathrm{sd}$ (reported as “$D_a$”) has been evaluated as a function of background concentration.⁸ The variation is seen to be small in all cases, and the data show no systematic dependencies on e.g. type of ion or density (i.e., at this level of accuracy, the variation is to be regarded as scatter).

The fact that $D_\mathrm{sd}$ basically is independent of background concentration strongly suggests that diffusion only occurs in a single domain, which by necessity must be interlayer pores. This conclusion is also fully in line with the basic observation that interlayer pores dominate in any relevant system.

Diffusion in the homogeneous mixture model

A more conceptually satisfying basis for describing compacted bentonite is thus a model that assigns all pore volume to the surface regions and discards the bulk water domain. I call this the homogeneous mixture model. In its simplest version, diffusive fluxes in the homogeneous mixture model is described by the familiar Fickian expression

\begin{equation} j = -\phi\cdot D_c \cdot \nabla c^\mathrm{int} \tag{2} \end{equation}

where the concentration of the species under consideration, $c^\mathrm{int}$, is indexed with an “int”, to remind us that it refers to the concentration in interlayer pores. The corresponding diffusion coefficient is labeled $D_c$. Notice that $c^\mathrm{int}$ and $D_c$ refer to macroscopic, averaged quantities; consequently, $D_c$ should be associated with the empirical quantity $D_\mathrm{macr.}$ (i.e. what we interpreted as $D_\mathrm{sd}$ in the previous section, and what many unfortunately interpret as $D_a$) — $D_c$ is not the short scale diffusivity within an interlayer.

For species that only “interact” with the bentonite by means of being part of an electric double layer (“simple” ions), diffusion is the only process that alters concentration, and the continuity equation has the simplest possible form

\begin{equation} \frac{\partial n}{\partial t} + \nabla\cdot j = 0 \end{equation}

Here $n$ is the total amount of diffusing species per volume porous system, i.e. $n = \phi c^\mathrm{int}$. Inserting the expression for the flux in the continuity equation we get

\begin{equation} \frac{\partial c^\mathrm{int}}{\partial t} = D_c \cdot \nabla^2 c^\mathrm{int} \tag{3} \end{equation}

Eqs. 2 and 3 describe diffusion, at the Fickian level, in the homogeneous mixture model for “simple” ions. They are identical in form to the Fickian description in a conventional porous system; the only “exotic” aspect of the present description is that it applies to interlayer concentrations ($c^\mathrm{int}$), and more work is needed in order to apply it to cases involving external solutions.

But for isolated systems, e.g. closed-cell diffusion tests, eq. 3 can be applied directly. It is also clear that it will reproduce the results of such tests, as the concentration evolution is known to obey an equation of this form (Fick’s second law).

Model comparison

We have now considered three different models — the traditional diffusion-sorption model, the (correct) surface diffusion model, and the homogeneous mixture model — which all can be fitted to closed-cell diffusion data, as exemplified here

three models fitted to diffusion data for Sr from Sato et al. (92)

The experimental data in this plot (from Sato et al. (1992)) represent the typical behavior of simple ions in compacted bentonite. The plot shows the resulting concentration profile in a Na-montmorillonite sample of density 600 $\mathrm{kg/m^3}$, where an initial planar source of strontium, embedded in the middle of the sample, has diffused for 7 days. Also plotted are the identical results from fitting the three models to the data (the diffusion coefficient and the concentration at 0 mm were used as fitting parameters in all three models).

From the successful fitting of all the models it is clear that bentonite diffusion data alone does not provide much information for discriminating between concepts — any model which provides a governing equation of the form of Fick’s second law will fit the data. Instead, let us describe what a successful fit of each model implies conceptually

The traditional diffusion-sorption model
The entire pore volume is filled with bulk water, in contradiction with the observation that bentonite is dominated by interlayer pores. In the bulk water strontium diffuse at an unphysically high rate. The evolution of the total ion concentration is retarded because most ions sorb onto surface regions (which have zero volume) where they become immobilized.
The “surface diffusion” model
The entire pore volume is filled with bulk water, in contradiction with the observation that bentonite is dominated by interlayer pores. In the bulk water strontium diffuse at a reasonable rate. Most of the strontium, however, is distributed in the surface regions (which have zero volume), where it also diffuse. The overall diffusivity is a complex function of the diffusivities in each separate domain (bulk and surface), and of how the ion distributes between these domains.
The homogeneous mixture model
The entire pore volume consists of interlayers, in line with the observation that bentonite is dominated by such pores. In the interlayers strontium diffuse at a reasonable rate.

From these descriptions it should be obvious that the homogeneous mixture model is the more reasonable one — it is both compatible with simple observations of the pore structure and mathematically considerably less complex as compared with the others.

The following table summarizes the mathematical complexity of the models ($D_p$, $D_s$ and $D_c$ denote single domain diffusivities, $\rho$ is dry density, and $\phi$ porosity)

Note that the simplicity of the homogeneous mixture model is achieved by disregarding any bulk water phase: only with bulk water absent is it possible to describe experimental data as pure diffusion in a single domain. This process — pure diffusion in a single domain — is also suggested by the observed insensitivity of diffusivity to background concentration.

These results imply that “sorption” is not a valid concept for simple cations in compacted bentonite, regardless of whether this is supposed to be an immobilization mechanism, or if it is supposed to be a mechanism for uptake of ions from a bulk water to a surface domain. For these types of ions, closed-cell tests measure real (not “apparent”) diffusion coefficients, which should be interpreted as interlayer pore diffusivities ($D_c$).

Footnotes

[1] Well, the subject was rather on “sorption” (with quotes), the point being that “sorbed” ions are not immobilized.

[2] Eq. 1 can be transformed to an equation for the “total” concentration by multiplying both sides by $\left (\phi + \rho\cdot K_d\right)$.

[3] Unfortunately, I called this quantity $D_\mathrm{macr.}$ in the previous post. As I here compare several different diffusion models, it is important to separate between model parameters and empirical parameters, and the diffusion coefficient in the “surface diffusion” model will henceforth be called $D_\mathrm{sd}$. $D_\mathrm{macr.}$ is used to label the empirically observed diffusion parameter. Since the “surface diffusion” model can be successfully fitted to experimental diffusion data, the value of the two parameters will, in the end, be the same. This doesn’t mean that the distinction between the parameters is unimportant. On the contrary, failing to separate between $D_\mathrm{macr.}$ and the model parameter $D_\mathrm{a}$ has led large parts of the bentonite research community to assume $D_\mathrm{a}$ is a measured quantity.

[4] It might seem silly to point out that the model should be “correct”, but the model which actually is referred to as the surface diffusion model in the literature is incorrect, because it assumes that diffusive fluxes in different domains can be added.

[5] There is a common alternative, implicit, and absurd definition of interlayer, based on the stack view, which I intend to discuss in a future blog post. Update (220906): This interlayer definition and stacks are discussed here.

[6] Note that, although $D_\mathrm{sd}$ is not given simply by a weighted sum of individual domain diffusivities in the surface diffusion model, it is some crazy function of the ion mobilities in the two domains.

[7] With this interpretation, the fraction of bulk water ions is given by $\frac{\phi}{\phi+\rho K_d}$.

[8] The plot may give the impression that such data is vast, but these are basically all studies found in the bentonite literature, where background concentration has been varied systematically. Several of these use “raw” bentonite (“MX-80”), which contains soluble minerals. Therefore, unless this complication is identified and dealt with (which it isn’t), the background concentration may not reflect the internal chemistry of the samples, i.e. the sample and the external solution may not be in full chemical equilibrium. Also, a majority of the studies concern through-diffusion, where filters are known to interfere at low ionic strength, and consequently increase the uncertainty of the evaluated parameters. The “optimal” tests for investigating the behavior of $D_\mathrm{macr.}$ with varying background concentrations are closed-cell tests on purified montmorillonite. There are only two such tests reported (Kozaki et al. (2008) and Tachi and Yotsuji (2014)), and both are performed on quite low density samples.

$D_a$ is not “apparent” (sorption, part I)

1 Reply

Disclaimer: The following discussion applies fully to ions that only interact with bentonite by means of being part of an electric double layer. Species with more specific chemical interactions will be discussed in separate blog posts.

The traditional diffusion-sorption model

A persistent idea, that has been around since the very start of research on compacted bentonite¹, is that “sorbed” ions are immobilized, and early diffusion models were based on two assumptions

The entire pore volume contains bulk water, where ions diffuse according to \begin{equation} j = -D_e \nabla c \tag{1} \end{equation} where $c$ is the ion concentration and $D_e = \phi D_p$, with $\phi$ being the porosity, and $D_p$ a so-called pore diffusivity (containing information on the ion mobility at the macroscopic scale of the porous system).
Ions in the porewater “sorb” on solid surfaces and become immobilized. The amount of “sorbed” ions per unit solid mass, $s$, is related to $c$ via a distribution coefficient $K_d$ \begin{equation} s = K_d\cdot c \tag{2} \end{equation}

With these assumptions, the total amount of ions in the bentonite (per volume porous medium) is

\begin{equation} n = \phi\cdot c + \rho\cdot s = \left (\phi+\rho K_d \right) c \tag{3} \end{equation}

where $\rho$ is the (dry) density of the clay. Applying the continuity equation, $\partial n/\partial t = -\nabla\cdot j$, gives

\begin{equation} \frac{\partial c}{\partial t} = \frac{D_e}{\phi + \rho K_d}\nabla^2c \tag{4} \end{equation}

This equation is usually written in terms of the so-called apparent diffusion coefficient $D_a$

\begin{equation} \frac{\partial c}{\partial t} = D_a \nabla^2c \tag{5} \end{equation}

with

\begin{equation} D_a = \frac{D_e}{\phi + \rho K_d} \tag{6} \end{equation}

The model defined by eqs. 5 and 6 — which we will refer to as the traditional model — has a distinct physical interpretation: ions diffuse relatively fast in a bulk water phase, while being retarded due to “sorption” onto the solid surfaces. Although eq. 5 has the form of Fick’s second law, it is clear that $D_a$ is not a real diffusion coefficient, but is influenced both by diffusion ($D_e$) and “sorption” ($K_d$).

The immobilization assumption is not valid

The traditional model has a major problem: it is not valid for compacted bentonite. There is — and was — strong evidence that the immobilization assumption does not hold for all types of “sorbing” ions. E.g Jenny and Overstreet (1939) reported transport of $\mathrm{Fe^{2+}}$, $\mathrm{H^+}$, $\mathrm{K^+}$, $\mathrm{Na^+}$, and $\mathrm{Ca^{2+}}$ as samples of H- K- Na- and Ca-bentonite was brought in contact with samples of Fe/H-bentonite (without gel contact there was no transport), and van Schaik et al. (1966) demonstrated steady-state fluxes of $\mathrm{Na^+}$ through Na-montmorillonite, corresponding to pore diffusivites significantly larger than the $\mathrm{Na^+}$ diffusivity in pure water. Moreover, ion mobility is an essential component of the most well-established theoretical concept used for describing montmorillonite — the electric double layer. The ionic part of the electric double layer is even referred to as the diffuse layer, for crying out loud!²

So it should come as no surprise that the pioneers in radwaste barrier research quickly came to the conclusion that the description underlying eq. 6 does not hold. Several research groups reported values of $K_d$ and $D_a$ that gives unreasonably large values of $D_e$ when evaluated using eq. 6 (even larger than the corresponding diffusivities in pure water, in several cases).

What is more surprising is that these conclusions did not lead to a complete reevaluation of the underlying assumptions, which is the preferred procedure for a model producing nonsense. Instead, the invalidity of the traditional model was only considered a problem for the specific ions for which it yield “unrealistic” parameter values (in particular, strontium and cesium), and an adopted “remedy” was the so-called surface diffusion model. The surface diffusion model keeps the description expressed by eqs. 1, 2, 5, and 6, while replacing the relation $D_e = \phi\cdot D_p$ with the following monster

\begin{equation} D_e = \phi\cdot D_p + \rho\cdot K_d\cdot D_s \tag{7} \end{equation}

where $D_s$ is a so-called surface diffusion coefficient.

This extension of the traditional model brings about several problems, the most glaring one being that eq. 7 is wrong. Therefore, to be able to continue the discussion, we first have to derive a correct “surface diffusion”³ model.

A correct “surface diffusion” model

Eq. 7 is not valid because it relies on the incorrect assumption that diffusive fluxes from different domains are additive. To derive the correct equations, we go back to the only conceptual modification made in the surface diffusion model: that the “sorbed” ions are no longer assumed immobile. Thus, the basic assumptions of the surface diffusion model are:

The entire pore volume contains bulk water, where ions diffuse.
Ions in the porewater “sorb” on the solid surfaces.
“Sorbed” ions also diffuse, along the surfaces.

There are consequently two domains (bulk water and surfaces) where diffusion is assumed to take place, and a correct equation for the macroscopic flux is in this case (for details, see here)

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

where $D_\mathrm{macr.}$ is the diffusion coefficient on the macroscopic scale, which conveys information on the mobility in both of the involved domains, as well as their geometrical configuration. $\bar{c}$ is the average concentration of all mobile entities; here, this include all ions — those in the porewater as well as those on the surfaces. The bulk water concentration is still labelled $c$. Using eq. 3, we can write

\begin{equation} \phi\bar{c} = n = \left (\phi + \rho K_d \right) c \tag{9} \end{equation}

and eq. 8 reduces to

\begin{equation} j = – D_\mathrm{macr.}\left (\phi + \rho K_d \right) \nabla c \tag{10} \end{equation}

Comparing with eq. 1 directly shows

\begin{equation} D_e = D_\mathrm{macr.}\left (\phi + \rho K_d \right) \tag{11} \end{equation}

Further, when using eq. 10 in the continuity equation, the factor $\left (\phi + \rho K_d \right)$ cancels, giving

\begin{equation} \frac{\partial c}{\partial t} = D_\mathrm{macr.} \nabla^2 c \tag{12} \end{equation}

A complete conceptual change

Eqs. 11 and 12 are the main result of a correctly derived “surface diffusion” model. Comparing with eqs. 5 and 6, we see that, although the equations are deceptively similar, a “surface diffusion” model brings about a complete conceptual change: with the immobilization assumption released, it is actually the real diffusion coefficient that appear in Fick’s second law (eq. 12). In fact, the surface diffusion model does not contain any apparent diffusivity!

Similarly, $D_e$ is not a real diffusion coefficient in the surface diffusion model, as clearly seen by eq. 11. Rather, $D_e$ is the product of the diffusion coefficient and a concentration factor, stemming from the fact that the general flux equation reads $j = -\frac{D}{RT} \bar{c}\cdot \nabla \mu$, where $\mu$ is the chemical potential of the diffusing ion.

In the same vein, $K_d$ does not quantify “sorption” — at least if “sorption” is supposed to refer to a process that retards the diffusive flux. On the contrary, $K_d$ acts as a kind of amplifier of the flux (eq. 10). This implication was already pointed out by Jahnke and Radke (1985):

Thus, $K_d$ values are of extreme importance during the steady release period. Counter to what may have been anticipated, larger $K_d$ values lead to higher release rates.

The conceptual changes of a “surface diffusion” model make a lot of sense when interpreting diffusion data:

$D_e$ is expected to vary with $K_d$ and can basically grow without limit, without becoming unphysical.
“$D_a$” for a whole bunch of cations is universally found to be insensitive to changes in $K_d$, even though $K_d$ may vary orders of magnitudes (by altering the background concentration). This is the expected behavior of a real diffusion coefficient, and it makes all sense to instead interpret this parameter as $D_\mathrm{macr.}$, using eq. 12.
When assuming all ions in the system to be mobile, we expect the governing equation (eq. 12) to be a real Fick’s second law, i.e. that it contains nothing but the actual diffusion coefficient.

Despite this sense-making, I have never seen these points being spelled out in the bentonite literature. Rather, it appears as the implication of releasing the immobilization assumption has not been fully comprehended. To this day, the traditional model is used as the basis for interpretation in basically all publications, to the extent that “$D_a$”, “$D_e$”, and “$K_d$” usually are considered to be experimental quantities. Similarly, most modern model development is focused on extending the traditional model, using flawed concepts regarding diffusivity in multi-porous systems.

I also find it symptomatic that “surface diffusion” not uncommonly is treated as being an optional mechanism, or even being completely discarded, while the dominant contribution of “sorbed” ions to the diffusive flux seems to have to be rediscovered every 25^th year, or so.

Note that the conclusions made here are not based on any model that I favor, but are the consequence of a correct conceptual treatment of what is stated in the contemporary bentonite literature. Although a “surface diffusion” model (correctly treated!) make much more sense as compared with the traditional model in interpreting cation tracer diffusion, it still involves an assumption that I don’t understand how anybody can endorse: it assumes the entire pore volume to consist of bulk water. I will discuss this point in a separate blog post. Update (210225): Letting go of the bulk water is discussed here.

Footnotes

[1] In the following I will write “bentonite”, but I mean any type of clay system with significant ion exchange properties.

[2] I guess a reason for not acknowledging interlayer diffusivity could have been the incorrect notions that interlayer water is crystalline, and that diffuse layers only “form” at higher water content.

[3] Since the established surface diffusion model is that expressed by eq. 7, I use quotation marks when speaking about models that don’t comply with this equation.

You cannot add fluxes from different domains in multiporosity models

2 Replies

Most models of compacted bentonite and claystone¹ 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.² 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,³ 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).

The through-diffusion set-up

The effective porosity diffusion model

Connection with experimentally accessible quantities

Additional comments on the effective porosity model

Footnotes

Muurinen et al. (1985)

Eriksen et al. (1999)

Gimmi & Kosakowski (2011)

Some thoughts

Footnotes

Matric suction (capillary suction)

Osmotic suction and total suction

Suction in bentonite

Footnotes

Overview

Chemical equilibrium

Illustrations

Equation for \(f_D\)

Comment on through-diffusion

Additional comments

Footnotes

The “surface diffusion” model is not suitable for compacted bentonite

Diffusion in the homogeneous mixture model

Model comparison

Footnotes

The traditional diffusion-sorption model

The immobilization assumption is not valid

A correct “surface diffusion” model

A complete conceptual change

Footnotes

Equation 1 does not make sense

Some illustrations

A more reasonable formulation

Footnotes