Category Archives: Multicomponent diffusion

Post-publication review: Tournassat and Steefel (2015), part VI

This is the sixth and final (!) part of the review of “Ionic Transport in Nano-Porous Clays with Consideration of Electrostatic Effects” (Tournassat and Steefel (2015) (referred to as TS15 in the following). For background and context, please check the first part.

Recap

Let’s begin by reminding ourselves where we are in the “story” of the article. TS15 begun by vaguely suggesting that diffusion in bentonite1 cannot be modeled by assuming “Fickian” diffusion. The argument, however, is based on the inability of the traditional diffusion-sorption model to be reasonably fitted to ion through-diffusion data. This inability is not primarily due to assuming Fickian diffusion, but due to assuming immobility of the “sorbed” ions and due to assuming that the pore volume contains nothing but bulk water.

After being promised to be presented “fundamental properties” of bentonite that apparently will demonstrate the problem with assuming “Fickian” diffusion, we were instead introduced to a large set of modeling concepts. As a description of actual compacted bentonite, these concepts are mostly irrelevant (e.g. the Gouy-Chapman model), unjustified (e.g. Stern layers on basal surfaces), or simply fantasies (e.g. the idea of “stacks”, or the idea of Donnan equilibrium on the microscopic scale).

Finally, we were presented the Nernst-Planck approach to ion diffusion, which indeed goes beyond a Fickian description. The presented model, however, which includes “bulk water”, “diffuse layer” and “interlayer” domains assumed to be in equilibrium locally, is not suitable for compacted bentonite. Furthermore, no description whatsoever is presented for how this equilibrium is supposed to be maintained. The model is therefore incapable, even principally, to account for swelling — the most profound feature of bentonite in contact with an aqueous solution. Note also that discrimination between “diffuse layer” and “interlayer” domains is based on the pure fantasy of “stacks” as fundamental structural units.

We also showed that the presented Nernst-Planck framework was incorrectly derived, and after having sorted out how it can be reinterpreted,2 we noted that the suggested “concentration” and “diffusion potential” contributions to the diffuse layer flux where completely unfounded.

“From diffusive flux to diffusive transport equations”

In this section I expected to be given an overview of how all the previously introduced model components are supposed to fit together in an overall description for bentonite. With so many parts, there are a large number of quantities to consider. Apart from all the species concentrations in the different “porosity domains”, TS15 have also introduced e.g. electrostatic potentials (eqs. 17, 18, 21, 22, 39, 48, 56 in TS15), adsorption sites on “basal surfaces” (eq. 19 in TS15), “exchanger” sites (eq. 24 in TS15), “layer edge amphoteric” sites (eq. 28 in TS15), models for activity coefficients and their spatial derivatives (eqs. 58, 60, 64 in TS15), and entertained the possibility of having varying sizes of the “porosity domains” (although no constitutive relations has been stated).

Instead, the only thing related to a model description here is a restatement of the mass conservation law. This was already presented in eq. 4 in TS15, but is here incorrectly called the “diffusion equation” (TS15 eq. 65), and stated with an incorrect sign (TS15 eqs. 65, 66 and 68). TS15 express it in its final form as (I have here corrected the sign error, and suppressed an index \(i\))

\begin{equation} \frac{\partial}{\partial t} \left ( \phi_\mathrm{bulk}c_\mathrm{bulk} + \phi_\mathrm{DL}c_\mathrm{DL} + \phi_\mathrm{exch}c_\mathrm{exch} \right) = -\frac{\partial}{\partial x} \left ( j_\mathrm{bulk} + j_\mathrm{DL} + j_\mathrm{exch} \right ) \tag{1} \end{equation}

where \(\phi_\mathrm{bulk}\), \(\phi_\mathrm{DL}\) and \(\phi_\mathrm{exch}\) denotes the bulk water, “diffuse layer” and “interlayer” porosites, respectively; \(c_\mathrm{bulk}\), \(c_\mathrm{DL}\), \(c_\mathrm{exch}\) are the concentrations of the considered species in these “porosity domains”; and \(j_\mathrm{bulk}\), \(j_\mathrm{DL}\), and \(j_\mathrm{exch}\) are the corresponding domain fluxes (presumably given from the previously presented Nernst-Planck framework?).

That’s it.

As the authors stress the importance of adsorption processes throughout the article (even on montmorillonite basal surfaces), I find it very odd to here present an equation that completely omits them. Somewhat ironically, adsorption was included previously, in the section titled “Diffusion basics” (eq. 6 — 8 in TS15). Instead, the above equation (eq. 1) implies mechanisms that may alter the sizes of the “porosity domains” (i.e. \(d\phi/dt \neq 0\), \( d\phi/dx \neq 0 \)). But the article has not presented any constitutive relations for such mechanisms! This whole model presentation feels like having been handed a heap of irregular puzzle pieces with a promise of a grand result, while the reference photo belongs to a different puzzle.

Although eq. 1 provides little concrete information on the actual model, it still manages to contain a fatal flaw: the three domain fluxes are simply added to form some sort of a “total” flux. But this is not how diffusive fluxes work! The first real post I wrote on this blog covered this issue. There we noticed that domain fluxes are added in essentially all publications where multi-porous models are applied to bentonite. To see that separate domain fluxes are not additive should be as simple as considering some specific domain configuration. A trivial example is two domain types connected in series

It is easy to see that the resulting flux in this case cannot be described as a sum of independent domain contributions (e.g, the flux vanishes if \(D_1 = 0\), independent on the value of \(D_2\)). More generally, we may instead refer to one of the favorite references on diffusion in the bentonite research field: Crank (1975) has an entire chapter on this topic (“Diffusion in heterogeneous media”).

It really feels silly to point out, but the knowledge that domain diffusivites generally don’t add was established long before any numerical reactive transport tools were developed. I have trouble expressing what an incredible own goal this procedure is: even if it made sense to partition compacted bentonite into various “porosity domains”, using eq. 1 will still lead to incorrect inferences.

“Applications”

Rather than providing the reader with an actual model overview, TS15 proceed with presenting three different cases, which were modeled with either “CrunchFlowMC” or “PHREEQC”. These case studies are motivated with that they

[…] are used to illustrate the importance of considering coupled diffusion/surface reaction effects in order to understand and to predict migration processes and associated parameters in charged porous media, especially clays.

Let’s keep this formulation in mind as we go through the cases. (Spoiler: we will not be illustrated the importance of much, nor will we be able to understand or predict anything.)

The three tests that have been modeled are:

  • Tachi and Yotsuji (2014)

    This is a diffusion and “sorption” study of cesium, sodium and iodide in pure Na-montmorillonite. Tracer through-diffusion has been performed for each of these ions in samples of 0.8 g/cm3, at three different background concentrations (NaCl): 0.01 M, 0.1 M, and 0.5 M.

  • Glaus et al. (2013)

    This is the “uphill” diffusion study. It presents the result of two through diffusion tests in pure Na-montmorillonite (1.9 g/cm3), where a Na tracer is observed to migrate from a 0.1 M NaClO4 reservoir to a 0.5 M NaClO4 reservoir, although the tracer concentration in the 0.1 M reservoir is lower than in the 0.5 M reservoir.

  • Tertre et al. (2015)

    This paper presents the results of performing Ca-to-Na ion exchange on initially pure Ca-vermiculite (1.8 g/cm3). Disc-shaped single crystal vermiculite samples were immersed in initially pure NaCl solutions, and the evolution of calcium concentration in these solutions was monitored. The samples were prepared in such a way that out-diffusion only occurred at the outer edge of the vermiculite discs. Tertre et al. (2015) present four such tests, performed with NaCl concentrations 0.003 M, 0.05 M, 0.1 M, and 1.0 M. In addition, three tests were performed with NaCl concentration 0.05 M, which were terminated after 1.94, 7.01, and 18.06 days, respectively. Subsequently, the radial distribution of calcium was measured in these samples, giving “snapshots” of the internal Ca profile.

A common theme of these studies is that they primarily demonstrate effects of exchangeable cations being mobile. Tachi and Yotsuji (2014) clearly show the commonly observed effect of cation “effective diffusivity” becoming increasingly large as the background concentration is lowered, which has been unambiguously demonstrated to be an effect of “interlayer diffusion being the dominant pathway”. Glaus et al. (2013) is a proof-of-concept demonstration of a consequence of the exchangeable cations dominating mass transfer, as I discussed in some detail here. Tertre et al. (2015) assert from the start that “For this type of model system, the interlayer region is the only pore space that contributes to the overall diffusion process […]”.

With this background information, we will go through how TS15 handle and present these modeling cases. We will do this with focus on what I think are some reasonable requirements on a modeling exercise, especially in a review article that is supposed be a “fully developed text which can be used for self-study, research, or as a text-book for graduate-level courses.” According to me, a modeling task should include

  • A clear statement of the exact model used
  • Some sort of quality analysis of the data being modeled, e.g. a separation of signal and noise. It should also be motivated why certain data is modeled while other is not.
  • A summary that relates the model results to actual physical processes, and discusses possible implications.

Tachi and Yotsuji (2014)

Without commenting on it, TS15 only model results obtained at background concentration 0.1 M and only in samples of length 10 mm. A reason for this may be that Tachi and Yotsuji (2014) only present reservoir concentration evolutions and final state concentration profiles for all ions at these conditions. However, TS15 anyway disregard the final state profiles and only compare their model with the concentration evolutions. Furthermore, for cesium, the original article contains plots also for background concentrations 0.01 M and 0.5 M. Tachi and Yotsuji (2014) actually present fitted model parameters for all tests. If TS15 were serious about exploring this type of diffusion they could consequently have made a systematic comparison rather than just fitting a model to some rather arbitrarily chosen test results. This is particularly remarkable, as the most prominent feature of tests like this — that the “effective diffusion coefficient” for cations increases essentially without bound with decreasing background concentration, while the same quantity for anions approaches zero — is not really explored, when they only model results from a single background concentration.

It is also clear that TS15 show no interest in the quality of the data. Looking at Tachi and Yotsuji (2014) in a bit more detail, the data has some flaws and inconsistencies. For example, for all background concentrations, Tachi and Yotsuji (2014) report two quite different values of the “anion accessible porosity”,3 depending on whether the sample is 5 or 10 mm long. (TS15 “conveniently” choose the value 0.421 and ignore the value 0.635.) A further inspection also indicates that filter restriction has influenced some of the data, although they have attempted to eliminate such artifacts.4

TS15 don’t clearly state what model components they include, and it is essentially impossible to figure out the exact model used for this exercise. Half way through the section, the model is referred to as the “bulk + diffuse layer water diffusion model”, so it is quite safe to state that they have only included a bulk water and a “diffuse layer” domain. As the description in TS15 is fundamentally based on the (flawed) idea that we can distinguish between “external” and “internal” surfaces in compacted bentonite, it is quite remarkable that the “interlayer” domain is here completely left out. To me, this suggests that TS15 are not themselves especially serious about differing between “diffuse layers” and “interlayers” (after all, this difference is a fantasy, as I have argued throughout the review). This sloppy approach to modeling is, in my opinion, a recipe for ending up with an overparameterized mess. It is also worth noting that TS15 simply fix the size of the the bulk water domain by requiring the model to be “in close agreement” with measured “anion accessible porosity”.5 Thus, after 34 pages of text containing a myriad of mathematical expressions and a specific “consideration of electrostatic effects,” we ultimately just end up with knob-twiddling.

Speaking of electrostatic effects, TS15 don’t comment how aspects of multi-component diffusion influence the results. It is not even clear that the Nernst-Planck framework has been used in this particular modeling. For example, only a single set of fitted tracer diffusion coefficients are reported6 (and no transport parameters are reported for the main electrolyte). Are these supposed to relate to the bulk or the “diffuse layer” domain? I find this lack of discussion quite remarkable, given the emphasis in the rest of the article on the idea that ion diffusion in bentonite is supposed to be fundamentally “non-Fickian”.

TS15 do not even manage to clearly communicate the principal reason as for why the model “makes it possible to calculate the diffusive flux of neutral species, anions and cations with the same conceptual model and with physically feasible parameters.”7 This reason is — of course — that the model allows for exchangeable ions to move and that it provides a mechanism (Donnan equilibrium) for altering the relative ion concentrations in the “diffuse layer” domain. In fact, a more reasonable analysis demonstrates that tracer through-diffusion in bentonite can be satisfactory modeled using a single homogeneous domain for the clay. A homogeneous description is also compatible with basically all other experimental findings in compacted bentonite, as we have discussed massively on the blog.

In contrast, TS15 assume, not only that the montmorillonite in this particular test contains 50% bulk water, but also the existence of a Stern Layer — with corresponding ion immobilization — on montmorillonite basal surfaces. As neither of these assumptions are justified or reasonable, it is clear that the resulting model is overparameterized. Any successful fit thus provides little real information or insight, which also becomes almost comically clear when it is concluded

Results are plotted in Figure 19. The tortuosity values follow the order \(\tau_\mathrm{I^-} < \tau_\mathrm{HTO} < \tau_\mathrm{Na^+} < \tau_\mathrm{Cs^+}\), indicating that the tortuous diffusion pathways are not the same for all of these species. Or alternatively, that the adsorbed species in the Stern layer, considered in the calculations as immobile, are in fact mobile.

The result of the modeling exercise indicates either one thing, or a completely different thing!

Note that TS15 here put forward the bizarre notion that each separate species diffuse in different pore structures, with different “tortuosities”. They apparently don’t dismiss the untenable idea that cesium (with fitted “tortuosity” 0.136) for some reason diffuses though a much less meandering pore network than water (with “tortuosity” 0.047). In reality, the much larger “tortuosity” (which here means a less tortuous network…) for cesium as compared with water indicates that the model underestimates the amount of mobile cesium in the “diffuse layer” and compensates for this by giving it a too large diffusivity. This is signature behavior of an overparameterized model.

How this treatment and way of reasoning is supposed to “illustrate the importance of considering coupled diffusion/surface reaction effects” should be a mystery to any reasonable reader.

Glaus et al. (2013)

Things get even worse, in my opinion, as TS15 next present a modeling of the “uphill” test. Providing a solid explanation for this test could have been a major contribution to the article, as it earlier has claimed that the “uphill” effect is an example of “non-Fickian diffusion processes” and has implied that the effect is caused by non-trival electrostatic couplings. But TS15 seem genuinely uninterested in analyzing how the “uphill” effect arises in their model. Again, the adopted model is not presented in full, and focus is mainly on irrelevant details, such as filter diffusivity.8 The only attempt at describing a mechanism for “uphill” diffusion is half a sentence that says

[…] the macroscopic driving force of the diffusion is the diffusion potential in the diffuse layer, as shown in Figure 12.

But “Figure 12” shows results of “example 3”, from the section on the Nernst-Planck framework (discussed here), rather than results from the actual “uphill” test modeling. The burden of connecting “Figure 12” with the actual “uphill” test is completely left to the reader. I find this very peculiar for several reasons. To begin with, as TS15 evidently have modeled the actual “uphill” test, why don’t they refer to and show the “diffusion potential in the diffuse layer” in the actual model? And although “Figure 12” shows an “uphill” cation flux, the Nernst-Planck transport equations are not solved in “example 3”; in these examples, flux contributions are calculated from imposed linear concentration profiles in a bulk water domain. But since the “uphill” test here has been modeled with a tool that apparently solves the Nernst-Planck transport equations, I find it incomprehensible that these results are not reported. It is also worth noting that TS15 don’t compare their model with the final-state concentration profiles, even though these are readily available in the original publication.

The earlier examples, moreover, assume a density of approximately 1.0 g/cm3 (porosity 0.64), and a bulk water domain that dominates the pore volume. The “uphill” test, in contrast, was conducted at nominal density 1.9 g/cm3 (porosity \(\sim\) 0.3). TS15 themselves explicitly comment that “[a]t this high degree of compaction, the actual presence of bulk water is not certain”.9 Other conditions also differ strongly between these systems. The “uphill” test was performed with NaClO4 as main electrolyte, while “example 3” uses NaCl. “Example 3” has boundary concentrations 0.1 M and 0.001 M for the main electrolyte, while the actual “uphill” test, in contrast, has 0.5 M and 0.1 M.

Disregarding these glaring differences between the actual “uphill” test and “example 3”, the explanation given by TS15 in the above quotation is that the cation tracer ions move “uphill” due to “the diffusion potential in the diffuse layer”. In an earlier part of this review we demonstrated that the way TS15 define the “diffusion potential” contribution to the flux is obvious nonsense, based on misunderstanding the electric potential. To remind ourselves, in “example 3” TS15 make the absurd claim that a tracer concentration profile in the “diffuse layer” that looks like this (this is part of “Figure 12”)

makes NO “concentration gradient” contribution to the flux (because TS15 define such contributions exclusively in terms of bulk water concentrations). Instead, the diffuse layer cation tracer flux in “example 3” is supposed to be driven by an electric field (the “diffusion potential” contribution) in the wrong direction.

We have shown that these preposterous definitions of flux contributions can be corrected, with the correction term involving the gradient of the Donnan potential.10 With reasonable flux definitions, we demonstrated that the “uphill” flux in “example 3” is essentially completely governed by the contribution from the “diffuse layer” concentration gradient (i.e. that the cation tracer flux is “Fickian”). Recently on the blog, we have also calculated the electrostatic potential within the homogeneous mixture model and specifically demonstrated that effects of electromigration on the tracer flux is negligible in the “uphill” test.

Making a way too generous interpretation of the above statement, it may perhaps be claimed to be “true” if we accept that “the diffusion potential” refers to the absurd definition given in TS15, rather than to an actual electric field. It is, however, clear that TS15 completely fail to clarify that the “uphill” effect is caused by exactly the same mechanism as is active in conventional cation tracer through-diffusion in bentonite. Even worse, the way the “uphill” effect is discussed here may fool a reader to believe that it has other causes.

Rather than identifying any of the above points, TS15 instead spend a large part of this section discussing that several models, which “are totally different from a conceptual view”, can be fitted to the “uphill” diffusion test data. But all models that can be fitted to the “uphill” test data have this ability because they (i) assume mobility of exchangeable ions and (ii) provide a mechanism for ion equilibrium with bulk water. From the point of view of fitting the data, it is secondary if a model treats the ion equilibrium via “exchange sorption” or Donnan equilibrium,11 or if the compartment containing the exchangeable ions is referred to as “interlayer”, “diffuse layer” or interlayer. That TS15 choose to emphasize that their model fits “equally good”, while it contains several additional and unjustified assumptions, simply indicates that they don’t consider the dangers with overparameterization.

The presentation of the modeling of the “uphill” test in TS15 does not “illustrate the importance of considering coupled diffusion/surface reaction effects in order to understand and to predict migration processes…”. Rather, it obscures process understanding.

Tertre et al. (2015)

As with the previous cases, the description of the modeling of Tertre et al. (2015) lacks sufficient detail, preventing the reader from determining exactly which model was used. Actually, the only model information given is that the data was “reinterpreted using the interlayer diffusion option of PHREEQC”. This must obviously mean that an “interlayer” domain is included in the model, in contrast to the previous cases. But are other domains included? (we have previously been informed, in passing, that “PHREEQC” requires a bulk water domain).

Furthermore, it is not clear if a Stern layer is included in this model, i.e. if some of the exchangeable ions are assumed immobilized, as in the previous cases. As not a word is spoken about this, we reasonably must assume that no Stern layer is included. If this is the case, the inclusion of an immobilization mechanism when modeling Tachi and Yotsuji (2014) and Glaus et al. (2013) becomes even more conspicuous. In fact, TS15 state about Tertre et al. (2015)

This experiment showed unambiguously that interlayer diffusion exists and it made it possible also to quantify the diffusion coefficient for \(Ca^{2+}\) in the interlayer.

I find this statement remarkable for a couple of reasons. As TS15 don’t say anything similar when discussing the other studies, I here get the impression that they haven’t comprehended that “interlayer diffusion” has been unambiguously demonstrated much earlier than in Tertre et al. (2015). Also, if TS15 now acknowledge that interlayer diffusion “exists” in bi-hydrated vermiculite, how come that they don’t assume this to be the primary mode of transport in the other models, but instead add an immobilization mechanism? Are ions supposed to become less mobile in less dense systems with lower layer charge density? We should also remember that earlier in the article, TS15 describe the “interlayer” as pure Stern layers. So, how is it, are Stern/interlayer ions supposed to be mobile or not? Again, I think this relaxed attitude towards when and how various domains are included, and how they are treated, demonstrates that TS15 themselves don’t take the concepts of “diffuse layer” and “interlayer” seriously.12

Furthermore, the “reinterpretation” of this experiment is entirely untenable. Tertre et al. (2015) put a lot of effort in demonstrating (using “Brownian dynamics simulations”) that the observed diffusion is strongly influenced by processes occurring at the interface between the sample and the external solution. This makes sense, as we may expect huge concentration gradients within the clay when two different cations are involved. Such gradients are difficult to maintain, and the resulting flux can instead be expected to be controlled by interface transfer resistance. This is similar to how filters limit the flux in cation through-diffusion tests at low ionic strength. In the “reinterpretation”, TS15 completely neglect any possible interface transfer resistance, claiming

It was possible to reproduce the data with the same level of quality as the Brownian dynamics calculations by considering an interlayer \(Ca^{2+}\) diffusion coefficient value of \(0.8\times 10^{-11}\) \(\mathrm{m}^2\cdot \mathrm{s}^{-1}\), together with an interlayer \(Na^{+}\) diffusion coefficient of \(4\times 10^{-11}\) \(\mathrm{m}^2\cdot \mathrm{s}^{-1}\), \(1\times 10^{-11}\) \(\mathrm{m}^2\cdot \mathrm{s}^{-1}\), \(0.1\times 10^{-11}\) \(\mathrm{m}^2\cdot \mathrm{s}^{-1}\), and \(0.1\times 10^{-11}\) \(\mathrm{m}^2\cdot \mathrm{s}^{-1}\), for the experiments at NaCl concentration of 1, 0.1, 0.05 and 0.003 \(\mathrm{mol}\cdot \mathrm{L}^{-1}\) respectively. The decrease of the \(Na^{+}\) interlayer diffusion coefficient with NaCl concentration is in agreement with MD results from Tertre et al. (2015), which show a decrease of its value with an increase of the \(Ca^{2+}\)/\(Na^{+}\) occupancy ratio in the interlayer (Fig. 21c).

The last sentence is false. Tertre et al. (2015) indeed report a decrease of the sodium self-diffusion coefficient, but this decrease is lower than a factor 3.13 In the model of TS15, in contrast, this decrease is a factor 40 (!), which certainly cannot be considered “in agreement”. In addition, experimentally, Kozaki et al.(2005) observed no systematic decrease of Na diffusivity with increasing Ca fraction, in bi-hydrated montmorillonite; if anything, they observe a slight decrease (25 — 60%) of Na diffusivity only in pure Ca-montmorillonite. TS15 anyway goes so far as to conclude this part by stating

Consequently, the apparent decrease in \(Ca^{2+}\) interlayer diffusion coefficient can also be interpreted as a result of the decrease of the \(Na^{+}\) interlayer diffusion coefficient that arises from the coupling between the diffusion of these two species through the diffusion potential term in Equation (57).

Let’s unpack the serious error of reasoning on display here: rather than taking into account interface transfer resistance — a process that we have every reason to believe is active, and that is expected to produce the observed behavior — the observed behavior is supposed to be caused by an unjustified 40-fold decrease of the sodium self-diffusivity, as the external concentration decreases. A mobility decrease that do not even correlate directly with the calcium content of the clay! Note also that TS15 state that the decrease in sodium diffusivity is caused by electromigration (“the diffusion potential term in Equation (57)”). But, as far as I understand, these sodium diffusivities were just put in by hand! (See previous quotation.) It is completely unreasonable that TS15 are here allowed to provide this “explanation” without showing a single result from the actual modeling of how this “diffusion potential term” is supposed to change under different conditions.14

It feels silly to point out that this modeling does not in any way “illustrate the importance of considering coupled diffusion/surface reaction effects in order to understand and to predict migration processes…”. Instead, it is a horrific example of “modeling” as a game, where the only aim is to be able to fit a model to some arbitrary data.

The overall impression of these modeling cases is that TS15 do not seem interested in understanding or communicating how compacted bentonite really works: models are not fully described, and components seem included or disregarded pretty much on a whim; no attention is paid to the quality of the data that is modeled, nor to how data is chosen; and discussions are completely absent regarding how the modeling can aid process understanding. I also find it curious how little these modeling cases connect with the main topics discussed in the rest of the article. For example, very little is said about the impact of multi-component diffusion, even though “non-Fickian” diffusion is a a major theme of the article (it is often not even clear when and how multi-component diffusion has been implemented). Likewise, although sorption processes are included in some of the model cases, we are not told anything about their significance. And not a word is spoken about evolving volumes of the various “porosity domains”, even though such processes seem to be taken seriously in other parts of the article.

“Summary and perspective”

We have finally arrived at the last page of the article. The first paragraph of this section mainly concerns the shortcomings of the traditional diffusion-sorption model (“the classical Fickian diffusion model”)

In this context, the classical Fickian diffusion model applied in the framework of a single pore diffusion model with linear adsorption processes (\(K_D\) model) appears to be inappropriate for describing diffusion data without making less than satisfying modeling assumptions, such as i) different definitions of the porosity as a function of the nature of the tracer of interest (see anion accessible porosity) and as a function of the conditions (anion accessible porosity change with time; ii) change of the adsorption parameters from batch to diffusion experiment; and iii) physically unrealistic pore diffusion coefficients (or tortuosity values).

