Category Archives: Numerical modeling

Solving the Nernst-Planck equations

I have started developing a tool that numerically solves the Nernst-Planck equations in bentonite on a 1D grid. It treats the bentonite within the homogeneous mixture model, accounting for the structural clay charge via the parameter \(c_\mathrm{IL}\). It also handles external solution reservoirs of finite size, and Donnan equilibrium is maintained dynamically at the interfaces between clay and external solutions. At the moment, the tool does not account for activity coefficients.

The expression for the flux of a specific species in the Nernst-Planck framework is

\begin{equation} j = j_\mathrm{conc} + j_\mathrm{E}\tag{1} \end{equation}

with the concentration concentration contribution (the “Fickan” part)

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

and the electromigration contribution (driven by the induced electric field)

\begin{equation} j_\mathrm{E} = -\phi\cdot z\cdot c\cdot D_c\cdot \nabla\tilde{\psi} \tag{3} \end{equation}

Here \(D_c\), \(z\) and \(c^\mathrm{int}\) are the diffusion coefficient, charge number and concentration, respectively, for the considered species. \(\phi\) denotes porosity. For convenience, we denote the spatial derivative (\(\partial/\partial x\)) by \(\nabla\). \(\tilde{\psi}\) is the dimensionless electric potential,1 whose gradient is

\begin{equation} \nabla\tilde{\psi} = -\frac{\sum z_i\cdot D_{c,i}\cdot \nabla c^\mathrm{int}_i} {\sum z_i^2\cdot D_{c,i}\cdot c^\mathrm{int}_i }\tag{4} \end{equation}

where the summations are over all involved charged species. Eq. 4 follows from requiring zero electric current.

The developed tool uses eqs. 1 — 4 and solves the corresponding continuity equations2

\begin{equation} \phi\cdot\frac{\partial c^\mathrm{int}_i}{\partial t} = -\nabla\cdot j_i\tag{5} \end{equation}

With this tool we can continue investigating e.g. the evolution of the electric potential in bentonite diffusion tests, as we begun in the latest blog post in the semi-permeability series.

Test case: CaCl2 through-diffusion

As we noted previously, compacted bentonite often only show minor effects of electromigration, with diffusion potentials below 1 mV. This is because substantial counter-ion concentrations reduce the size of the electric field required to counteract effects from differences in ion mobilities (the denominator in eq. 4 is always large). To see some effect when exploring the full solution of the Nernst-Planck equation, let’s therefore focus on one of the hypothetical, “extreme” cases considered in the previous post. This test case considers diffusion of CaCl2 through bentonite of quite low density (\(c_\mathrm{IL}\) = 1.5 M), with a 10 times lower calcium mobility relative to chloride, and has a 0.25 M CaCl2 reservoir on the left side and a 0.05 M CaCl2 reservoir on the right side.

In the previous consideration we did not solve the full Nernst-Planck problem (eqs. 1 — 5), but imposed a linear concentration gradient across the clay and calculated the corresponding electric potential (eq. 4) and contributions to the flux (eqs. 2 and 3). For the “extreme” CaCl2-case, the results of this approach looked like this (relevant model parameters are summarized below the figures)

Model parameters: \(c_\mathrm{IL}\)= 1.5 M, \(D_\mathrm{Ca}^\mathrm{clay}\) = 0.203⋅10-10 m2/s, \(D_\mathrm{Cl}^\mathrm{clay}\) = 2.03⋅10-10 m2/s, c(1) = 0.25 M CaCl2,c(2) = 0.05 M CaCl2.

The concentration gradient contributions to the fluxes are constant throughout the clay domain (dashed lines), signifying the imposed linear concentration profile. In contrast, the total fluxes varies significantly (full lines), indicating that the system is not in steady-state. As the calcium mobility is 10 times lower than the chloride mobility, it is no surprise that the calcium flux is dominated by the electromigration contribution (dotted lines). Also the chloride flux is seen to be significantly influenced (retarded) by the electric field. The diffusion potential was evaluated to -12.1 mV, while the Donnan contribution to the membrane potential is +18.8 mV, giving an evaluated membrane potential of +6.7 mV (for a fuller discussion on these potentials, see the previous post). But note that this “membrane potential” is not completely accurately evaluated, as the linear concentration profile does not correspond to steady-state.

By solving the full Nernst-Planck framework (eqs. 1 — 5), we can investigate the size of the error in the above treatment. The full solution looks like this

Model parameters: \(c_\mathrm{IL}\) = 1.5 M, \(D_\mathrm{Ca}^\mathrm{clay}\) = 0.203⋅10-10 m2/s, \(D_\mathrm{Cl}^\mathrm{clay}\) = 2.03⋅10-10 m2/s, c(1) = 0.25 M CaCl2,c(2) = 0.05 M CaCl2.

