• Part I
• Part II
• Part III
• Part IV

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.¹ 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” (Fickan 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,²

\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.³ 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.⁴

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/m².

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.

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.

• Part I
• Part II
• Part III

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,¹ 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.²

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.³

\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 coefficient⁴, 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,⁷ 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”,⁸ 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 explicitly⁹

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 domains¹⁰

\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

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 as¹¹

\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.¹² 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 by¹³ (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.¹⁴ (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 this¹⁵

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 different¹⁶

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 m²/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 m²/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 this¹⁶

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.
  

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.

[12] Again, this is related to Donnan equilibrium between the bulk and “diffuse layer” domains, that we will discuss further in the next part.

[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.