I agree with that the traditional diffusion-sorption model is not suitable for bentonite. But this has been known for a long time, as I have argued from the beginning of this review. That TS15 use this model as a starting point — and argues for modifying it rather than dismissing it — is in my head a recipe for disaster.

I want to emphasize that all empirical evidence points to that compacted water saturated bentonite is a “single pore” material. The main problem with the traditional diffusion-sorption model is not that it assumes a “single pore”, but that this “single pore” is assumed to contain bulk water. Letting go of the bulk water allows for a satisfying bentonite model using a single pore type. I’ve elaborated on these points here and here.

This first paragraph also says

Diffusion processes through clay materials is the result of a complex interplay of transport and non-linear adsorption processes under the influence of electrostatic fields.

This formulation must be interpreted as summarizing what has been covered in the article as a whole. But I don’t agree with this description, and throughout this review we have exposed flaws and shortcomings in much of the presented material. Of course, ion diffusion is generally influenced by electric fields, since the diffusing entities carry charge. But we have shown that the definition of a “diffusion potential” made by TS15 makes no sense and is based on misunderstanding how electric potentials function. With corrected definitions we have further shown that processes described by TS15 as being due to a “diffusion potential” is actually mainly driven by concentration gradients, i.e. the diffusion processes are “Fickian”. I see it as a major problem that many researchers don’t recognize diffusion in compacted bentonite as simpler than usually described.

This sentence also claims that “non-linear adsorption processes” are involved generally in ion diffusion in bentonite. But this is not the case! The ions that the article mostly focuses on — e.g. sodium, calcium, chloride and iodide — don’t show any signs of adsorption15 at all, let alone “non-linear” adsorption.

The second paragraph of this section is fully devoted to “non-linearity of adsorption processes for strongly adsorbed species”, with particular focus on cesium. I find this peculiar, as non-linear sorption is only brought up in a single small part in the article, and as nothing is said about non-linearity in the the reported “application” that involves cesium (Tachi and Yotsuji, 2014)

The third paragraph begins

Recent developments of selected RT [reactive transport] codes that can handle diffusion processes in diffuse layer and interlayer porosities made it possible to model the diffusion data of neutral, anionic and cationic species within the same conceptual framework.

Here is presumed that partitioning of the pore space into “diffuse layer” and “interlayer” domains makes sense. But we have thoroughly discussed that this whole notion — which is based on the fantasy concept of “stacks” as fundamental structural units of compacted bentonite — makes essentially no sense. And, as we noted earlier, TS15 themselves employ these “porosity domains” for convenience, rather than from being constrained by empirical evidence.

Moreover, it is not the ability to define both “diffuse layer” and “interlayer” that makes it possible to treat any type of species “within the same conceptual framework”16 Rather, a coherent model for compacted bentonite is made possible by realizing that exchangeable ions are mobile and by treating the corresponding ion equilibrium when required. This is (and has been) most adequately made in a homogeneous model.

The third paragraph continues

Also, the contribution of the diffuse layer and/or the interlayer to the overall diffusion of ions makes it possible to explain the origin of the apparent acceleration of cation diffusion as compared to water, which otherwise would require unrealistic tortuosity values for cations in classical Fickian diffusion models.

It is not true generally that cation diffusion is “accelerated”. Speaking of “accelerated” cation fluxes is mostly relevant for steady-state in conventional through-diffusion tests at low background concentration. Writing like this makes me suspect that TS15 have not identified this distinction, especially since they do not recognize the importance of interface equilibrium. As I discussed in the first part of the review, e.g. sodium and chloride show almost identical diffusive behavior in bentonite. The additional effects observed in through-diffusion is completely due to Donnan equilibrium at the interfaces to the external solutions, as has been explained in detail long before the publication of the present article. Again, TS15 contrast with “classical Fickian diffusion models” — but it has been fully clear from the start that the traditional diffusion-sorption model does not apply to compacted bentonite. The way TS15 here use the phrasing “diffuse layer and/or interlayer” again indicates that they don’t take their own concepts seriously. The third paragraph ends

RT modeling has also helped to improve our understanding of anomalous diffusion behavior such as that of up-hill diffusion.

This is a horrible statement. That they refer to the “uphill” effect as “anomalous” tells me that they have neither understood this experiment, nor how bentonite works.17 And if they by “RT modeling” mean the previously presented modeling of the “uphill” effect, nothing can be further from the truth than claiming that it has improved understanding. I also find it significant that they don’t connect the “uphill” test to the previous statements in the same paragraph; the “uphill” test was designed to demonstrate that cation diffusion is not simply “accelerated” (it is, in a sense, conventional tracer through-diffusion that is “anomalous”).

The fourth paragraph reads

Despite the successes of these new RT modeling approaches, it must be stressed that the model and its parameters derived from diffusion experiments are not always unique. Two examples given above highlight the fact that several different conceptual models can provide equally good fits of the data. As such, the modeling effort is typically under constrained, a fact that explains the multitude of conceptual and numerical models available in the literature that describe the ionic transport properties of clay media (Leroy et al. 2006; Appelo and Wersin 2007; Gonçalvès et al. 2007; Birgersson and Karnland 2009; Gimmi and Kosakowski 2011; Tachi et al. 2014).

I strongly disagree with this whole line of reasoning. The multitude of models is not explained by them being “under constrained”; it is explained by most of them being overparameterized (which is poison).

Serious modeling work requires that all included components be evaluated and justified on physical grounds. I see very little of this in the contemporary bentonite scientific literature. We can use the “uphill” test modeling to illustrate this point. The model of this test in TS15 includes, among other things, a mechanism for sodium adsorption and immobilization on basal surfaces (a Stern layer). The test, however, has already been demonstrated to be well described by models that don’t include such a mechanism. This test can therefore not in itself be used to motivate the inclusion of a Stern layer. On the contrary, if there are no other reasonable arguments for why sodium is supposed to adsorb and become immobilized in montmorillonite, it is bad modeling practice to include a Stern layer component.

To make an absurd analogy: if I, as a modeler, don’t have to motivate the included components, I may claim that bentonite contains invisible unicorns that eat anions (or whatever), as long as my model can be fitted to some arbitrary data. This may sound funny, but at least in my head, contemporary bentonite models are quite full of “unicorns”, e.g.

  • They include significant amounts of bulk water
  • They rely on the fantasy concept of stacks
  • Equilibrium between various domains is maintained without a mechanism
  • They postulate regions completely devoid of anions
  • Diffusive fluxes from various domains are added
  • Simple ions, e.g. sodium, are assumed to bond covalently (or whatever) to montmorillonite basal surfaces

In the last paragraphs of this section, TS15 argue for that more research is needed to overcome the claimed problem of models being “under constrained”, using e.g. molecular dynamics simulations and developing new observation techniques. In contrast, I firmly believe that for modeling activities to be meaningful, researchers instead need to spend much more effort scrutinizing and analyzing their models — and the data they are fitting them to. This is much of what I’m trying to do on this blog.