This figure shows snapshots of profiles of the chloride concentration, the electric potential and the chloride flux at 0.25, 0.5, 1, 2 and 8 days; the initial condition is a chloride free calcium clay. In the flux diagram are plotted both the full flux (solid lines) and the concentration gradient contribution3 (dashed lines). Although only barely noticeable, the steady-state has a slightly non-linear concentration profile, resulting in a varying “Fickian” flux contribution (\(j_\mathrm{conc}\)). However, compared with the case with imposed linear concentration profiles, the diffusion potential changes only slightly, to -12.2 mV, giving a true steady-state membrane potential of +6.6 mV.

Below are compared the steady-state fluxes for both calcium and chloride with the fluxes evaluated from imposing linear concentration profiles.

The — now constant — total Cl flux is 3.64⋅10-6 mol/m2/s, to be compared with a variation between 3.0⋅10-6 mol/m2/s and 4.5⋅10-6 mol/m2/s in the previous calculation. Likewise, the previously constant chloride concentration gradient flux of 4.88⋅10-6 mol/m2/s now varies between approximately 6⋅10-6 mol/m2/s and 4⋅10-6 mol/m2/s.

We should remind ourselves that we have considered an “extreme” case, and also take a look at a more realistic/relevant model, to illustrate the often quite insignificant contribution of electromigration in bentonite. In the plot below we are solving the full Nernst-Planck framework for CaCl2 through diffusion, as previously, but now with a denser clay (\(c_\mathrm{IL}\) = 3.0 M) and by setting the Cl/Ca diffusivity ratio equal to the corresponding value in bulk water (2.57).

Model parameters: \(c_\mathrm{IL}\) = 3.0 M, \(D_\mathrm{Ca}^\mathrm{clay}\) = 0.79⋅10-10 m2/s, \(D_\mathrm{Cl}^\mathrm{clay}\) = 2.03⋅10-10 m2/s, c(1) = 0.25 M CaCl2,c(2) = 0.05 M CaCl2.

We considered also this case in the previous treatment. The diffusion potential is here -1.1 mV, and the membrane potential +18.8 mV, which we also evaluated previously. Also the steady-state fluxes are essentially the same as previously evaluated: the chloride steady-state flux is 3.6⋅10-6 mol/m2/s, while in the case with an imposed linear concentration profile, the total flux varies between 3.5⋅10-6 mol/m2/s and 3.6⋅10-6 mol/m2/s. All of these similarities are of course a manifestation of that the linear concentration profiles are very close to being the true steady-state profiles.

Note how similar the concentration gradient contribution is to the total flux for chloride during the entire process, in contrast to the behavior in the “extreme” case. This indicates that chloride here essentially show “Fickian” diffusion, at the rate set by its individual diffusion coefficient. It should be remembered, however, that the full transport process is still influenced by electromigration, since charge neutrality is required (the slower calcium ions are “boosted” by the electric field). But such effects are in principle no different from how salts generally diffuse in ordinary solutions.

The examples in Tournassat and Steefel (2015)

Let’s also take a new look at “example 2” and “example 3” in Tournassat and Steefel (2015), which we discussed quite thoroughly in the review of this paper.

In the original publication, these examples consider a clay component consisting of parallel bulk water and “diffuse layer” domains, and linear concentration profiles are imposed in the bulk water; a main NaCl concentration decreases from 0.1 M to 0.001 M in both examples, while different 22Na and 36Cl concentrations are maintained. In “example 2”, the tracer concentrations decrease in proportion with the main bulk electrolyte, from 10-9 M to 10-11 M. In “example 3” the bulk water tracer concentrations are kept at a constant value 10-9 M.

Here we instead consider only the “diffuse layer” domain, and impose the bulk water concentrations as boundary conditions. The “diffuse layer” domain is characterized by the parameter \(c_\mathrm{IL}\) = 0.33 M (Tournassat and Steefel (2015) call this quantity \(q\)) and is assumed initially free of chloride and tracers. Solving the full Nernst-Planck equations gives the following evolution of the electric potential (this is the same for both examples)

Model parameters: \(c_\mathrm{IL}\)= 0.33 M, \(D_\mathrm{Na}^\mathrm{clay}\) = 1.33⋅10-10 m2/s, \(D_\mathrm{Cl}^\mathrm{clay}\) = 2.03⋅10-10 m2/s, c(1) = 0.1 M NaCl, c(2) = 0.001 M NaCl.