Final comments

My main interest in doing this review has not been in the specific article that we have dissected, but in how it reflects the state of a broader part of the bentonite research field. After finishing, my initial opinion is strengthen: this research field is in a terrible state. I am now fully convinced that reasonable support is nowhere to be found for several of the concepts frequently adopted in contemporary bentonite descriptions. This applies in particular to the ideas of “stacks” as fundamental structural units, the presence of bulk water within the bentonite, and a general adsorption/immobilization mechanism for ions on montmorillonite basal surfaces (Stern layer). To me, it is also clear that this research field lacks knowledge of how salt exclusion and swelling actually work in compacted bentonite, as well as insight into how central these phenomena are for an adequate description.

In today’s bentonite scientific literature, you can simply close your eyes and point at random to find articles that peddle much of the type of nonsense listed above. An illustrative (literally!) example of this legacy is found in Chen et al. (2024)

In the past decades, scholars have carried out a lot of research on the pore structure of compacted bentonite, and it is now generally accepted that there are three different types of pores, including inter-aggregate pores, inter-particle pores and interlayer (or inter-laminar) pores (Lloret and Villar, 2007; Wang et al., 2014), as shown in Fig. 1.

[…]

In existing research, water existing inter-aggregate pores, interparticle pores and interlayer pores are usually referred to as free water, double layer water and interlayer water. In general, anions only diffuse in free water because anions are expelled from the interlayer water and double layer water owing to exclusion effect (Glaus et al., 2007; Wu et al., 2020). Adversely [sic] cation diffusion is supposed to occur in all three types of water due to positive adsorption, corresponding to free water diffusion, surface diffusion and interlayer diffusion, respectively (Melkior et al., 2009; Tachi and Yotsuji, 2014; Yang and Wang, 2019)

“Fig. 1” in Chen et al. (2024) looks very similar to this18

The saddest aspect of this quotation is that it is partly correct: “Scholars have carried out a lot of research” for “decades”, and this view is “generally accepted”.

I cannot say that doing this review has answered my question as to why this is the state of the bentonite research field, but maybe a hint is given by the attitude expressed in TS15 towards modeling and model development. I believe this field will have trouble progressing as long as “modeling” is synonymous with getting arbitrary mathematical descriptions to fit arbitrarily chosen data, padded with empty talk.

Footnotes

[1] As I have commented in the earlier parts: TS15 are fond of using the general terms “clays” and “clay minerals”, while it is clear that the publication mainly focus on systems with substantial ion exchange capacity and swelling properties. Here we will continue to use the term “bentonite” for these systems, and ignore the frequent references in TS15 to more general terms.

[2] Such reinterpretation involves the gradient of the bulk water electrical potential — a potential that TS15 put strictly equal to zero. It also involves the gradient of the Donnan potential — a potential that TS15 mistake for the electric potential in the diffuse layer, when it actually describes the potential difference between bulk and diffuse layer.

[3] Tachi and Yotsuji (2014) report “rock capacity factors” for iodide, which erroneously often are interpreted as anion-accessible porosities. TS15, as we have seen, go “all in” on this concept, while simultaneously admitting that it is fictitious. (?!)

[4] We have noted earlier that it seems very difficult to completely eliminate the effects of interface transfer resistance.

[5] Note that the “correct” value is the arbitrarily selected experimental result 0.42, while the other experimental result of 0.64 is ignored. Also remember that the way this parameter varies with ionic strength is not at all explored by TS15, as they only model results from tests with 0.1 M background concentration.

[6] TS15 report “tortuosities”, rather than actual diffusion coefficients. In a bentonite context, the only way such quantities make sense is as a measure of the diffusion coefficients in terms the corresponding bulk water value, as we discuss here.

[7] I object to that the model parameters are “feasible”. It is in my head a completely ludicrous idea that half the pore volume in pure sodium montmorillonite at 0.8 g/cm3 should contain bulk water — the swelling pressure for a 0.1 M NaCl background concentration is expected to be above 0.4 MPa! Also, the idea that the diffusivity of cesium in montmorillonite is much larger than for sodium, iodide or water is completely unreasonable. In fact, TS15 state that “feasible” simply means a cesium diffusivity that is smaller than for pure bulk water. This shows that they are under the impression that “\(D_e\)”, rather than “\(D_a\)”, quantifies the actual diffusivity in the clay. For a further discussion on this, see here.

[8] Including filter diffusivity in this modeling is certainly not irrelevant per se. However, TS15 merely assign a seemingly random value to this quantity and do nothing more about it. All subsequent inferences are thus contingent on this particular value, yet the authors have neither demonstrated if nor how it influences the overall process.

[9] Although I agree with that these systems don’t contain any appreciable amount of bulk water, we note that no reference is given for this claim. I find this annoying, as TS15 simultaneously don’t hesitate to ascribe 50% bulk water to pure Na-montmorillonite at 0.8 g/cm3. As commented earlier, such a system is expected to have swelling pressure above 0.4 MPa.

[10] We still found serious problems with the model formulation, as further discussed in part 5 of the review.

[11] Also remember that ion exchange is Donnan equilibration!

[12] As always, we should also remember that the distinction betweeen “diffuse layers” and “interlayers” (with quotation marks), or similarly, between “inner” and “outer” basal surfaces, is pure fantasy.

[13] It’s not completely straightforward to extract exact MD values from Tertre et al. (2015). But, discussing them in terms of “surface mobility”, they write: “the mobility for Na is variable from approximately 0.04 to approximately 0.015 when the equivalent fraction of Ca in the solid increases from 0 to 1.” This is a decrease by a factor 2.7.

[14] Note, for instance, that sodium is set to be faster than calcium for the cases of 1.0 M and 0.1 M NaCl, while it is set to be slower than calcium for the cases of 0.05 M and 0.003 M NaCl. Is calcium supposed to be able to retard the sodium flux to such an extent that the sodium ion in the end diffuses slower than the calcium ion?! Note also that the sodium mobility in the TS15 model is lowered by an order of magnitude as the external concentration is only lowered from 0.1 M to 0.05 M! (Or, if the corresponding figure rather than the text is correct, the mobility is even lowered by a factor 40 (the later sections of the paper contain an unacceptable amount of typos of this kind).)

[15] Of course, if ion exchange is supposed to be an “adsorption” process, any cation do “adsorb”. But using such a label for a process that actually may enhance the diffusive flux is, in my opinion, a very bad idea.

[16] To me, the “conceptual framework” promoted in TS15 is “add stuff as you please (but keep bulk water at all costs)”.

[17] Calling the “uphill” effect “anomalous” also seems to contradict the previous general statement that “Diffusion processes through clay materials is the result of a complex interplay of transport and non-linear adsorption processes under the influence of electrostatic fields.”

[18] This only represents a part of “Fig. 1” in Chen et al. (2024), which is explicitly stated to represent compacted water saturated bentonite. I will not elaborate on what is wrong with this figure, but here are some questions to ask: 1. Does this figure represent a realistic size distribution of montmorillonite TOT-layers? 2. Does it represent a typical size distribution of accessory minerals? 3. What is the effective montmorillonite dry density of the depicted system? 4. How is swelling pressure maintained? 5. What keeps “particles” and “aggregates” together? 6. How did this illustration pass peer-review?

Semi-permeability, part III: The membrane potential

Lately on the blog we have had a certain focus on the electric potential in compacted bentonite systems. In the ongoing post-publication review of Tournassat and Steefel (2015), for example, we have deep-dived into how electric potentials are treated in multi-porous models. We have also discussed the seeming “uphill” diffusion effect in different contexts. Although this effect is principally explained by the same mechanisms as for traditional cation tracer through-diffusion, we should address the potential difference that is present across the sample in an “uphill” test.

Here is presented a somewhat more general description of the electric potential difference across a bentonite sample, within the homogeneous mixture model. This treatment is essentially equivalent to how this so-called membrane potential is conventionally evaluated within the ion-exchange membrane scientific discipline. As usual on the blog, the message here is that compacted bentonite should not be modeled with various “porosities” and sorption “sites” (for ion exchange). Rather, from the perspective of ion transport,1 compacted bentonite appears to be well described as a conventional charged ion-exchange membrane.2

Multi-component diffusion in the homogeneous mixture model

Consider a sample of compacted bentonite sandwiched between two external reservoirs.

The external solutions (labeled “1” and “2”) are specified by general species compositions, with the only constraint that they should be electrically neutral. We consequently consider a situation with possible chemical gradients — and corresponding mass transfer — across the bentonite sample, while the hydrostatic pressures in the two reservoirs are assumed equal.3 The values of the electric potential in the external solutions have been labeled \(\psi^{(1)}\) and \(\psi^{(2)}\), respectively.

Treating the bentonite as a single homogeneous phase, the electric potential has discontinuities at the external solution/clay interfaces, due to Donnan equilibrium. In addition, the potential may vary within the bentonite, due to diffusion of charged species with different mobilities (electromigration). Schematically, the steady-state electric potential profile can be illustrated like this

Here, \(\Delta \psi_\mathrm{Donnan}^{(1)}\) and \(\Delta \psi_\mathrm{Donnan}^{(2)}\) denote the potential changes due to Donnan equilibrium at the two interfaces.4 \(\Delta \psi_\mathrm{diff}\), which conventionally is named the diffusion potential, denotes the total change in potential across the sample due to electromigration. The membrane potential, \(\Delta \psi_\mathrm{m}\), can be expressed

\begin{equation} \Delta \psi_\mathrm{m} = \psi^{(2)} – \psi^{(1)} = \Delta \psi^{(1)}_\mathrm{Donnan} + \Delta \psi^{(2)}_\mathrm{Donnan} + \Delta \psi_\mathrm{diff} \tag{1} \end{equation}

Given the reservoir concentrations, the Donnan contributions can be calculated using the general framework, if also a value of \(c_\mathrm{IL}\) — which expresses the montmorillonite structural charge as a monovalent interlayer species concentration — is specified for the bentonite component. Such calculations also provide the corresponding sets of internal interface concentrations, i.e. \(\left \{c_{1}^\mathrm{int}(0), c_{3}^\mathrm{int}(0), c_{3}^\mathrm{int}(0), \ldots \right \}\), and \(\left \{c_{1}^\mathrm{int}(L), c_{3}^\mathrm{int}(L), c_{3}^\mathrm{int}(L), \ldots \right \}\), where we denote by \(c_i^\mathrm{int}(x)\) the concentration for species \(i\) at position \(x\) within the clay, and we assume the clay domain located between \(x=0\) and \(x=L\).

The sets of internal interface concentrations, in turn, constitute boundary conditions for solving e.g. the Nernst-Planck equation in the clay domain. Within the Nernst-Planck framework, the gradient of the electrostatic potential in the clay is expressed in terms of species concentrations and diffusivities as

\begin{equation} \nabla \psi(x) = – V_T \cdot \frac{\sum z_i\cdot D_i \cdot \nabla c_i^\mathrm{int} (x) }{\sum z_i^2\cdot D_i \cdot c_i^\mathrm{int} (x)} \tag{2} \end{equation}

where \(z_i\) and \(D_i\) denote, respectively, the charge number and diffusion coefficient in the clay domain, for species \(i\), and \(V_T = RT/F\) is the thermal voltage (we assume \(V_T\) = 25.7 mV here). From eq. 2 the diffusion potential can be evaluated.

Note that the problem of ending up with a too constrained model when both the Donnan equilibrium and the Nernst-Planck frameworks are invoked does not occur here. This is in contrast to the problems we recently discussed in the model of Tournassat and Steefel (2015). That model requires Donnan equilibrium at each point within the clay, while here it is only imposed at the two interface points.

In the following we will evaluate the membrane potential numerically in some different cases. We start, however, by analyzing the 1:1 system analytically, in order to get a better feel for these processes in compacted bentonite.

1:1 system

The pure 1:1 system has a single type of monovalent cation and a single type of monovalent anion (e.g. Na-montmorillonite in contact with NaCl solutions). The requirement of charge neutrality dictates equal concentrations of cations (\(c_+\)) and anions (\(c_-\)) in each external solution, and we write

\begin{equation} c_+^{(n)} = c_-^{(n)} = c^{(n)} \tag{3} \end{equation}

where \(c^{(n)}\) is the electrolyte concentration in reservoir \(n\) (\(n\) = 1 or 2).

In this analytical treatment, we assume the often relevant condition of having small external concentrations as compared with the concentration of structural charge in the clay, i.e.

\begin{equation} c^{(n)} \ll c_\mathrm{IL},\tag{4} \end{equation}

where the amount of structural charge is quantified by

\begin{equation} c_\mathrm{IL} = \frac {CEC\cdot\rho_w}{F\cdot w}. \tag{5} \end{equation}

Here \(CEC\) is the cation exchange capacity, \(w\) the water-to-solid mass ratio, \(\rho_w\) is water density and \(F\) the Faraday constant.

Under the condition expressed in eq. 4, the so-called Donnan factors (\(f_D\)) can be approximated as

\begin{equation} f_D^{(n)} \equiv e^{\psi^{\star, (n)}/V_T}\approx \tag{6} \frac{c^{(n)}}{c_\mathrm{IL}} \end{equation}

where \(\psi^{\star, (n)}\) is the Donnan potential. In eq. 6 we have ignored the activity coefficient correction factor; contributions to the membrane potentials from activity coefficient differences are expected to be small, and here we will not consider them.5 Note that Donnan potentials are negative, and that \(\Delta \psi_\mathrm{Donnan}^{(1)} = \psi^{\star, (1)}\) and \(\Delta \psi_\mathrm{Donnan}^{(2)} = -\psi^{\star, (2)}\). Utilizing eq. 6, we can express the total Donnan contribution to the membrane potential as6

\begin{equation} \Delta \psi_\mathrm{Donnan} = \Delta \psi^{(1)}_\mathrm{Donnan} + \Delta \psi^{(2)}_\mathrm{Donnan} \approx V_T\ln{\frac{c^{(1)}}{c^{(2)}}} \tag{7} \end{equation}

For the clay domain, the requirement of charge neutrality can be written

\begin{equation} c_+^\mathrm{int}(x) = c_\mathrm{IL}+c_-^\mathrm{int}(x) \tag{8} \end{equation}

from which it follows that the gradients of anion and cation concentrations are equal

\begin{equation} \nabla c_+^\mathrm{int}(x) = \nabla c_-^\mathrm{int}(x) \tag{9} \end{equation}

We can therefore approximate the gradient of the electric potential in the clay (eq. 2) as

\begin{equation} \nabla \psi(x) \approx V_T \frac{\left ( D_-/D_+-1 \right)}{c_\mathrm{IL} } \nabla c_-^\mathrm{int}(x) \tag{10} \end{equation}

where we also have continued assuming low external concentrations (eq. 4).

If we impose a constant internal concentration gradient, i.e. \(\nabla c_-^\mathrm{int} (x) = \Delta c_-^\mathrm{int}/L\), where \(\Delta c_-^\mathrm{int} \equiv c_-^\mathrm{int}(L) – c_-^\mathrm{int}(0)\), the diffusion potential can be written

\begin{equation} \Delta \psi_\mathrm{diff} \approx V_T \frac{\left ( D_-/D_+-1 \right )}{c_\mathrm{IL} } \Delta c^\mathrm{int}_-\tag{11} \end{equation}

Combining eqs. 7 and 11 gives an approximation for the membrane potential for a 1:1 system, valid for low external concentrations

\begin{equation} \psi_\mathrm{m} \approx \left ( \ln{\frac{c^{(1)}}{c^{(2)}}} + \frac {\left ( D_-/D_+-1 \right )}{c_\mathrm{IL} } \Delta c^\mathrm{int}_- \right ) V_T \tag{12} \end{equation}

The first term — the ideal Donnan contribution — depends only on the ratio of the concentration in the two reservoirs, and usually dominates over the diffusion potential (second term).

The diffusion potential is expected to be small under most circumstances, as it is seen to explicitly depend on \(1/c_\mathrm{IL}\). This behavior reflects that potential variations within the bentonite are reduced due to the presence of counter-ions (bentonite is a decent conductor). But the factor \(\Delta c_-^\mathrm{int}\) in eq. 11 contributes with an additional factor \(1/c_\mathrm{IL}\) to the diffusion potential, because of Donnan equilibrium (see eq. 6)

\begin{equation} \Delta c_-^\mathrm{int} = f_D^{(2)}c^{(2)} – f_D^{(1)}c^{(1)} \approx \frac {\left ( c^{(2)} \right )^2 – \left ( c^{(1)} \right )^2} {c_\mathrm{IL}} \tag{13} \end{equation}

Note also that the diffusion potential only depends on the relative values of the diffusivities, and consequently vanishes if there is no difference in mobility between the cation and the anion.

Comparisons with experiment

We can compare the present theory and approximations with measurements reported in Yaroshchuck et al. (2007). In this study, the membrane potential (called the “concentration potential”) was measured across sodium montmorillonite of density 2.0 g/cm3, with one external solution kept at 0.1 M NaCl and the other raised to 0.2 M, 0.3 M, or 0.5 M. The very high density corresponds to a value of \(c_\mathrm{IL} \approx\) 6.0 M,7 and thus eq. 12 applies (see eq. 4).

In order to also investigate the influence of a diffusion potential we should adopt values for the individual diffusivities of chloride and sodium. Experimental tracer diffusion results indicate that these ions have quite similar mobility. E.g., in pure Na-montmorillonite at 1.8 g/cm3 and 298 K, Kozaki et al. (1999, 2001, 2005) report values in the range \(1.5\cdot 10^{-11}\) — \(2.5\cdot 10^{-11}\) m2/s for chloride, and \(2.1\cdot 10^{-11}\) — \(2.4\cdot 10^{-11}\) m2/s for sodium. From these values we can estimate a range for the factor (\(D_-/D_+-1\)) of -0.4 — 0.2. Rather than adopting such a value, however, we here assume the diffusivities to have the corresponding bulk water ratio (1.53), giving \((D_-/D_+ -1)\) = 0.53. This is similar to the treatment in the examples in Tournassat and Steefel (2015), which we have discussed recently. Theoretical and experimental results are compared in the table below

Evaluations using cIL = 6.0 M, c(2) = 0.1 M, (D/D+ – 1) = 0.53:

This table compares calculations performed using eq. 12 (columns 2 — 6), calculations performed using the full frameworks for Donnan equilibrium and the Nernst-Planck potential (columns 7 — 11), and experimental values (last column). We note that, for this very high value of \(c_\mathrm{IL}\), the calculated Donnan contributions (\(\Delta \psi_\mathrm{D}\)) are essentially identical whether or not the full framework is adopted. We also see that the calculated diffusion potential contributions are completely negligible: the largest value calculated is still below 0.1 mV! Finally, we note excellent agreement between calculated and measured membrane potentials.

The agreement with experiment suggests, not only that the membrane potential is completely determined by Donnan equilibrium in these systems, but that bentonite essentially play by the same rules as conventional ion-exchange membranes.2 In particular, these results8 are even more evidence for that all exchangeable cations are mobile in bentonite.

To explore the NaCl system further, let’s investigate how the results change for lower clay density. In the tables below are shown corresponding results as above, but for \(c_\mathrm{IL}\) = 3.0 M and \(c_\mathrm{IL}\) = 1.5 M, which roughly correspond to densities 1.6 g/cm3 and 1.1 g/cm3, respectively.

Evaluations using cIL = 3.0 M, c(2) = 0.1 M, (D/D+ – 1) = 0.53:

Evaluations using cIL = 1.5 M, c(2) = 0.1 M, (D/D+ – 1) = 0.53:

For \(c_\mathrm{IL}\) = 3.0 M, the difference between using the full framework or eq. 12 is seen to still only affect the evaluated membrane potentials on the order of tenths of mV. Only for the higher external concentrations in the system with \(c_\mathrm{IL}\) = 1.5 M are the differences becoming above 1 mV. Similarly, the diffusion potential is seen to be essentially negligible for all these cases.

Further investigation of potentials and fluxes

It should be emphasized that we are not actually solving the full Nernst-Planck transport framework in these examples. Rather, we calculate the corresponding potential difference as we impose a linear internal concentration profile. This is similar to the evaluation made by Tournassat and Steefel (2015) — although they impose concentration profiles in a presumed bulk water phase within the clay. Similar to that publication, we will here visualize the corresponding fluxes in some cases of interest.

The total flux of a charged species in the Nernst-Planck framework can be written as the sum of a “Fickian” contribution (i.e. proportional to the corresponding concentration gradient) and a contribution from electromigration (i.e. flux induced by the presence of an electric field).9

\begin{equation} j_\mathrm{tot} = j_\mathrm{conc} + j_\mathrm{E} \tag{14} \end{equation}

where the “Fickian” flux is given by

\begin{equation} j_\mathrm{conc} = -D\cdot\nabla c^\mathrm{int}(x) \tag{15} \end{equation}

and the electromigration term is given by

\begin{equation} j_\mathrm{E} = -z\cdot D\cdot c^\mathrm{int}(x)\frac{\nabla \psi(x)}{V_T} \tag{16} \end{equation}

NaCl

Below are plotted imposed concentration profiles, Nernst-Planck fluxes, and the electric potential profile, for the NaCl-system with \(c_\mathrm{IL}\)=3.0 M, \(c^{(1)}\) = 0.5 M, and \(c^{(2)}\) = 0.1 M (for which we previously calculated \(\Delta \psi_\mathrm{m}\) = 40.4 mV and \(\Delta \psi_\mathrm{diff}\) = -0.3 mV). To calculate fluxes we must choose absolute values of the diffusion coefficients, which we take to be \(D_\mathrm{Na}=1.33\cdot 10^{-10}\) m2/s and \(D_\mathrm{Cl}=2.03\cdot 10^{-10}\) m2/s.10

The total NaCl flux is essentially constant throughout the sample (we have adopted a sample length of 10 mm), demonstrating that the imposed linear concentration profile in essence is the steady-state solution. We see that the flux of chloride — which has a much smaller concentration than sodium — is completely dominated by the concentration gradient contribution (eq. 15), while the sodium flux has significant contributions from both the concentration gradient and the electric field (eq. 16). The picture is clear: the faster chloride ions diffuse with a rate set by their individual diffusion coefficient, while the induced electric field boosts the less mobile (but more abundant) sodium ions in order to “keep up”. The behavior is similar to what would be experienced in a bulk solution with a NaCl gradient. In that case, however, the chloride is as much retarded as the sodium is boosted, resulting in a diffusivity for the full salt that is between the two individual values. In compacted bentonite, on the other hand, the electromigration flux is mainly significant for sodium (the counter-ions). Therefore the full salt diffuse with a rate essentially set by the chloride mobility.

We can quantify the relative influcence of electromigration for a specific species by taking the ratio of eqs. 15 and 16, and also utilizing the approximation for \(\nabla \psi (x)\) in eq. 10

\begin{equation} \frac{j_\mathrm{E}}{j_\mathrm{conc}} \approx \frac{z\cdot c^\mathrm{int}(x)} {c_\mathrm{IL}} \left (D_\mathrm{Cl}/D_\mathrm{Na}- 1 \right ) \tag{17} \end{equation}

This expression explicitly shows that the electromigration flux is negligible for chloride (because \(c^\mathrm{int}_\mathrm{Cl}(x) \ll c_\mathrm{IL}\)), while the \(j_\mathrm{E}/j_\mathrm{conc}\)-ratio for sodium is essentially determined by the factor \(\left (D_\mathrm{Cl}/D_\mathrm{Na} – 1 \right )\) (because \(c^\mathrm{int}_\mathrm{Na}(x) \approx c_\mathrm{IL}\)).

For a similar NaCl-system at lower density, the picture is quite similar (here \(c_\mathrm{IL}\)=1.5 M, with the other parameters as previously; \(\Delta \psi_\mathrm{m}\) = 37.8 mV and \(\Delta \psi_\mathrm{diff}\) = -1.1 mV)

Let’s also play with the relative mobility between anion and cation. In the plot below we have decreased the chloride diffusivity to \(D_\mathrm{Cl}=0.133\cdot 10^{-10}\) m2/s, while keeping with the same values as before for the other parameters (for \(c_\mathrm{IL}\) = 3.0 M).

In this case, the diffusivity factor for the electromigration is \((D_-/D_+ -1) = -0.9\). The negative value indicates an electric field directed in the opposite direction as compared to the previous cases (chloride is now the slower ion). The diffusion potential is \(\Delta \psi_\mathrm{diff}\) = +0.6 mV, which still is quite small, given that the chloride is assumed to have ten times lower mobility as compared with sodium. The evaluated membrane potential is \(\Delta \psi_\mathrm{m}\) = 41.3 mV. We note that the flux is constant throughout the domain and that chloride is completely governed by the concentration gradient contribution, while the sodium flux is composed of two oppositely directed much larger contributions of similar magnitude.

If we instead make the sodium diffusivity 10 times smaller than the chloride diffusivity (\(D_\mathrm{Na}=0.203\cdot 10^{-10}\) m2/s, \(D_\mathrm{Cl}=2.033\cdot 10^{-10}\) m2/s), the effect is more dramatic

With \((D_-/D_+ -1)\) = 9.0, the diffusion potential becomes considerably larger than in any of our previously investigated cases: \(\Delta \psi_\mathrm{diff}\) = -5.2 mV (\(\Delta \psi_\mathrm{m}\) = 35.5 mV). This is reasonable, since a significantly smaller mobility of the counter-ions requires a larger electric field to produce similar fluxes. The sodium flux is therefore dominated by electromigration, and we also note a non-negligible contribution from electromigration for chloride. The overall flux varies quite significantly throughout the clay domain; we’re not really in steady-state.

“Uphill” diffusion

With a solid understanding of the principles for diffusion of a 1:1 electrolyte in compacted bentonite, let’s turn our attention to the “uphill” test (Glaus et al., 2013). The conditions for the main electrolyte (NaClO4) and clay in this study are the same as in our previous case: one reservoir is at concentration 0.1 M, and one at 0.5 M, and the clay is pure high-density Na-montmorillonite. For a further discussion of this study, see this previous post. Here we adopt the value \(c_\mathrm{IL}\) = 4.5 M and set \(c^{(1)}\) = 0.1 M and \(c^{(2)}\) = 0.5 M (in order to have a positive tracer flux). In addition to a main electrolyte, the “uphill” test has a sodium tracer (\(^{22}\mathrm{Na}\)) as a third component, with equal concentration in the two reservoirs (\(10^{-11}\) M). Here we choose sodium diffusivity \(D_\mathrm{Na}=3\cdot 10^{-11}\) m2/s, similar to what is evaluated in Glaus et al. (2007) and Birgersson and Karnland (2009). For perchlorate we choose \(D_\mathrm{ClO_4}=4\cdot 10^{-11}\) m2/s, in order to keep the same diffusivity ratio as for bulk (the perchlorate bulk value is \(1.79\cdot 10^{-9}\) m2/s (Heil, 1995)).

Since the only difference between the “uphill” study and the pure 1:1 case is the presence of a tracer, it is no surprise that concentrations, fluxes and potential behave very similarly

The membrane potential is here \(\Delta \psi_\mathrm{m} = -41.0\) mV, which is essentially identical to what we have evaluated for the NaCl system (the minus sign simply reflects that the left reservoir here has the lower concentration). With a smaller mobility difference between anion and cation, the diffusion potential is even more negligible than for NaCl, \(\Delta \psi_\mathrm{diff} = 0.1\) mV. The tracer flux is not resolved in the above plot, but looks like this

An analytical evaluation of the \(j_\mathrm{E}/J_\mathrm{conc}\)-ratio becomes a bit more complicated for the tracer as compared with the main salt (eq. 17). At the two interfaces, it can be shown that11

\begin{equation} \frac{j_\mathrm{E,^{22}Na}}{j_\mathrm{conc,^{22}Na}} \approx -\frac{c^\mathrm{(2)} \left ( c^{(1)} + c^{(2)} \right )} {c_\mathrm{IL}^2} \left (D_\mathrm{ClO_4}/D_ \mathrm{Na}- 1 \right ),\;\;\;\; x=0 \tag{18} \end{equation}

and

\begin{equation} \frac{j_\mathrm{E,^{22}Na}}{j_\mathrm{conc,^{22}Na}} \approx -\frac{c^\mathrm{(1)} \left ( c^{(1)} + c^{(2)} \right )} {c_\mathrm{IL}^2} \left (D_\mathrm{ClO_4}/D_ \mathrm{Na}- 1 \right ),\;\;\;\; x=L \tag{19} \end{equation}

We see that this ratio is expected to be small as long as the reservoir concentrations are small as compared with \(c_\mathrm{IL}\). This is also confirmed in the above plot; the electromigration contribution is completely negligible.

This result is the ultimate support for the claim I made in the “uphill” blog post: the phenomena is first and foremost an effect of ordinary Donnan equilibrium, just as conventional cation tracer diffusion is explained. Here we have explicitly included the electric potential in the analysis, and not only showed that electromigration is negligible for the tracer component, but qualitatively described how and where possible effects due to ion mobility differences occur in the set-up as a whole.

Tertre et al. (2024) — the publication critisized in the “uphill” blog post — states that the “uphill” study

[…] demonstrated the marked influence of background electrolyte concentration gradient on tracer diffusion, and thus the necessity to understand the couplings between diffusion of several charged species present at contrasting concentrations and experiencing different concentration gradients.

In contrast to Tertre et al. (2024), we have here presented this understanding and — in contrast to the statement — demonstrated that the “uphill” effect does not depend critically on any additional mechanisms or couplings as compared with conventional cation tracer diffusion.

As for the previous cases, we can play around with the parameters to explore more “extreme” conditions in an “uphill” test. Here we have increased the perchlorate diffusivity to \(D_\mathrm{ClO_4}=9\cdot 10^{-11}\) m2/s (i.e. three times the sodium diffusivity), and lowered \(c_\mathrm{IL}\) to \(c_\mathrm{IL} = \) 2.0 M (corresponding roughly to a density of 1.25 g/cm3).

Predictably, the main counter-ion transport is now dominated by electromigration, and we see some influence of electromigration also on the perchlorate and the sodium tracer. Nevertheless, even for this extreme case, the “uphill” tracer flux is mainly “Fickian” and is principally explaind as an effect due to Donnan equilibrium at the sample interfaces. For this case we have \(\Delta \psi_\mathrm{m}\) = -37.4 mV and \(\Delta \psi_\mathrm{diff}\) = 2.6 mV.

CaCl2

Let’s also examine the pure 2:1 system, which corresponds e.g. to CaCl2 solutions in contact with Ca-montmorillonite. As we discussed in a recent blog post, the Donnan equilibrium for this system is quite different as compared with the 1:1 system. This will naturally also impact the membrane potential. Moreover, experiments indicate that calcium diffuses several times slower than chloride in bentonite (Kozaki et al., 2000), which suggests a larger effect on the diffusion potential. But, similar to previous cases, let’s begin by assuming the same diffusivity ratio as for bulk. With \(D_\mathrm{Cl}=2.03\cdot 10^{-10}\) m2/s, \(D_\mathrm{Ca}=0.79\cdot 10^{-10}\) m2/s, \(c_\mathrm{IL} = \) 3.0 M, and having \(c^{(1)}\) = 0.25 M and \(c^{(2)}\) = 0.05 M12 we get the following picture

The membrane potential is here \(\Delta \psi_\mathrm{m} = 18.8\) mV and the diffusion potential \(\Delta \psi_\mathrm{diff} = -1.1\) mV. Thus, even if the relative influence of the diffusion potential is larger han for NaCl, the main contribution to the membrane potential is still due do Donnan equilibrium; the considerably lower membrane potential in this case (about half of the NaCl value with the same amount of chloride in the reservoirs) is primarily a consequence of Donnan equilibrium: the Donnan factors have in this case a square-root dependence on external concentration. However, as a consequence of calcium mobility being considerably lower than sodium mobility (\(D_-/D_+ – 1\) = 1.57), the magnitude of the diffusion potential is about three times larger than for the NaCl case. The resulting flux is also seen to be not completely constant throughout the clay domain. Nevertheless, the same qualitative picture of the transport mechanism as for NaCl is valid here: chloride transport is dominated by its “Fickian” contribution, diffusing with the rate set by the individual mobility, while the less mobile calcium is boosted by an electromigration contribution, which here is actually larger than its concentration gradient contribution. The salt as a whole diffuses with the rate essentially set by the chloride diffusivity.