The diffusion potential is -1.0 mV, and is comparable to the potential difference we evaluated when reviewing Tournassat and Steefel (2015). The membrane potential is +115.2 mV.

For “example 2”, the concentration and flux evolution for the tracers looks like this

Model parameters: \(c_\mathrm{IL}\)= 0.33 M, \(D_\mathrm{Na}^\mathrm{clay}\) = 1.33⋅10-10 m2/s, \(D_\mathrm{Cl}^\mathrm{clay}\) = 2.03⋅10-10 m2/s, c(1) = 0.1 M NaCl, c(2) = 0.001 M NaCl. c(1)tracers = 10-9 M, c(2)tracers = 10-11 M

In this example, the steady-state flux of both the chloride and the sodium tracer mimics the flux of the main electrolyte; the steady-state flux is the same for both species (5.6⋅10-15 mol/m2/s) and a factor 108 smaller than for the main electrolyte (5.6⋅10-7 mol/m2/s). As the clay incorporates much more 22Na than 36Cl in the transient stage, the initial fluxes are considerably larger for 22Na. The 36Cl flux is in essence completely “Fickian”, i.e. differences between total fluxes and the concentration gradient contributions are negligible. The steady-state 22Na flux, in contrast, has a non-negligible contribution from electromigration, in the same proportion as the stable sodium (this is not resolved in the figure above, but looks qualitatively the same as for the stable sodium, displayed below). These results are completely in line with what we earlier have concluded for the flux of these tracers.

For “example 3”, the tracer concentration and flux evolution looks like this

Model parameters: \(c_\mathrm{IL}\) = 0.33 M, \(D_\mathrm{Na}^\mathrm{clay}\) = 1.33⋅10-10 m2/s, \(D_\mathrm{Cl}^\mathrm{clay}\) = 2.03⋅10-10 m2/s, c(1) = 0.1 M NaCl, c(2) = 0.001 M NaCl. c(1)tracers = 10-9 M, c(2)tracers = 10-9 M

The evolution for 36Cl is basically the same as for “example 2”. This is expected, as the right hand side boundary conditions set the clay concentration at this interface to effectively zero in both cases. The behavior of 22Na, on the other hand, is very different as compared with “example 2”. Now a huge 22Na concentration (relatively speaking) is established at the right hand side boundary, which results in a flux in the opposite direction with much larger magnitude (the steady-state flux is -4.3⋅10-12 mol/m2/s). This illustrates the mechanism behind the seeming uphill diffusion effect. Note that the 22Na flux in this case is essentially completely “Fickian” during the entire simulation.

These simulations thus confirm what we concluded in the review: all “diffuse layer” tracer fluxes in the examples in Tournassat and Steefel (2015) are essentially “Fickian”, regardless of what is stated and concluded in that publication. Likewise, the “uphill” diffusion effect does NOT rely on any significant flux contributions from electromigration.

The only significant effect of electromigration in this model is, as we have concluded previously, for the transport of the main electrolyte sodium (or the corresponding 22Na transport in “example 2”, which mimics the main sodium flux). The concentration and flux evolution for main chloride and sodium looks like this (this is the same in both examples).

Model parameters: \(c_\mathrm{IL}\) = 0.33 M, \(D_\mathrm{Na}^\mathrm{clay}\) = 1.33⋅10-10 m2/s, \(D_\mathrm{Cl}^\mathrm{clay}\) = 2.03⋅10-10 m2/s, c(1) = 0.1 M NaCl, c(2) = 0.001 M NaCl.

The total flux are the same for chloride and sodium at all times, as a consequence of the requirement of charge neutrality.4 To achieve this, the sodium flux is “boosted” by the induced electric field (electromigration). The effect is mainly noticeable on the cation, because it has a much larger concentration.

Footnotes

[1] \(\tilde{\psi}\) is normalized by the thermal voltage as \(\tilde{\psi} = \psi/V_T\), where \(V_T = RT/F\). Here we adopt \(V_T\) = 25.7 mV.

[2] The corresponding description for Fickian diffusion in the homogeneous mixture model is found here.

[3] The concentration gradient contributions to the fluxes (\(j_\mathrm{conc}\)) in the diagrams presented here are not color coded (but black), in order for them to be visible also when \(j \approx j_\mathrm{conc}\). I hope it is still reasonably easy to understand which flux profile corresponds to which time.

[4] Strictly, the main chloride and sodium fluxes differ ever so slightly in “example 3”, because there is also an 22Na/stable Na exchange process ongoing. But this difference is not noticeable for the main ions, as the reservoir tracer concentrations are a factor 106 — 108 smaller.