Let’s also explore an “extreme” case. Below we have set the Ca mobility to be a factor 10 lower than Cl (\(D_\mathrm{Cl}=2.03\cdot 10^{-10}\) m2/s, \(D_\mathrm{Ca}=0.203\cdot 10^{-10}\) m2/s) and lowered \(c_\mathrm{IL}\) to \(c_\mathrm{IL} = \) 1.5 M.

Here the membrane potential is only \(\psi_\mathrm{m} = 6.7\) mV due to a large contribution from the diffusion potential (\(\Delta \psi_\mathrm{diff} = -12.0\) mV), The calcium flux is seen to be strongly dominated by electromigration, and chloride transport is also strongly affected. We are far from steady-state.

Mixing sodium and calcium

As a final case we consider systems containing both a mono- and a di-valent type of cation. Since cation concentrations are enhanced within the clay, a significant diffusion potential may develop, even for relatively small variations of the external solutions. For instance, with a pure 0.1 M NaCl solution on one side (2), and a mixture of 0.003 M CaCl2 and 0.094 M NaCl on the other (1), we get the following result

Here we have assumed \(c_\mathrm{IL} = \) 3.0 M, \(D_\mathrm{Cl}=2.03\cdot 10^{-10}\) m2/s, \(D_\mathrm{Ca}=0.79\cdot 10^{-10}\) m2/s, and \(D_\mathrm{Na}=1.33\cdot 10^{-10}\) m2/s.

Note that the Donnan equilibrium at the left interface acts as to make the clay contain 50% sodium and 50% calcium (charge wise) even though the calcium content in the reservoir is only 6%! A quite minor “disturbance” in the external solution consequently induce huge transport processes.13 This transport is essentially independent of the presence of an anion, and is governed by mutual diffusion of calcium and sodium. The mobility difference for these two ions thus results in a large diffusion potential: \(\Delta \psi_\mathrm{diff}\) = -5.0 mV. The membrane potential is quite substantial: \(\Delta \psi_\mathrm{m}\) = 11.4 mV. The calcium transport is dominated by the concentration gradient flux, while the oppositely directed sodium flux is more profoundly retarded by the electric field (although it is still dominated by the concentration gradient). We conclude that the full Na-Ca transport process occurs approximately at the rate set by the individual mobility for calcium. This process is essentially charge neutral, i.e. as much Ca charge is transported to the right as Na charge is transported to the left.

It may also be worth noticing that although the chloride concentration is the same in the two external solutions, an internal gradient — and a corresponding chloride flux — is induced

This chloride flux, which is directed to the right and seen to be quite influenced (retarded) by electromigration, is more than two orders of magnitude smaller than the Ca-Na flux. This result implies that also anions may experience seeming “uphill” diffusion (though I’m not sure how easy it would be to detect experimentally).

A lesson here is that calcium concentrations often must be very small in order to be treated as a tracer when added to a sodium system. For example, Tinnacher et al. (2016) incorrectly treats calcium as a tracer when adding 0.001 M CaBr2 to an otherwise pure 0.1 M NaCl-montmorillonite system at 0.8 g/cm3.14 At that density, we may approximate \(c_\mathrm{IL} = 1.0\) M, which gives the following situation (\(D_\mathrm{Cl} = 2.03 \cdot 10^{-10}\) m2/s, \(D_\mathrm{Br}=2.01\cdot 10^{-10}\) m2/s, \(D_\mathrm{Ca}=0.79\cdot 10^{-10}\) m2/s, and \(D_\mathrm{Na}=1.33\cdot 10^{-10}\) m2/s; CaBr2 in reservoir 1)

The amount of calcium in the clay at the left interface is in this case predicted to be about 15% — which certainly cannot be regarded a trace level. Consequently, the main mode of transport in this test is oppositely directed sodium and calcium fluxes that essentially does not involve the anions. In addition, we predict a non-zero membrane potential of \(\Delta \psi_\mathrm{m} = 2.5\) mV, with \(\Delta \psi_\mathrm{diff} = -1.5\) mV.

The bromide, in contrast, can be treated as tracer, and its flux look like this

This flux is about two orders of magnitude smaller than the dominating Ca-Na flux,15 and is seen to be essentially unaffected by the electric field (i.e. mainly “Fickian”).

Additional comments

Here we have investigated transport processes that occur in the presence of chemical gradients between two external reservoirs sandwiching a compacted bentonite component. For this case, compacted bentonite appears to function similar to a conventional charged ion exchange membrane.2 However, compacted bentonite used e.g. for waste storage is quite different from ordinary ion exchange membranes, in that we are more interested in what is going on inside it, rather than directly utilizing its semi-permeable properties. The results discussed in this blog post are useful when interpreting tests that involve bentonite sandwiched between external solutions. For a description of what is going on inside the bentonite component, such results and tests may be very useful, but do not provide a full picture.

As with many topics discussed on the blog, I find the lack of relevant experimental results quite extraordinary. I mean that many labs should routinely measure membrane potentials as part of exploring bentonite systems. Instead, as far as I understand, the study by Yaroschchuck et al. (2007) is the only one published on actual compacted bentonite.16 Moreover, this study appears to be completely forgotten!8

Yaroschchuck et al. (2007) use extremely high density, and it would be interesting with data also at lower densities (1.6 — 1.3 g/cm3, say); the present theory predicts only minor changes in membrane potentials. A “simple” further test of the present description would be to measure the membrane potential in e.g. a pure CaCl2-montmorillonite system. We end this post by presenting such potentials as predicted from the present theory (the calculation uses \(D_\mathrm{Cl} = 2\cdot 10^{-10}\) m2/s and \(D_\mathrm{Ca} = 0.5\cdot 10^{-10}\) m2/s).

Footnotes

[1] The focus here is on what I call simple ions, i.e. ions that only interact with the bentonite by beeing part of a diffuse layer. In bentonite, more complex processes are active for certain species that interact with e.g. montmorillonite edges. Here is presented an extension of the homogeneous mixture model for “truly” sorbing species.

[2] What is presented in this blog post is, as far as I understand, the “Teorell- Meyer-Sievers” model for the membrane potential, applied to compacted bentonite. This model, which was originally developed in the 1930s, is used routinely within the ion-exchange membrane scientific discipline. I’m currently reading up on this literature, and is here using my “usual” terminology (for bentonite). I have not yet come across sources with specific focus on bentonite/montmorillonite, and I’m not sure how well established these are as a membrane materials within this discipline. In a well-cited review by Xu (2005), bentonite is mentioned once — misspelled! Perhaps this is indicative of a missing link between the ion exchange membrane and compacted bentonite research fields.

[3] Although different water chemistries will induce different swelling pressure repsonses, and thus couple to mechanical processes in the bentonite, we here ignore pressure altogether (but we certainly will have reason to return to this).

[4] Note that Donnan potentials generally are negative, while \(\Delta \psi_\mathrm{Donnan}^{(1)}\) and \(\Delta \psi_\mathrm{Donnan}^{(2)}\) are negative or positive, depending on whether the contribution decreases or increases the potential as we move to the right.

[5] Including activity coefficients, eq. 6 becomes

\begin{equation} f_D^{(n)} \approx \gamma^{\mathrm{ext},(n)}/\gamma^\mathrm{int}\cdot c^\mathrm{(n)}/c_\mathrm{IL} \end{equation}

The expression for the membrane potential will contain the ratio of the Donnan factors at each interface

\begin{equation} V_T\ln(f_D^{(1)}/f_D^{(2)}), \end{equation}

and the interlayer activity coefficients (\( \gamma^\mathrm{int}\)) will therefore cancel (to the extent that the this quantity can be considered independent of external solution concentration). The contribution from Donnan equilibrium to the membrane potential thus become

\begin{equation} \Delta \psi_\mathrm{Donnan} \approx V_T(\ln{\frac {c^{(1)}}{c^{(2)}}} + \ln{\frac {\gamma^{\mathrm{ext},(1)}}{\gamma^{\mathrm{ext},(2)}}}). \end{equation}

The primary approximation made in the present treatment is to ignore the second term.

[6] This equation is often referred to as the Nernst equation. Note, however, that this quantity fundamentally involves two semi-permeable components. I have previously referred to a similar relation for a single interface as the Nernst equation.

[7] The value of \(c_\mathrm{IL}\) (eq. 5) depends on the cation exchange capacity of the material: montmorillonite “from Milos”. Different sources report this value in the range 0.80 — 0.85 eq./kg.

[8] The obscurity of Yaroshchuck et al. (2007) amazes me. This paper is part of a “trilogy” published in 2007, by the same authors (to a large extent). The other two studies are featured massively on the blog: Glaus et al. (2007) demonstrates that interlayer diffusion dominates mass transfer in compacted bentonite (see here and here), and Van Loon et al. (2007) is the only study that passed our assessment of published chloride equilibrium concentrations in compacted bentonite (see here, here, and here). Glaus et al. (2007) has at the moment 178 citations on Google Scholar, and Van Loon et al. (2007) has 333 citations. Yaroshchuck et al. (2007), in contrast, has 7 citations… Perhaps even more weirdly, Yaroshchuck et al. (2007) is not cited in Van Loon et al. (2007), although the former study presents the framework for Donnan equilibrium between compacted bentonite and a 1:1 electrolyte, and relates it to the findings of Glaus et al. (2007). It should thus had been an “easy win” for Van Loon et al. (2007) to provide a proper explanation for anion exclusion in bentonite, and to make the correct connection between anion and cation tracer through-diffusion. Instead, and for reasons I can’t get my head around, van Loon et al. (2007) get lost in ideas on “anion-accessible porosity”.

[9] Generally, the flux may also have an advective component, which we don’t consider here.

[10] The absolute values of the diffusivities are chosen rather arbitrarily; in this example we have adopted values ten times smaller than the corresponding values in bulk. The absolute value of the corresponding flux is therefore also arbitrary. Here, however, we are mainly interested in how constant the flux is throughout the clay, and in the relative contributions from concentrations and electric fields, rather than the absolute values. In some cases in this section we have adopted empirical values for diffusivities in order to compare also the absolute value of the flux.

[11] The electromigration and concentration gradient fluxes for the tracer in the “uphill” test can be written (see eqs. 15 and 16)

\begin{equation} j_\mathrm{E,tr} \approx -D_\mathrm{Na} \cdot c_\mathrm{tr}^\mathrm{int}(x) \cdot \left (D_\mathrm{ClO_4}/D_ \mathrm{Na}- 1 \right) \frac{\nabla c^\mathrm{int}_\mathrm{ClO_4}}{c_\mathrm{IL}} \end{equation}

\begin{equation} j_\mathrm{conc,tr} = -D_\mathrm{Na} \cdot \nabla c_\mathrm{tr}^\mathrm{int}(x) \end{equation}

Their ratio is

\begin{equation} \frac{j_\mathrm{E,tr}}{j_\mathrm{conc,tr}} \approx \frac{c_\mathrm{tr}^\mathrm{int}(x)}{c_\mathrm{IL}} \cdot \left (D_\mathrm{ClO_4}/D_ \mathrm{Na}- 1 \right) \cdot \frac{c^\mathrm{int}_\mathrm{ClO_4}(L) – c^\mathrm{int}_\mathrm{ClO_4} (0)}{c_\mathrm{tr}^\mathrm{int}(L) – c_\mathrm{tr}^\mathrm{int}(0)} \end{equation}

where we have explicitly written the gradients in terms of the internal interface concentrations. Using the relation between Donnan factors and concentrations, this concentration difference ratio can in turn be rewritten

\begin{equation} \frac{c^\mathrm{int}_\mathrm{ClO_4}(L) – c^\mathrm{int}_\mathrm{ClO_4} (0)}{c_\mathrm{tr}^\mathrm{int}(L) – c_\mathrm{tr}^\mathrm{int}(0)} = \frac{f_D^{(2)}c^{(2)} – f_D^{(1)} c^{(1)}} {c_\mathrm{tr}/f_D^{(2)} – c_\mathrm{tr}/f_D^{(1)}} = \end{equation} \begin{equation} \frac{1}{c_\mathrm{tr}c_\mathrm{IL}^2} \frac{(c^{(2)})^2 – (c^{(1)})^2} {1/c^{(2)} – 1/c^{(1)}} = -\frac{c^{(1)}c^{(2)}}{c_\mathrm{tr}c_\mathrm{IL}^2} \left (c^{(1)} + c^{(2)} \right ) \end{equation}

where we have written \(c_\mathrm{tr}\) for the tracer concentration in the reservoirs (it’s equal in the two solutions). The flux ratio can thus be written

\begin{equation} \frac{j_\mathrm{E,tr}}{j_\mathrm{conc,tr}} \approx -\frac{c_\mathrm{tr}^\mathrm{int}(x)c^{(1)}c^{(2)}}{c_\mathrm{tr}\cdot c_\mathrm{IL}^3} \cdot \left ( c^{(1)} + c^{(2)} \right ) \cdot \left (D_\mathrm{ClO_4}/D_ \mathrm{Na}- 1 \right) \end{equation}

At interface (1), \(c_\mathrm{tr}^\mathrm{int}(0) =c_\mathrm{tr}/f_D^{(1)} = c_\mathrm{tr}\cdot c_\mathrm{IL}/c^{(1)}\), giving

\begin{equation} \frac{j_\mathrm{E,tr}}{j_\mathrm{conc,tr}} \approx -\frac{c^{(2)}\cdot \left (c^{(1)} + c^{(2)} \right )} {c_\mathrm{IL}^2} \cdot \left (D_\mathrm{ClO_4}/D_ \mathrm{Na}- 1 \right), \;\;\;\; x=0 \end{equation}

At interface (2), \(c_\mathrm{tr}^\mathrm{int}(L) =c_\mathrm{tr}/f_D^{(2)} = c_\mathrm{tr}\cdot c_\mathrm{IL}/c^{(2)}\), giving

\begin{equation} \frac{j_\mathrm{E,tr}}{j_\mathrm{conc,tr}} \approx -\frac{c^{(1)}\cdot \left (c^{(1)} + c^{(2)} \right )} {c_\mathrm{IL}^2} \cdot \left (D_\mathrm{ClO_4}/D_ \mathrm{Na}- 1 \right), \;\;\;\; x=L \end{equation}

[12] These reservoir concentrations correspond to having the same chloride concentration (0.5 M and 0.1 M) as in the NaCl case.

[13] Such transport may be difficult to maintain, as it may be restricted by the confining filters.

[14] We have discussed this study previously, in the blog post on transport limitations in confining filters.

[15] The actual measured bromide and calcium fluxes in Tinnacher et al. (2016) differ only by about one order of magnitude. This may be in part due to filter limitations, and in part be an indication that a completely homogeneous treatment begins to break down at these densities.

[16] Some results on bentonite samples of somewhat lower density are presented by Heister (2005).

Post-publication review: Tournassat and Steefel (2015), part V

This is the fifth part of the review of “Ionic Transport in Nano-Porous Clays with Consideration of Electrostatic Effects” (Tournassat and Steefel (2015) (referred to as TS15 in the following). For background and context please check the first part. Here we make some further comments on electric potentials, and this part is more of an appendix to the previous part, which discussed the proposed model for ion transport in bentonite.1 It is therefore recommended to read up on that part before continuing here.

Even worse problems?

So far in this review we have showed that the model for bentonite proposed in TS15 has both conceptual and mathematical flaws. In the previous part we showed that the “diffuse layer” flux is not derived adequately, although the resulting expression nevertheless can make some sense, if the involved electric potentials are completely reinterpreted. We also noted that the suggested “contributions” to the flux — related to “concentration gradients” (Fickian contribution) and “diffusion potential” (electric field contribution) — make no sense. These “contributions” can, however, be corrected, with the correction term (\(j^\star \)) involving the gradient of \(\Psi^\star\), the electric potential difference between bulk and “diffuse layer”

\begin{equation} j_\mathrm{conc,corr} = j_\mathrm{conc,TS15} + j^\star \end{equation}

\begin{equation} j_\mathrm{E,corr} = j_\mathrm{E,TS15} – j^\star \end{equation}

\begin{equation} j^\star \equiv \frac{A c_\mathrm{bulk} D_\mathrm{DL}zF}{RT} \nabla \Psi^\star \end{equation}

Here \(c_\mathrm{bulk}\) is the bulk water concentration of the considered species (generally a function of the model coordinate \(x\)), \(D_\mathrm{DL}\) is the corresponding diffusion coefficient in the “diffuse layer” domain, \(z\) the charge number, \(F\) the Faraday constant, and \(RT\) the usual thermal energy factor. \(A\) is what TS15 call the “DL enhancement factor”, and is given by (when correctly derived)

\begin{equation} A = e^{-\frac{zF}{RT}\Psi^\star} \end{equation}

The electric potentials in the bulk and “diffuse layer” domain are denoted \(\Psi_\mathrm{bulk}\) and \(\Psi_\mathrm{DL}\), respectively, giving the relation

\begin{equation} \Psi^\star \equiv \Psi_\mathrm{DL} – \Psi_\mathrm{bulk} \tag{1} \end{equation}

With these corrections and reinterpretations it may appear as if the transport model presented in TS15 actually has some merit. Here, however, I would like to point out what I think is a deeper problem, related to the assumption of requiring the various “porosity domains” to be locally in equilibrium. As in the previous part, we focus on the equilibrium between the bulk and “diffuse layer” domains.

In the previous part we did not pay detailed attention to the electric potential difference \(\Psi^\star\). TS15 suggest that \(\Psi^\star\) is to be calculated using the Donnan equilibrium framework. This makes some sense, from the perspective that the model assumes bulk water and “diffuse layer” to be in equilibrium locally, and in the provided examples \(\Psi^\star\) is calculated using the Donnan formula for a 1:1 electrolyte,2

\begin{equation} e^\frac{F\Psi^\star}{RT} = – \frac{q}{2c_\mathrm{bulk}} + \sqrt{\frac{q^2}{4c_\mathrm{bulk}^2} + 1} \end{equation} where \(q\) is a measure of the structural charge in the “diffuse layer”, in the examples set to \(q\) = 0.33 M.

But the requirement of Donnan equilibrium constrains the overall model quite heavily. For example, if we — as in the provided examples — impose concentration profiles in the bulk water, these determine the electric potential profile in this domain, via the Nernst-Planck framework. At the same time, the bulk water concentrations, together with a specified value of \(q\), also completely determine \(\Psi^\star\) (as a function of \(x\)), as well as all “diffuse layer” ion concentrations, via the Donnan equilibrium framework. These “diffuse layer” concentrations will, in turn, determine the electric potential (up to a constant) in this domain, via the Nernst-Planck framework.

But note that the electric potentials in the bulk and “diffuse layer” domains determined in this way, cannot in general also fulfill the requirement for Donnan equilibrium, i.e. \(\Psi^\star \neq \Psi_\mathrm{DL} – \Psi_\mathrm{bulk}\) (cmp. eq. 1). We consequently end up with a contraction, where we have used eq. 1 to determine the “diffuse layer” concentrations, while eq. 1 no longer applies after invoking the Nernst-Planck condition of requiring zero electric current!

This incompatibility is illustrated in the figure above.3 Note that there is nothing special with that the contradiction is here expressed in terms of \(\Psi_\mathrm{bulk}\) and \(\Psi_\mathrm{DL}\) not fulfilling eq. 1. Rather, all of the four relations illustrated in the above figure cannot in general be fulfilled simultaneously. In other words, as far as I see, for an imposed set of concentration profiles in one of the domains, the TS15 model cannot in general simultaneously fulfill these conditions

  • Zero electric current in the bulk water domain
  • Zero electric current in the “diffuse layer” domain
  • Donnan equilibrium between bulk and “diffuse layer”

This flaw is well illustrated in examples 2 and 3 in TS15 (the electric potential is the same for these cases). The electric potential in the bulk water can be calculated from the imposed concentration gradients of the main electrolyte, and is plotted in the figure below (blue line)

Here we have chosen \(\Psi_\mathrm{bulk}(0) = 0\) as reference point. (Note that TS15 are under the false impression that \(\Psi_\mathrm{bulk}(x)\) is zero everywhere.)

In the figure is also plotted the electric potential in the “diffuse layer”, as calculated from the Nernst-Planck framework and the concentration profiles in this domain (orange line; also shown here). As we noted previously, the reason for the much smaller variation of \(\Psi_\mathrm{DL}(x)\) as compared with \(\Psi_\mathrm{bulk}(x)\) is the ever-present counter-ions in the “diffuse layer” domain. Notice further that TS15, in contrast, are under the impression that the electrostatic potential in the “diffuse layer” equals the Donnan potential (red line), while such a profound potential variation is not supposed to affect the electromigration, for unclear reasons.

For \(\Psi_\mathrm{DL}\) we have here chosen the reference point \(\Psi_\mathrm{DL}(0) = \Psi^\star(0)\), i.e. we dictate that the bulk and “diffuse layer” potentials should differ by the Donnan potential when \(x=0\). But this is the only point where this condition is fulfilled: The difference between the electric potentials as calculated in this way (green line) does not at all resemble the Donnan potential!

At the moment, I don’t have the energy to sort out if there is any way out of this paradox, but my guess is that the problem is related to the very different treatments of transport in the \(x\)- and \(y\)-directions in the TS15 model. Remember that the assumption of equilibrium between all domains for a given \(x\)-value is equivalent to assuming infinite mobility in the \(y\)-direction of all species. To me, it appears as the TS15 model attempts to squeeze an intrinsically two-dimensional problem into a one-dimensional form. Note that this problem is not unique for the TS15 model, but arise in any multi-porous, multi-component description which assumes local equilibrium between domains (e.g. Appelo and Wersin (2007)).

Appendix: Comments on Figure 8

Disclaimer: This part functions as an appendix to this already appendix-like part of the review. The reader has been warned.

When considering electric potentials in the TS15 model I have developed a need to further comment on “Figure 8”. I alluded to this figure in the previous part of the review when discussing the term “Pseudo 2-D Cartesian system” used by TS15 (a term I don’t understand). “Figure 8” looks like this

The caption reads

Pseudo-2-D Cartesian system with diffusion along the x-axis and electrostatic potential developing along the y-axis due to the negative charge at clay mineral surfaces.

The only reference to this figure in the article text is at the beginning of the section presenting the transport model (Nernst-Planck equation in several “porosity domains”)

In the following, we will consider a pseudo 2-D Cartesian system in which diffusion takes place along the x axis only (Fig. 8).

The more I think about how this figure is included in the article, the more peculiar I find it. The figure clearly shows a quantitative result, as it displays specific values of a two-dimensional electric potential function outside “negative charge at clay mineral surface”. Yet, any real information about this calculation is nowhere to be found, which I see as a major problem (especially as the text is supposed to be a “fully developed text which can be used for self-study, research, or as a text-book for graduate-level courses”). Here we will first discuss information that is not provided, before continuing with speculating about how this figure may have been produced.

Stuff we are not told anything about

The ionic strength

Note that the “ionic strength gradient in bulk water” is only indicated in the figure, but not at all mentioned in the caption or in the text snippet that refers to this figure. Apparently, the physical situation depicted involves “bulk water” where the ionic strength, i.e. the concentration, falls of with the coordinate \(x\). But what is the actual value of this gradient? What are the boundary concentrations? And how has this ionic strength been used when calculating the electrostatic potential? We are not even told what type of electrolyte is considered! Is it 1:1?

Pseudo-2D and length scales

The only two places in TS15 where “Figure 8” is referenced are also the only places where the concept of a “pseudo-2D Cartesian system” is brought up. But what is meant by this term?! To me, it is more than a little strange to base an entire model presentation on a concept that is not really explained.

Taken at face value, the presented figure does not seem to show any “pseudo-“, but a “real” 2D coordinate system. Guided by the variation of the potential in the \(y\)-direction and by the clue that we are outside a negatively charged mineral surface, it’s quite safe to say that we — in the \(y\)-direction — are looking at an actual diffuse layer, as calculated e.g. in the Gouy-Chapman model. The \(y\)-dimension is thus reasonably on the nanometer scale. An ionic strength gradient, however, only makes sense as being on the macroscopic scale. The scale of the \(x\)-dimension is thus reasonably comparable to e.g. that of a laboratory sample, i.e. centimeters or even meters. Although “Figure 8” consequently is oddly scaled (the \(x\)-dimension is at least a million times larger than the \(y\)-dimension), I don’t see any reason to refer to it as “pseudo-2D”. It should of course be completely mandatory to state the length scales when presenting such a figure in a peer-reviewed scientific journal.

Further, we note that the figure seems to depict two copies facing each other of the same potential function, with one copy being flipped vertically (this explains the two oppositely pointing \(y\)-axes). Is this the reason for calling the coordinate system “Pseudo-2D”? But this manipulation (two copies facing each other) seems only to be for illustrative purposes, and has no impact on what actually seems to have been calculated.4

Surface charge

The potential is supposed to be “developing […] due to the negative charge at the clay mineral surfaces”. But then a surface charge density must have been specified when calculating the potential. Needless to say, this value should have been stated.

Stuff that is incorrect or odd

“Developing” potential

The figure caption “explains” that an electrostatic potential “develop[s] along the \(y\)-axis”. But what is significant is that an electric field is present near the surface. And, as we have brought up earlier, an electric field is equivalent to a potential variation. To put focus on that a potential “develops”, rather than that the potential varies, mostly misses the point. Also, why is the word “develop” used here? It suggests some sort of ongoing (unexplained) process.

“Developing” potential along the y-axis

A clear illustration that something is missed when speaking of potentials rather than fields is that TS15 implies that this “development” occurs exclusively in the \(y\)-direction. Note that an electric field is directed normal to the lines of constant electric potential. It is thus completely clear from the figure itself that the electric field is not in general directed solely in the \(y\)-direction!

This type of manifestation of an electric field in the \(x\)-direction relates to the complications discussed earlier of having a Donnan potential that varies with \(x\).

The potential has not resulted from solving the Nernst-Planck equation

Even though the purpose of the section is a treatment of the Nernst-Planck framework for transport of ionic charges, “Figure 8” does not present the result of a Nernst-Planck calculation. A main objective of such a treatment is to calculate the gradient of the electric potential (in the \(x\)-direction), i.e. an electric field. In contrast, the figure clearly shows that the potential is essentially zero, some distance away from the surface.

We note that a zero electric potential in the region that the authors presumably classify as “bulk” (some nanometers away from the surface…) is in line with the erroneous procedure of setting the bulk electric potential identically zero in the derivation of the Nernst-Planck flux.

What I think is illustrated/has been calculated

My guess is that “Figure 8” has been produced by calculating a set of potential profiles (in the \(y\)-direction) from the Gouy-Chapaman model for a given set of values of the ionic strength. These profiles have then simply been “stacked” in the \(x\)-direction. Below are some attempts to recreate “Figure 8” using this approach.

Since no parameters are given in TS15, we have to make several assumptions. For the figure below is assumed a linear profile of a 1:1 electrolyte that falls from 1000 mM at \(x=0\) to 100 mM at \(x=12\). We furthermore assume a surface charge density of 0.04 C/m2.

The magnitude and variation of this potential is comparable to that in “Figure 8” (the unit for the electric potential is here mV). We may note, however, that the linear ionic strength profile results in lines of constant electric potential that more and more “bends away” from the surface. This is because the Debye length in a 1:1 electrolyte depends on ionic strength as \(1/\sqrt{I}\). In contrast, the lines of constant potential in “Figure 8” are seen to “bend back”, suggesting a leveling off of the ionic strength, like this

If the present speculations are correct, why on earth has such a concentration profile been chosen? More broadly, why is this type of figure at all presented? If the authors simply want to show the quantities being calculated, why not provide a schematic illustration rather than a quantitative calculation? And if they want to present an actual calculation — which for unclear reasons does not involve the Nernst-Planck framework — they must of course provide a reader with enough information to make it understandable.

Update (260824): Part VI of this review is found here.

Footnotes

[1] As I have commented in the earlier parts: TS15 are fond of using the general terms “clays” and “clay minerals”, while it is clear that the publication mainly focus on systems with substantial ion exchange capacity and swelling properties. Here we will continue to use the term “bentonite” for these systems, and ignore the frequent references in TS15 to more general terms.

[2] TS15 are under the impression that \(\Psi_\mathrm{bulk}= 0\), and they believe that they calculate \(\Psi_\mathrm{DL}\) rather than \(\Psi^\star\). They furthermore emphasize that a more accurate calculation involves integrating the Poisson-Boltzmann equation (and claim falsely that the Donnan equilibrium “model” is based on an approximation of the Poisson-Boltzmann equation). But under all circumstances, the TS15 model requires specifying an arbitrary size of the “diffuse layer”; in the examples is used a width of 3 nm, which approximately corresponds to 3, 1 and 0.3 Debye lengths, respectively, for 1:1 background concentrations 0.1 M, 0.01 M, and 0.001 M. As we have discussed in the blog post on why multi-porosity models cannot be taken seriously, such a partitioning requires that a mechanism is provided for how equilibrium is maintained (i.e. why the “diffuse layer” does not occupy the full pore volume). TS15 do not supply such a mechanism, and neither does any other author promoting multi-porous models.

[3] In this figure, \( \left \{ c_{\mathrm{bulk} ,i} \right \} \) denotes the complete set of imposed concentrations in the bulk water domain. The corresponding set of species concentrations in the “diffuse layer” is denoted \(\left \{ c_{\mathrm{DL},i}\right \}\).

[4] Even if this choice is only “cosmetic”, I actually don’t understand why the figure is presented like this. Is it supposed to depict an (insanely long) interlayer? Perhaps with a “bulk water” phase squeezed in the middle? The same authors have actually presented such types of figures in later publications. As I discuss here, such a representation of a “bulk water” phase is nonsense.

Post-publication review: Tournassat and Steefel (2015), part IV

This is the fourth part of the review of “Ionic Transport in Nano-Porous Clays with Consideration of Electrostatic Effects” (Tournassat and Steefel, 2015) (referred to as TS15 in the following). For background and context please check the first part. This part covers the two sections “Constitutive equations for diffusion in bulk, diffuse layer, and interlayer water” and “Relative contributions of concentration, activity coefficient and diffusion potential gradients to total flux”.

“Constitutive equations for diffusion in bulk, diffuse layer, and interlayer water”

This section presents a mathematical formulation of ion diffusion in bentonite,1 based on the material descriptions in the earlier sections. As we have previously noted, these descriptions are fundamentally flawed in several respects. In particular, compacted bentonite is presented as consisting of stacks (called “particles”), where it is supposed to make sense to differ between external and internal interface water. TS15 also mean that compacted bentonite (sometimes?) is supposed to contain a bulk water phase.

As I have commented on in earlier parts, the only reason I can see to provide this nonsensical material description is as an attempt to to motivate a macroscopic, multi-porous model of bentonite. Here, TS15 make this claim quite explicit, as they write

Still it is possible to define three porosity domains, or water domains, that can be handled separately: the bulk water, the diffuse layer water and the interlayer water, the properties for which can be each defined independently.

This is in essence what I have referred to as “the mainstream view” of bentonite. It is basically “possible” to define anything, but the real question is if provided definitions are relevant and useful. And, as we have already discussed in detail, there is no rationale for introducing these “porosity domains” when modeling water saturated, compacted bentonite.2

Here we will first comment on the conceptual aspects of the provided mathematical description. Thereafter, we will delve into the mathematical formulations, as I’m quite convinced that these are not correct. Unfortunately, this latter part will be quite burdened with equations and notation, but for the motivated reader I think it may be worth going through.

Conceptual aspects

TS15 choose the Nernst-Planck description of ion diffusion, and begin by commenting that this is more rigorous than using Fick’s law. I certainly agree with that a general description of ion diffusion in bentonite requires treating electrostatic couplings between the various system components (TOT-layers, ions). I don’t think, however, that putting up a massively complex description of multi-component diffusion in “three porosity domains” is the appropriate starting point for including such couplings. Since we have every reason to believe that e.g. no bulk water phase is present, I mean that this type of treatment only leads us astray from understanding the actual processes involved (we will return to this aspect in later parts of the review).

Also, as the “Fickian” aspect was the focus of the earlier section on diffusion, a reader of TS15 could here understandably get the impression that a Nernst-Planck treatment will “fix” the “issues” addressed there. But, as we have already discussed in some detail, the shortcomings of the traditional sorption-diffusion model are not solved by including multi-component diffusion in a bulk water phase. They are solved by removing the bulk water phase.

Although the above quotation states that the various “porosity domains” can be handled separately, and that their properties can be defined independently, this is not what is done in TS15. Rather, the treatment of any “porosity domain” assumes equilibrium with the corresponding bulk water phase. The entire description in TS15 is thus fully centered around the bulk water phase.

TS15 insist on treating their model quantities as functions of two spatial coordinates (\(x\) and \(y\)), in what they refer to as a “pseudo 2-D Cartesian system” (I don’t fully understand what that means). Diffusive flux is only assumed to take place in the \(x\)-direction, while the “\(y\)”-dimension is used for stacking the different “porosity domains”. The description can be schematically illustrated like this

Here we have for illustrative purposes discretized the various components in \(x\)- and \(y\)-directions. The bulk water domain is colored blue, the “interlayer” domain pink, and the “diffuse layer” domain green. For a given \(x\)-position, the “diffuse layer” and the “interlayer” domains are assumed to always be in equilibrium with the corresponding bulk water phase. TS15 nowhere consider the length scale in the \(y\)-direction (is it therefore the coordinate system is referred to as “pseudo 2-D”?), which in practice makes the model a collection of 1-dimensional domains that are in equilibrium locally. Note that even though diffusion only is accounted for in the \(x\)-direction, transport occurs also in the \(y\)-direction, as a consequence of equilibration between the “porosity domains”.

This description is exactly what we have investigated in the blog post on why multi-porous models cannot be taken seriously. To summarize what was said there, without properly defining the length scales, it makes no sense to “short-cut” the model in the \(y\)-direction (to assume equilibrium for all domains at the same \(x\) is in a sense equivalent to assuming infinitely high mobility of all components in the \(y\)-direction). And even if we assume that such an assumption is valid — which would mean that we consider a thin strip of stacked parallel domains, where the extension in \(y\) is negligible in comparison to the extension in \(x\) — the resulting model has really nothing to do with actual bentonite. As we concluded in the multi-porosity blog post, the only way to make sense of this type of description is as a set of macroscopic continua that are assumed to be locally in equilibrium. How this equilibrium is supposed to be maintained has never been suggested by any proponent of this description. Note that this description (in particular the existence of a bulk water phase in equilibrium) disqualifies the model for describing swelling and swelling pressure.

Incorrect application of the Nernst-Planck framework

While the presented model makes little sense conceptually, TS15 also fail in applying the Nernst-Planck framework. The problem arises, as far as I can see, from that they don’t fully recognize the role of the electric potential.

As we now begin scrutinizing the details of the formulation, we will suppress the variables \(x\) and \(y\) in order to, hopefully, make the equations a little more readable. It should be understood that any quantity is evaluated for some specific value of \(x\), and that all “porosity domains” are supposed to be in equilibrium at the same value of \(x\).

The electro-chemical potential

In most standard thermodynamic text books we learn that the chemical potential governs the equilibrium associated with mass transfer. Just as e.g. pressure and temperature (which govern mechanical and thermal equilibrium, respectively), the chemical potential is defined by a specific derivative of a thermodynamic potential, e.g.3

\begin{equation} \bar{\mu} = \left ( \frac{\partial G}{\partial n} \right )_{T,p} \tag{1} \end{equation}

where \(G\) is the Gibbs free energy, \(n\) number of moles, \(p\) pressure, and \(T\) temperature. The corresponding mass flux is generally written

\begin{equation} j = -\frac{cD}{RT}\nabla \bar{\mu} \tag{2} \end{equation}

where \(c\) is concentration, \(D\) the diffusion coefficient4, and \(RT\) the usual absolute temperature factor. Here, and in the following, we use the symbol \(\nabla\), which denotes the general gradient operator, but since the model is effectively one-dimensional, it can simply be seen as a neat way of writing \(\partial/\partial x\).

For charged species, it is common to refer to the quantity defined in eq. 1 as the electro-chemical potential, and write it as composed of an “ordinary” and a purely “electrical” part

\begin{equation} \bar{\mu} = \mu + zF\Psi \tag{3} \end{equation}

where \(z\) is the charge number of the considered species, \(F\) is the Faraday constant and \(\Psi\) is the electric potential. The “ordinary” chemical potential \(\mu\) (without bar) is, perhaps a bit confusingly, also often referred to as the chemical potential. I will here continue to refer to this part as “ordinary”. The “ordinary” chemical potential is furthermore conventionally expressed in terms of a reference potential (\(\mu^0\)) and an activity \(a\)

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

A lot can be said about the decomposition in eq. 3, but it is clear that singling out an electric potential term is useful in e.g. electrochemistry or for describing charged clay. It should, however, be kept in mind that mass transfer is fundamentally governed by gradients in \(\bar{\mu}\); always keeping eqs. 1 and 2 in mind will avoid us from making mistakes, because the mass transfer rate relates to the “total” (i.e. electrochemical) potential, and for charge neutral species the description reduces to gradients in the “ordinary” chemical potential.

Nernst-Planck flux

Combining eqs. 3 and 4 gives

\begin{equation} \bar{\mu} = \mu^0+RT \ln a + zF\Psi \tag{5} \end{equation}

with the corresponding flux (eq. 2)

\begin{equation} j = -cD\nabla \ln a -\frac{cDzF}{RT} \nabla \Psi \end{equation}

Expressing the activity in terms of an activity coefficient, \(a = \gamma c\), the flux can also be written (TS15 are quite fond of including activity coefficients explicitly)

\begin{equation} j = -D\nabla c -cD\nabla \ln \gamma -\frac{cDzF}{RT} \nabla \Psi \tag{6} \end{equation}

Considering an arbitrary set of diffusing charged species (using the index \(i\)), and utilizing that the electric current is zero, lead to an expression for the electric potential gradient

\begin{equation} \nabla \Psi = -\frac{RT}{F}\frac{\sum z_iD_i \left ( \nabla c_i + c_i \nabla \ln \gamma_i \right )}{\sum z_i^2 D_i c_i} \tag{7} \end{equation}

Misunderstanding the electric potential

For the bulk water phase, TS15 indeed provide an expression for the flux that is essentially the same as eq. 6 (their eq. 37), and which they refer to as the Nernst-Planck equation. They claim, however, that the electrochemical potential in this case lack an electric potential term (my emphasis)5,6

In absence of an external electric potential, the electrochemical potential in the bulk water can be expressed as (Ben-Yaakov 1981; Lasaga 1981) \begin{equation} \bar{\mu}_\mathrm{bulk} = \mu^0+RT \ln a_\mathrm{bulk} \end{equation}

But even without an externally applied electric field,7 a zero bulk electric potential cannot be assumed, of course, if the goal is to treat individual ion mobilities; as just shown, the gradient in electric potential that appears in eq. 6 is a result of a corresponding term in the electrochemical potential (eq. 5). Oddly, TS15 seem to treat the electric potential term in the flux as a quantity unrelated to the electrochemical potential, giving it a separate symbol, \(^\mathrm{b}\Psi_\mathrm{diff}\), and writing

\(^\mathrm{b}\Psi_\mathrm{diff}\) is the diffusion potential that arises because of the diffusion of charged species at different rates.

It may be natural for a reader at this point to simply assume that TS15 have missed writing out the term \(zF^\mathrm{b}\Psi_\mathrm{diff}\) when stating the electrochemical potential. But this seems to be a genuine misunderstanding rather than a mistake/typo, because the pattern repeats in the derivation of the flux in the other “porosity domains”.

For e.g. the “diffuse layer”,8 TS15 recognize the presence of an electric potential in the expression for the electrochemical potential, writing it (this is more awkwardly expressed in eq. 42 in TS15)

\begin{equation} \bar{\mu}_\mathrm{DL} = \mu^0 + RT\ln a_\mathrm{DL} + zF\Psi_\mathrm{DL} \tag{8} \end{equation}

where index “DL” refers to quantities in the “diffuse layer”.

However, the corresponding flux expression contains a different potential, labelled \(^\mathrm{DL}\Psi_\mathrm{diff}\) (eq. 46 in TS15)

\begin{equation} j_\mathrm{DL} = -\frac{c_\mathrm{DL}D_\mathrm{DL}}{RT}\nabla \bar{\mu}_\mathrm{DL} -\frac{c_\mathrm{DL}D_\mathrm{DL}zF}{RT} \nabla ^{\mathrm{DL}}\Psi_\mathrm{diff} \tag{9} \end{equation}

TS15 don’t further comment what \(^{\mathrm{DL}}\Psi_\mathrm{diff}\) is supposed to represent, but it must reasonably be understood as “the diffusion potential that arises because of the diffusion of charged species at different rates”, in analogy with what was claimed for the bulk water phase. Note that when eq. 8 is combined with eq. 9, the flux expression contains two different electric potential gradients! (TS15 never address this oddity.)

It is thus quite clear that TS15 misunderstand the function of the electric potential in the Nernst-Planck framework. When presenting the expression for the “diffuse layer” flux (eq. 9), they also refer to Appelo and Wersin (2007), who, in turn, express the misconception explicitly9

The gradient of the electrical potential [in the expression for the flux] originates from different transport velocities of ions, which creates charge and an associated potential. This electrical potential may differ from the one used in [the expression for the electro-chemical potential], which comes from a charged surface and is fixed, without inducing electrical current.

I cannot understand this passage in any other way than that Appelo and Wersin (2007) are under the impression that different electric potentials can simultaneously act independently in a given point. And it seems like TS15 are under some similar impression.

This ignorance leads to more errors in the description of the “diffuse layer” in TS15. We should remember that the promoted model requires the “diffuse layer” and bulk water domains to be in equilibrium (for the same coordinate value \(x\)). When TS15 express this condition, i.e. \(\bar{\mu}_\mathrm{DL} = \bar{\mu}_\mathrm{bulk}\), they again leave out the electric potential in the bulk water (eq. 42 in TS 15)

\begin{equation} \mu^0 + RT\ln a_\mathrm{DL} + zF\Psi_\mathrm{DL} = \mu^0 + RT\ln a_\mathrm{bulk}\;\:\;\;\;\;\;\mathrm{(WRONG)} \tag{10} \end{equation}

eq. 10 can be rewritten

\begin{equation} a_\mathrm{DL} = a_\mathrm{bulk}\cdot e^{-\frac{zF}{RT}\Psi_\mathrm{DL}} \;\:\;\;\;\;\;\mathrm{(WRONG)} \tag{11} \end{equation}

TS15 utilize a simplified version of eq. 11, expressed in terms of concentrations rather than activities, by assuming identical activity coefficients in the two domains10

\begin{equation} c_\mathrm{DL} = c_\mathrm{bulk}\cdot e^{-\frac{zF}{RT}\Psi_\mathrm{DL}} \;\:\;\;\;\;\;\mathrm{(WRONG)} \tag{12} \end{equation}

Note that the exponential in eqs. 11 and 12 actually should contain the electric potential difference between “diffuse layer” and bulk (see below).

As TS15 have not included any electric potential in the bulk water phase, they continue by incorrectly substituting \(RT\nabla \ln a_\mathrm{bulk}\) for \(\nabla\bar{\mu}_\mathrm{DL}\) in eq. 9 (i.e. they use the incorrect relation in eq. 10), giving (TS15 eq. 47)

\begin{equation} j_\mathrm{DL} = -c_\mathrm{DL}D_\mathrm{DL}\nabla \ln a_\mathrm{bulk} -\frac{c_\mathrm{DL}D_\mathrm{DL}zF}{RT} \nabla ^{\mathrm{DL}}\Psi_\mathrm{diff} \;\;\;\;\mathrm{(WRONG)} \tag{13} \end{equation}

Note that this additional error “solves” the earlier pointed out problem of having two electric potential gradients.

By utilizing the requirement of zero electric current, eq. 13 gives

\begin{equation} \nabla ^{\mathrm{DL}}\Psi_\mathrm{diff} = -\frac{RT}{F}\frac{\sum z_iD_{\mathrm{DL},i}c_{\mathrm{DL},i} \nabla \ln a_{\mathrm{bulk},i} }{\sum z_i^2 D_{\mathrm{DL},i} c_{\mathrm{DL},i}} \;\:\;\;\;\;\;\mathrm{(WRONG)} \tag{14} \end{equation}

By substituting eq. 12 into this expression, we end up with the formula for the gradient of the mysterious potential \(^{\mathrm{DL}}\Psi_\mathrm{diff}\) (TS15 eq. 48)

\begin{equation} \nabla ^{\mathrm{DL}}\Psi_\mathrm{diff} = \end{equation} \begin{equation} -\frac{RT}{F}\frac{\sum z_iD_{\mathrm{DL},i} e^{-\frac{zF}{RT}\Psi_\mathrm{DL}} \left ( \nabla c_{\mathrm{bulk},i} + c_{\mathrm{bulk},i}\nabla \ln \gamma_{\mathrm{bulk},i} \right )} {\sum z_i^2 D_{\mathrm{DL},i} e^{-\frac{zF}{RT}\Psi_\mathrm{DL}}c_{\mathrm{bulk},i}} \;\mathrm{(WRONG)} \tag{15} \end{equation}

At face value, eq. 15 is a quite weirdly looking equation, as it relates two electric potentials — \(^{\mathrm{DL}}\Psi_\mathrm{diff}\) and \(\Psi_\mathrm{DL}\) — that both are supposed to be associated with the “diffuse layer”. But, as we will see below, there is actually a way to make some sense of eq. 15, by completely reinterpreting what these potentials represent.

A “correct” formulation

Update (260803): The electrostatic potential within the homogeneous mixture model is treated here.

Most of the errors pointed out above are corrected by including the electric potential in the bulk water and writing the condition for equilibrium as (compare eq. 10)

\begin{equation} \mu^0 + RT\ln a_\mathrm{DL} + zF\Psi_\mathrm{DL} = \mu^0 + RT\ln a_\mathrm{bulk} + zF\Psi_\mathrm{bulk} \tag{16} \end{equation}

Writing the electric potential difference between “diffuse layer” and bulk water as11

\begin{equation} \Psi^\star \equiv \Psi_\mathrm{DL} – \Psi_\mathrm{bulk} \tag{17} \end{equation}

eq. 16 can be rewritten

\begin{equation} a_\mathrm{DL} = a_\mathrm{bulk}\cdot e^{-\frac{zF}{RT}\Psi^\star} \tag{18} \end{equation}

Note that, when correctly derived, eq. 18 naturally contains the difference in electric potential between “diffuse layer” and bulk.

The flux in the “diffuse layer” is (eq. 6)

\begin{equation} j_\mathrm{DL} = -c_\mathrm{DL} D_\mathrm{DL} \nabla \ln a_\mathrm{DL} – \frac{c_\mathrm{DL} D_\mathrm{DL} zF}{RT} \nabla \Psi_\mathrm{DL} \tag{19} \end{equation}

But if we now plug in eq. 18 in eq. 19 we of course get

\begin{equation} j_\mathrm{DL} = -c_\mathrm{DL} D_\mathrm{DL} \nabla \ln a_\mathrm{bulk} + \frac{c_\mathrm{DL} D_\mathrm{DL} zF}{RT} \nabla \Psi^\star – \frac{c_\mathrm{DL} D_\mathrm{DL} zF}{RT} \nabla \Psi_\mathrm{DL}, \end{equation} which can be simplified to \begin{equation} j_\mathrm{DL} = -c_\mathrm{DL} D_\mathrm{DL} \nabla \ln a_\mathrm{bulk} – \frac{c_\mathrm{DL} D_\mathrm{DL} Fz}{RT} \nabla \Psi_\mathrm{bulk}, \tag{20} \end{equation} and, by identifying the electro-chemical potential in the bulk \begin{equation} j_\mathrm{DL} = -\frac{c_\mathrm{DL} D_\mathrm{DL}}{RT} \left ( RT \nabla \ln a_\mathrm{bulk} + Fz \nabla \Psi_\mathrm{bulk} \right ) = -\frac{c_\mathrm{DL} D_\mathrm{DL}}{RT} \nabla \bar{\mu}_\mathrm{bulk} \end{equation}

This whole “derivation” leads back to the rather trivial result that the flux in the diffuse layer is given by eq. 2, which we could have written down from the start! (because the model assumes \(\bar{\mu}_\mathrm{bulk} = \bar{\mu}_\mathrm{DL}\); eq. 16)

As TS15 have established the expression for the gradient of the electrochemical potential in the bulk water phase (which is implicit in their eq. 40), there should strictly be no need to consider a new expression for the same quantity in any other phase. Rather, they could simply have used the bulk water expression in all “porosity domains”, as a consequence of the assumption that these are all supposed to be in equilibrium. In a sense, this is actually what is done in TS15 — mainly by chance! — by establishing eq. 13 (their eq. 47).

Comparing with eq. 20, we see that the incorrect eq. 13 can be “saved” by reinterpreting \(^{\mathrm{DL}}\Psi_\mathrm{diff}\) as \(\Psi_\mathrm{bulk}\). Similarly, as TS15 assume the bulk electric potential to be zero, eq. 15 can be “saved” by also reinterpreting \(\Psi_\mathrm{DL}\) as \(\Psi^\star\) in that expression.12 I find this quite hilarious: By making several errors in its derivation, eq. 15 is in a sense a correct expression for the electric potential gradient in the bulk water — a potential that TS15 has put identically equal to zero.

But even if the total flux in the “diffuse layer” is correctly given by combining eqs. 12, 13 and 15 (and by completely ignoring what TS15 mean \(\Psi_\mathrm{DL}\) and \(^{\mathrm{DL}}\Psi_\mathrm{diff}\) represent), TS15 continue by defining the separate terms in eq. 13 as contributions from the “concentration gradient”, and the “diffusion potential”. As we will explore next, this interpretation fails miserably.

“Relative contributions of concentration, activity coefficient and diffusion potential gradients to total flux”

According to TS15, the “concentration gradient” and the “diffusion potential” contributions to the “diffuse layer” flux are given by, respectively (TS15 eq. 50 and below)

\begin{equation} j_\mathrm{conc,TS15} = -D_\mathrm{DL}A\nabla c_\mathrm{bulk} \tag{21} \end{equation}

and

\begin{equation} j_\mathrm{E,TS15} = zD_\mathrm{DL}c_\mathrm{bulk}A\frac{\sum z_iD_{\mathrm{DL},i}A_i\left ( \nabla c_{\mathrm{bulk},i} + c_{\mathrm{bulk},i} \nabla \ln \gamma_{\mathrm{bulk},i} \right )} {\sum z_i^2 D_{\mathrm{DL},i} A_i c_{\mathrm{bulk},i}} \tag{22} \end{equation}

Here we use the index “conc” for the “concentration gradient” contribution, and “E” for the “diffusion potential” contribution. \(A\) is referred to as a “DL enrichment factor”, and is essentially defined as the concentration ratio \(c_\mathrm{DL}/c_\mathrm{bulk}\). Using the incorrect relation in eq. 12, TS15 write these as \(A = e^{-\frac{zF}{RT}\Psi_\mathrm{DL}}\), but, as we see from eq. 18, they are really given by13 (we continue assuming identical activity coefficients in the two domains)

\begin{equation} A = e^{-\frac{zF}{RT}\Psi^\star} \tag{23} \end{equation}

TS15 also define a third contribution, related to the gradient of the bulk water activity coefficient. Here we will not further discuss this contribution, as it does not give any additional insight. Moreover, since TS15 anyway derive their model under the unjustified assumption that activity coefficients in the “diffuse layer” and the bulk water are identical, I cannot see the use of including their spatial variation in the description.14 (TS15 spend a couple of pages on activity coefficient models that we will ignore.)

Examples

To explore the various couplings in the presented model, TS15 apply the Nernst-Planck framework in three examples. We can see immediately from the presented graphs that their partitioning of the total flux in “concentration gradient” and “diffusion potential” contributions makes no sense.

“Example 2” imposes constant concentration gradients in the bulk water of NaCl and corresponding \(^{22}\mathrm{Na}^+\) and \(^{36}\mathrm{Cl}^-\) tracers; the NaCl concentration drops from 0.1 M to 0.001 M, and the tracer concentrations drop from 10-9 M to 10-11 M (domain length is 10 mm).

The corresponding sodium and chloride tracer concentrations in the “diffuse layer” look like this15

These profiles make sense: bulk water ionic strength decreases with distance, but so do the tracer concentrations. For the process of accumulating \(^{22}\mathrm{Na}^+\) in the diffuse layer, these two effects oppose each other, resulting in a quite flat profile. We thus expect the corresponding “concentration gradient” contribution to the flux to be quite moderate, and to fall off with distance (as the profile flattens with distances). The corresponding flux graph presented in TS15, however, looks completely different16

This plot makes no sense: The “concentration gradient” contribution is seen to increase quite dramatically with distance, rather than falling off. The value of this contribution is also orders of magnitude too large, given the imposed sodium diffusion coefficient of 1.33⋅10-10 m2/s. Moreover, the “concentration gradient” contribution is “compensated” by an equally nonsensical “diffusion potential” contribution. Note, for instance, that the “diffusion potential” contribution is negative, which implies that the corresponding electric field is supposed to be directed towards higher concentrations. This can certainly not be the case, as the electric potential gradient is caused by the negative ion having higher mobility than the positive ion (chloride diffusivity is set to 2.03⋅10-10 m2/s).

In “example 3”, the tracer concentrations in the bulk is set to a constant value (1⋅10-9 M), while the same concentration gradient as in “example 2” is maintained for the main NaCl electrolyte (from 0.1 M to 0.001 M). We thereby expect the corresponding \(^{22}\mathrm{Na}^+\) concentration in the “diffuse layer” to strongly increase with distance, which is also what is presented in TS15

while the corresponding flux plot looks like this16

This plot is almost comically absurd. According to TS15, the highly skewed concentration profile above is supposed to give no (zero, nil, 0) contribution to the flux (we see from eq. 21 that this is a consequence of that this “contribution” is directly proportional to the concentration gradient in the bulk). Instead, the huge flux is supposed to be caused entirely by an electric field that has the wrong direction! I can’t even really begin to imagine how these two plots have ended up next to each other in a peer-reviewed published article.

Note that the flux associated with a concentration gradient is what we may reasonably call a “Fickian” contribution. If TS15 mean (and they do) that these examples demonstrate how ion diffusion in bentonite works, we can understand the focus on the “Fickian” aspect at the beginning of the article (covered here). But the only reasonable response to these outlandish results is that they demonstrate that the definitions of eqs. 21 and 22 simply make no sense.

The real concentration gradient and electric field contributions

The only reasonable way to define “concentration gradient” and “diffusion potential” contributions to the “diffuse layer” flux is as the two terms in eq. 19, respectively. To rewrite these, we utilize eq. 16 (or 18), giving for the “concentration gradient” contribution (we continue ignoring activity coefficients)

\begin{equation} j_\mathrm{conc, corr} = -c_\mathrm{DL} D_\mathrm{DL} \nabla \ln c_\mathrm{DL} = \end{equation} \begin{equation} -c_\mathrm{DL} D_\mathrm{DL} \nabla \ln c_\mathrm{bulk} + \frac{c_\mathrm{DL} D_\mathrm{DL}zF}{RT} \nabla \Psi^\star = \end{equation} \begin{equation} -D_\mathrm{DL} A \nabla c_\mathrm{bulk} + \frac{A c_\mathrm{bulk} D_\mathrm{DL}zF}{RT} \nabla \Psi^\star = j_\mathrm{conc, TS15} + j^\star \tag{24} \end{equation}

where we have defined

\begin{equation} j^\star \equiv \frac{A c_\mathrm{bulk} D_\mathrm{DL}zF}{RT} \nabla \Psi^\star. \tag{25} \end{equation}

In the same manner, the correct “diffusion potential” contribution is

\begin{equation} j_\mathrm{E, corr} = – \frac{c_\mathrm{DL} D_\mathrm{DL} zF}{RT} \nabla \Psi_\mathrm{DL} = \end{equation} \begin{equation} – \frac{D_\mathrm{DL}Ac_\mathrm{bulk} zF}{RT} \nabla \Psi_\mathrm{bulk} – \frac{D_\mathrm{DL}Ac_\mathrm{bulk} zF}{RT} \nabla \Psi^\star = \end{equation} \begin{equation} D_\mathrm{DL}Ac_\mathrm{bulk} z \frac{\sum z_iD_iA_i\left ( \nabla c_{\mathrm{bulk},i} + \nabla c_{\mathrm{bulk},i} \ln \gamma_{\mathrm{bulk},i} \right )} {\sum z_i^2 D_i A_i c_{\mathrm{bulk},i}} – \frac{D_\mathrm{DL}Ac_\mathrm{bulk} zF}{RT} \nabla \Psi^\star = \end{equation} \begin{equation} j_\mathrm{E, TS15} – j^\star \tag{26} \end{equation}

where we have utilized that \(\nabla \Psi_\mathrm{bulk}\) is actually what is expressed in eq. 15 (where \(\Psi_\mathrm{DL}\) should be replaced by \(\Psi^\star\)).

We note that, to compensate the nonsensical expressions given in TS15, we should add the term \(j^\star\) (eq. 25) to the “concentration concentration” contribution (eq. 21), and subtract the same term from the “diffusion potential” contribution (eq. 22). Making these corrections gives the following components of the tracer fluxes in “example 2”

This is an infinitely more reasonable situation than what is depicted in TS15. Although the sodium flux has a non-negligible contribution from the electric field, the larger contribution is still from the concentration gradient (and none of these are gigantic terms that cancel). The concentration contribution also falls off with distance, in accordance with the shape of the concentration profile.

For chloride, the field contribution to the flux is negligible, i.e. this flux is essentially fully governed by the concentration gradient. The electric field contributions for both ions are also seen to have the correct signs: the electric field is directed from high to low concentration, and mainly functions to boost the sodium transport, in order to “keep up” with the faster chloride ions.

For “example 3” we get the following picture

The corrected \(^{22}\mathrm{Na}^+\) flux is essentially fully due to the concentration gradient, in absolute contrast to what is concluded in TS15, who mean that this flux is completely governed by an electric field in the wrong direction. Also the \(^{36}\mathrm{Cl}^-\) transport is basically solely governed by the concentration gradient, rather than by an incorrectly directed electric field (as stated in TS15). In conclusion, most of the “diffuse layer” diffusion in these examples can actually be classified as “Fickian”.

We may also investigate the electric potential profile in the “diffuse layer” in both of these examples (this is the same in the two cases, as the main electrolyte distribution does not change)

Here we have chosen the reference \(\Psi_\mathrm{DL}(0) = 0\). The total potential drop is only about 1 mV. Such a relatively small drop is reasonable because the denominator in the Nernst-Planck expression for the electric potential gradient (eq. 7) will always be large due to the ever-present counter-ions in the “diffuse layer”. The electric potential gradient — and thus the corresponding electric potential drop — is therefore suppressed. Physically, this means that since many (equally charged) charge carriers are always present, smaller potential differences are required to cancel electric currents caused by differences in mobility (a “diffuse layer” is a quite good conductor).

Even worse problems?

Even though some sense can be made out of the derived expression for the flux in the “diffuse layer” domain — by completely reinterpreting the electric potentials involved — it seems as the overall model is too constrained. Specifically, for an imposed set of concentration profiles in one domain it is not possible, as far as I can see, to simultaneously have zero current in all domains, while also maintaining (Donnan) equilibrium. As this blog post is already quite massive, I will elaborate on this point in the next part of the review.

Summary

Here is an attempt to sum up the main messages of this blog post.

  • Conceptually, the clay model presented in TS15 is exactly what was discussed in the blog post on multi-porous models, and the same issues that are identified there are present here. In particular, no attention is paid to length scales (perhaps that is why TS15 call the coordinate system “pseudo-2D”…), and no mechanism whatsoever is suggested for how the different diffusing domains are supposed to maintain equilibrium.
  • Mathematically (or perhaps physically), the presented Nernst-Planck flux expressions are incorrectly derived. The source of the error, as far as I can see, appears to be a misunderstanding of how electric potentials function.
  • TS15 define “contributions” to the “diffuse layer” flux, claimed to be related to the concentration gradient and the “diffusion potential” (i.e. the electric field), respectively. It is, however, quite obvious that these “contributions” are completely nonsensical: highly skewed concentration profiles are claimed to not have any concentration gradient contributions, and several “diffusion potential” “contributions” have the electric field in the wrong direction. We have shown that these “contributions” can be corrected, where the correction term involves the gradient of the Donnan potential. With these corrections, fluxes in the provided examples must be interpreted completely differently (they’re basically “Fickian”).
  • As far as I can see, the proposed model has even larger problems, related to the imposed Donnan equilibrium. We will address this issue in the next part.

Update (260526): Part V of this review is found here.

Footnotes

[1] As I have commented in the earlier parts: TS15 are fond of using the general terms “clays” and “clay minerals”, while it is clear that the publication mainly focus on systems with substantial ion exchange capacity and swelling properties. Here we will continue to use the term “bentonite” for these systems, and ignore the frequent references in TS15 to more general terms.

[2] It is of course crucial to include a component that represents compartments where the exchangeable ions reside. This is done in the TS15 model by both the “diffuse layer water” and the “interlayer water” domains. But the distinction made between these domains is based on the flawed “stack” concept.

[3] This equation assumes a single component. The formulation of the Nernst-Planck framework naturally involves several different charged species. When several species are involved, we will indicate this with an index \(i\) in the equations.

[4] In some of their equations, TS15 use (electrical) mobility, \(u\), rather than diffusivity, \(D\). These quantities are related via the Einstein relation \(D = uRT/(F|z|)\). I don’t see the point in involving \(u\), as it typically makes expressions even more cluttered, and since we here ultimately are interested in diffusion coefficients anyway.

[5] In order to not cause too much confusion, and to try to simplify a bit, I use slightly different mathematical notation than what is actually used in the quotation. In particular, I use the notation \(\bar{\mu}\) for the electro-chemical potential, while TS15 don’t use a bar (\(\mu\)). I also try to avoid the index \(i\) as much as possible.

[6] Fun fact: this statement is nowhere found in neither (Ben-Yaakov, 1981) nor (Lasaga, 1981) (at least I can’t find any).

[7] I whined about electrostatics being poorly understood in the bentonite research field in an earlier part of this review, but here is more fuel for my argument. The statement “absence of an [external] potential” has no physical meaning, as we are free to choose the reference point (the absolute value of a potential has no physical meaning). What TS15 must mean in the quote is “the absence of an external electric field”. The electric field relates to the potential as \(E = -\nabla \Psi\). Thus, all gradients of electric potentials that occur in this text are synonymous with electric fields (electric fields drive electric currents).

[8] This post focus almost entirely on the “diffuse layer” domain, but a similar analysis can be made for the “interlayer” domain. This is left as an exercise for the reader.

[9] It should of course also rather read “…which creates a charge separation and an associated potential gradient.”, or simply “…which induces an electric field.” (showing that this part of the sentence is redundant). See also footnote 7.

[10] TS15 write cryptically that equating the activity coefficients (and the reference potentials) in bulk and “diffuse layer” is assumed “by following the [Modified Gouy-Chapman] model”. But I don’t see why this model has to be alluded to here, these assumptions can just be made.

[11] Yes, this is a Donnan potential. We will discuss this more in the next part part of the review. Update (260526): Part V of this review is found here.

[12] Again, this is related to Donnan equilibrium between the bulk and “diffuse layer” domains, that we will discuss further in the next part. Update (260526): Part V of this review is found here.

[13] This is \(f_D^{-z} \), where \(f_D\) is the Donnan factor.

[14] Rather, I would argue for that the activity coefficients in a “diffuse layer” domain will be quite insensitive to the imposed external (bulk) concentration, for details see Birgersson (2017).

[15] In producing these graphs we have used the Donnan equilibrium framework to calculate the “diffuse layer” concentrations. These are given from eq. 23, where \(\Psi^\star\) is calculated from

\begin{equation} f_D = e^\frac{F\Psi^\star}{RT} = – \frac{q}{2c_\mathrm{bulk}} +
\sqrt{\frac{q^2}{4c_\mathrm{bulk}^2} + 1} \end{equation}

where \(q\) is a measure of the structural charge in the “diffuse layer”, in the examples set to \(q\) = 0.33 M.

[16] Note that I have not included activity coefficient gradients when producing the plots in this section. They may therefore differ slightly from the published plots. This does not in any way influence the conclusions drawn here.