*Article* **Time Evolution Study of the Electric Field Distribution and Charge Density Due to Ion Movement in Salty Water**

**Vasileios Bartzis <sup>1</sup> and Ioannis E. Sarris 2,\***


**Abstract:** Desalination and water purification through the ion drift of salted water flow due to an electric field in a duct is perhaps a feasible membrane-free technology. Here, the unsteady modulation of ion drift is treated by employing the Poison–Nernst–Plank (PNP) equations in the linear regime. Based on the solution of the PNP equations, the closed-form relationships of the charge density, the ion concentration, the electric field distribution and its potential are obtained as a function of position and time. It is found that the duration of the ion drift is of the order of one second or less. Moreover, the credibility of various electrical circuit models is examined and successfully compared with our solution. Then, the closed form of the surface charge density and the potential that are calculated without the linear approximation showed that the compact layer is crucial for the ion confinement near the duct walls. To test this, nonlinear solutions of the PNP equations are obtained, and the limits of accuracy of the linear theory is discussed. Our results indicate that the linear approximation gives accurate results only at the fluid's bulk but not inside the double layer. Finally, the important issue of electric field diminishing at the fluid's bulk is discussed, and a potential method to overcome this is proposed.

**Keywords:** electric field; salt ion drift; water duct flow; diffuse layer thickness

## **1. Introduction**

Research on desalination and ion removal methods from a water solution, in general, is of crucial importance for freshwater sustainability. Except for the classical methods, such as the multi-stage flash distillation (MSF), multiple effect distillation and vapor compression distillation [1–11], other methods based on special membranes or electrodes, such as the electrodialysis and the reverse osmosis (RO) electro-deionization, electrophoresis, electroosmosis and capacitive deionization [12–36], have recently been developed. The bottleneck of all the above methods, which are based on the ion drift due to electric fields, is the attempt to resolve the formation of a double layer, i.e., a boundary layer adjacent to the electrodes. Upon its formation, due to ion confinement at the higher potential region of the electrodes, no electric force remains to attract other ions from the fluid's bulk.

Recently, we proposed the electric ion drift (EID) method, which works similarly to the capacitive deionization method but without the use of membranes or special electrodes [37]. Our method is based on the ion drift due to the application of an electric field from a capacitor externally placed from the solution. A similar configuration to the EID method is studied by molecular dynamics simulations and proved that the proposed method can successively drift salt ions in nanochannels [27]. This method works similarly to the capacitive deionization method [38–45], but in its original concept, it is more elegant since no porous electrodes are required. The EID method is of the same energy consumption as the capacitive deionization method of about 0.2 kWh/m<sup>3</sup> and can only be efficient in low salt concentrations as it is presented in our following work [46] for solution equilibrium but without the use of membranes or special electrodes. There, the spatial distributions of

**Citation:** Bartzis, V.; Sarris, I.E. Time Evolution Study of the Electric Field Distribution and Charge Density Due to Ion Movement in Salty Water. *Water* **2021**, *13*, 2185. https:// doi.org/10.3390/w13162185

Academic Editors: Cristina Palet and Julio Bastos-Arrieta

Received: 15 July 2021 Accepted: 8 August 2021 Published: 10 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the electric field and potential, the ion concentration and also the electric charge confined adjacent to the walls are presented for various applied external fields and duct widths.

An in-depth investigation is performed here for the salty water solution flow in a duct that is submerged into an external electric field using the Poisson–Nernst–Planck (PNP) equations and considering the double-layer formation bottleneck. Both the temporal and the spatial distributions are studied. The model equations and details are described in Section 2. By considering the Stern model, the closed-form solutions of charge density, ion concentration, electric field intensity and potential are presented in the linear regime in Section 3. The spatial and temporal distributions found are compared to existing electric circuit models and the equilibrium distributions of Reference [46].

The linear approximation is incapable of providing reasonable results neither for the surface charge density nor for the electric field or potential in the electrodes nearby area (adjacent to the duct walls), i.e., inside the double layer. This analysis is made in Section 4, and a novel method to calculate the surface charge density and the other field quantities is proposed. Lastly, the non-linear PNP equation is deployed in Section 5, and the conditions under which the linear model is valid are investigated. Finally, in Sections 6 and 7, we proceed to compare our results with any experimental data available and to draw the most important conclusions.

#### **2. Model Equations and Mathematical Evaluation**

It is assumed that positive and negative ions of equal concentrations are diluted in a water solution. The salted water flows in a duct, while the duct is subjected to a steadystate and homogeneous electric field created by a pair of plate electrodes, which are under voltage *V* (Figure 1). The external electric field intensity *E* has a direction normal to the plate electrodes from positive to negative along the *y*-axis. The fluid is flowing in the duct perpendicular to the external electric field with velocity <sup>→</sup> *υ* . Thus, ions are moving due to the combined motion of the water flow along the streamwise direction and along the *y*-axis due to the electric force with velocity <sup>→</sup> *υ <sup>y</sup>*. The *y*-axis ionic flux of the positive ions is given by the relation *J*<sup>+</sup> = *C*+*υy*, where *C*<sup>+</sup> is the concentration of the positive ions. From the analysis that is presented in Appendix A and Equation (A5), we have

$$J\_{+} = -\mathbb{C}\_{+} z e \mu\_{+} \frac{\partial \varrho}{\partial y} - D\_{+} \frac{\partial \mathbb{C}\_{+}}{\partial y} \tag{1}$$

where *<sup>z</sup>* is the number of overflow protons; *<sup>e</sup>* <sup>=</sup> 1.6×10−<sup>19</sup> Cb; *<sup>µ</sup>*<sup>+</sup> <sup>≡</sup> <sup>1</sup> 6*πνr* (see Appendix A) is the mobility of the positive ions, with *ν* being the dynamic viscosity of water and *r* being the effective radius of the ions; *ϕ* is the electric potential in the fluid's bulk; and *D*+ = *µ*+*kT* is the diffusion coefficient of the positive ion. *<sup>k</sup>* <sup>≡</sup> *<sup>R</sup> N<sup>A</sup>* is the Boltzmann constant; *R* = 8.314 J/(mol K); *N<sup>A</sup>* is the Avogadro constant; and *T* is the absolute temperature, which is considered *T* = 300 K throughout this work.

In the above Equation (1), the electric field is substituted by *E* = − *∂ϕ ∂y* . Following the same procedure, the ionic flux of the negative ions is read as *J*<sup>−</sup> = +*C*−*zeµ*<sup>−</sup> *∂ϕ <sup>∂</sup><sup>y</sup>* − *D*<sup>−</sup> *∂C*<sup>−</sup> *∂y* , where *z* is the number of overflow electrons. Considering that both positive and negative ions have the same *z*, the same mobility *µ*<sup>+</sup> = *µ*<sup>−</sup> = *µ* and, thus, *D*<sup>+</sup> = *D*<sup>−</sup> = *D*, the ionic fluxes can be read more simply as

$$J\_{\pm} = \mp \mathbb{C}\_{\pm} z e \mu \frac{\partial \varphi}{\partial y} - D \frac{\partial \mathbb{C}\_{\pm}}{\partial y} \tag{2}$$

where the sign of the ionic fluxes indicates the direction of motion and *z* is considered positive.

The concentration conservation equation may read as

$$\frac{\partial \mathbb{C}\_{\pm}}{\partial t} = -\frac{\partial}{\partial y} (\mathsf{J}\_{\pm}) = -\frac{\partial}{\partial y} (\mp \mathsf{C}\_{\pm} z e \mu \frac{\partial \rho}{\partial y} - D \frac{\partial \mathsf{C}\_{\pm}}{\partial y}) \tag{3}$$

where the concentration *<sup>C</sup>*<sup>±</sup> is in (mol/*m*<sup>3</sup> ) and the ionic fluxes *J* in mol/ *m*2 ·*s* . The charge density due to the above ion concentration reads as state and homogeneous electric field created by a pair of plate electrodes, which are under voltage *V* (Figure 1). The external electric field intensity *E* has a direction normal to the plate electrodes from positive to negative along the *y*‐axis. The fluid is flowing in the duct

electric circuit models and the equilibrium distributions of Reference [46].

$$
\rho = \mathbf{C}\_{+} zF - \mathbf{C}\_{-} zF = zF(\mathbf{C}\_{+} - \mathbf{C}\_{-}) \tag{4}
$$

and *F* = 96,485.34 C/mol is the Faraday constant. axis due to the electric force with velocity ⃗௬. The *y*‐axis ionic flux of the positive ions is

**2. Model Equations and Mathematical Evaluation**

draw the most important conclusions.

*Water* **2021**, *13*, x FOR PEER REVIEW 2 of 20

duct widths.

efficient in low salt concentrations as it is presented in our following work [46] for solution equilibrium but without the use of membranes or special electrodes. There, the spatial distributions of the electric field and potential, the ion concentration and also the electric charge confined adjacent to the walls are presented for various applied external fields and

An in‐depth investigation is performed here for the salty water solution flow in a duct that is submerged into an external electric field using the Poisson–Nernst–Planck (PNP) equations and considering the double‐layer formation bottleneck. Both the temporal and the spatial distributions are studied. The model equations and details are described in Section 2. By considering the Stern model, the closed‐form solutions of charge density, ion concentration, electric field intensity and potential are presented in the linear regime in Section 3. The spatial and temporal distributions found are compared to existing

The linear approximation is incapable of providing reasonable results neither for the surface charge density nor for the electric field or potential in the electrodes nearby area (adjacent to the duct walls), i.e., inside the double layer. This analysis is made in Section 4, and a novel method to calculate the surface charge density and the other field quantities is proposed. Lastly, the non‐linear PNP equation is deployed in Section 5, and the conditions under which the linear model is valid are investigated. Finally, in Sections 6 and 7, we proceed to compare our results with any experimental data available and to

It is assumed that positive and negative ions of equal concentrations are diluted in a water solution. The salted water flows in a duct, while the duct is subjected to a steady‐

By applying the 1D Poisson equation in the *y*-direction of the duct, it is found that given by the relation ା ൌ ା௬, where ା is the concentration of the positive ions. From

the combined motion of the water flow along the streamwise direction and along the *y*‐

$$\frac{\partial^2 \varphi}{\partial y^2} = -\frac{\rho}{\varepsilon} \rightarrow \frac{\partial^2 \varphi}{\partial y^2} = -\frac{zF(\mathbb{C}\_+ - \mathbb{C}\_-)}{\varepsilon} \tag{5}$$

The above Equations (3) and (5) constitute the so-called Poisson–Nernst–Planck (PNP) equations system, which is solved in the next Section. where *<sup>z</sup>* is the number of overflow protons; ൌ 1.6 ൈ 10ିଵଽ Cb; ା <sup>≡</sup> <sup>ଵ</sup> గఔ (see Appendix

It should be noticed that the present study is in accordance with the Debye–Huckel theory [47–50], also used in References [37,46], where water is considered as a continuous dielectric medium. Thus, the external electric field effect in water is to initiate change in its electric permittivity in the solution to *ε* = *εrε*0, where *ε<sup>r</sup>* ≈ 80 is the relative permittivity of the water and *<sup>ε</sup>*<sup>0</sup> <sup>=</sup> 8.85 <sup>×</sup> <sup>10</sup>−<sup>12</sup> *<sup>F</sup>*/*m*. A) is the mobility of the positive ions, with being the dynamic viscosity of water and being the effective radius of the ions; *φ* is the electric potential in the fluid's bulk; and ା ൌ ା is the diffusion coefficient of the positive ion. ≡ ோ ேಲ is the Boltzmann constant; *R* = 8.314 J/(mol K); ௮ is the Avogadro constant; and *T* is the absolute temperature, which is considered *T* = 300 K throughout this work*.*

**Figure Figure 1. 1.** Flow Flow configuration ( configuration (**aa**) and indicative salt ion distribution (**b**) of the present model. ) and indicative salt ion distribution (**b**) of the present model.

#### **3. PNP Equations Solution for the Stern Model**

#### *3.1. Boundary Conditions*

For the study of the model, the general Stern model can be used where the ions have finite dimensions and thus can only approach to a small distance from the wall of the duct. The contribution of the Stern model is analytically discussed in Appendix B.

Before discussing the solution of the PNP equations, we must define the appropriate boundary conditions. In the center of the duct, at *y* = *<sup>L</sup>* 2 , the electric potential is kept equal to zero, *ϕ y* = *<sup>L</sup>* 2 = 0, while *ϕ*(0) and −*ϕ*(0) are the potential at *y* = 0 and *y* = *L*, respectively. Thus, by considering a linear relation of distance for the potential at the compact part of the double layer adjacent to the two electrodes, we have

$$
\varphi = \pm \varphi(0) \pm \lambda\_s \frac{\partial \varphi}{\partial y} \quad \text{for} \ y = 0, L \tag{6}
$$

where *λ<sup>s</sup>* is the effective thickness of the compact part of the double layer [51,52] (see Appendix B).

Moreover, since the process is a non-Faradaic one [52,53], because no charge is transferred through the electrodes, the ion fluxes and the current density should be zero at the boundaries, i.e.,

$$J\_{\pm} = \mp \mathcal{C}\_{\pm} z e \mu \frac{\partial \varphi}{\partial y} - D \frac{\partial \mathcal{C}\_{\pm}}{\partial y} = 0 \text{ for } y = 0, L \tag{7}$$

$$\text{si } = \text{zeN}\_A(f\_+ - f\_-) = 0 \text{ for } y = 0, L \tag{8}$$

Finally, the initial conditions of the model are such that a uniform ion distribution is applied at *t* < 0, together with a zero potential.

#### *3.2. Linearized Solution of the PNP Equation*

We considered that the ion concentration in the fluid's bulk is supposed to slightly change linearly along the *y*-direction from 0 to *L*. As the ion concentrations are linearly increased or decreased to opposite walls, the total concentration is constant; i.e., the total charge density is lower than the density of each ion (either positive or negative) [52,54–56]. Thus, most of the solution bulk is quasi-neutral, except in the double layer, i.e., a boundary layer of width of the order of 1 µm or less. Thus, from Equation (3), it is found that

$$\frac{\partial \mathcal{C}\_{+}}{\partial t} = +D \frac{\partial^2 \mathcal{C}\_{+}}{\partial y^2} + ze\mu \mathcal{C}\_{+} \frac{\partial^2 \varrho}{\partial y^2} \tag{9}$$

$$\frac{\partial \mathcal{C}\_{-}}{\partial t} = +D \frac{\partial^2 \mathcal{C}\_{-}}{\partial y^2} - ze\mu \mathcal{C}\_{-} \frac{\partial^2 \mathcal{q}}{\partial y^2} \tag{10}$$

Considering that *C*<sup>+</sup> + *C*<sup>−</sup> = 2*CM*, where *C<sup>M</sup>* is the concentration of each ion in the center of the duct (equal for the two ions due to symmetry), subtracting the above equations and by recalling that the charge density is given by the relationship:

$$\rho = zF(\mathbb{C}\_{+} - \mathbb{C}\_{-}) \tag{11}$$

it is found that

$$\frac{1}{D}\frac{\partial\rho}{\partial t} = \frac{\partial^2\rho}{\partial y^2} - \kappa^2\rho\tag{12}$$

where

$$\kappa = \sqrt{\frac{2z^2 \mathcal{C}\_M F^2}{\varepsilon RT}}\tag{13}$$

Thus, the width *κ* −1 corresponds to the so-called Debye length due to the capacitor that may simulate the diffuse layer according to the Gouy–Chapman theory at the edges of the duct. The analytical solution of Equation (12) is given in Appendix C and results in the following relations:

$$\rho = a \sinh\left(\kappa \left(y - \frac{L}{2}\right)\right) \left(1 - e^{-\frac{L}{\tau}}\right) \tag{14}$$

$$\frac{\partial \varphi}{\partial y} = -\frac{\kappa}{\varepsilon \kappa} \cosh \left( \kappa \left( y - \frac{L}{2} \right) \right) \left( 1 - e^{-\frac{t}{\tau}} \right) - \alpha \frac{1}{\tau} e^{-\frac{t}{\tau}} \frac{\cosh \left( \kappa \frac{L}{2} \right)}{\varepsilon D \kappa^3} \tag{15}$$

$$\varphi = -a \frac{\sinh\left[\kappa \left(y - \frac{L}{2}\right)\right]}{\varepsilon \kappa^2} \left(1 - e^{-\frac{L}{\tau}}\right) - a \frac{1}{\tau} e^{-\frac{L}{\tau}} \frac{\left(y - \frac{L}{2}\right)}{\varepsilon D \kappa^3} \cosh\left(\kappa \frac{L}{2}\right) \tag{16}$$

where

$$\alpha = \frac{\wp(0)\varepsilon\kappa^2}{\cosh\left(\kappa\frac{L}{2}\right)\left(\lambda\_s\kappa + \tanh\left(\kappa\frac{L}{2}\right)\right)}\tag{17}$$

$$\tau = \frac{L\left(1 + \frac{2\lambda\_s}{L}\right)}{2D\kappa\left(\lambda\_s \kappa + \tanh\left(\kappa \frac{L}{2}\right)\right)}\tag{18}$$

and

$$
\kappa^2 = \frac{2z^2 \mathcal{C}\_M F^2}{\varepsilon RT} \tag{19}
$$

Thus, from the charge density, it is found that the concentrations of the positive and negative ions are varied as

$$\mathcal{C}\_{+} = \mathcal{C}\_{M} + \frac{a \sinh\left(\kappa \left(y - \frac{L}{2}\right)\right)}{2zF} \left(1 - e^{-\frac{L}{\tau}}\right) \tag{20}$$

$$\mathcal{C}\_{-}=\mathcal{C}\_{M}-\frac{a\sinh\left(\kappa\left(y-\frac{L}{2}\right)\right)}{2zF}\left(1-e^{-\frac{L}{\tau}}\right)\tag{21}$$

Validation of the Electric Circuit Models for Long Ducts

In case of a duct of width of mm order or wider, *<sup>L</sup>* <sup>≥</sup> <sup>10</sup>−<sup>3</sup> m, and since *<sup>λ</sup><sup>s</sup>* <sup>∼</sup> <sup>1</sup> <sup>−</sup> <sup>10</sup> . *A*, it is found for the rate *<sup>λ</sup><sup>s</sup> <sup>L</sup>* <sup>≤</sup> <sup>10</sup>−<sup>6</sup> and, thus, may be considered almost equal to zero. Moreover, for such length of the duct and since *<sup>κ</sup>* <sup>∼</sup> <sup>10</sup><sup>8</sup> *<sup>m</sup>*−<sup>1</sup> , it can also be considered that tanh *κ L* 2 ≈ 1. Thus, after these simplifications:

$$\pi = \frac{L}{2D\kappa(\lambda\_s \kappa + 1)}\tag{22}$$

$$\alpha = \frac{\varphi(0)\varepsilon\kappa^2}{\cosh\left(\kappa\frac{L}{2}\right)(\lambda\_s\kappa + 1)}\tag{23}$$

As it may be observed, the quantities in Equations (14)–(21) are exponentially varied, and *τ* is a characteristic time, the so-called time constant. Here, the above theory fits to the Gouy–Chapman and Stern theories that consider the solution as an RC circuit (see Appendix B), because the characteristic time of the above analysis is exactly the same as the time of the RC circuit (compare Equations (22) and (A10)). Since the fluid bulk is in quasi-*L*(*κ* −1 ) 2

neutral condition, its resistance per unit surface is equal to *R* = *Dε* (Equation (A7)), and for the time constant effective capacitor *τ* = *RCtot* (due to the charge confinement at the edges), we have

$$\mathcal{C}^{tot} = \frac{\varepsilon}{2(\lambda\_S + \kappa^{-1})} \tag{24}$$

As it may be observed, Equation (24) is exactly the same as the Stern's model (see Appendix B) as would be expected. Thus, we have shown that the electric model in long ducts is equal to the linear PNP equation solution, and it is valid under the same assumptions for the solution. The time duration to the end of the phenomenon is plotted in Figure <sup>2</sup> for the case of Na ions, for the salted water solution of NaCl (*D*Na <sup>=</sup> 1.3 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>m</sup><sup>2</sup> s ), as a function of the initial concentration *C<sup>M</sup>* and various duct widths. As always for the exponentially decay variations, equilibrium is expected after 5*τ*. Thus, here it will also be achieved in less than a second at the most as it is observed in Figure 2.

#### *3.3. Long Time Behavior*

As the time increases, then (1 − *e* − *t <sup>τ</sup>* ) → 1 and *e* − *t <sup>τ</sup>* → 0, and Equations (14)–(21) have their equilibrium forms as

$$\rho = \kappa \sinh\left(\kappa \left(y - \frac{L}{2}\right)\right) \tag{25}$$

$$\frac{\partial \varphi}{\partial y} = -\frac{a}{\varepsilon \kappa} \cosh \left( \kappa \left( y - \frac{L}{2} \right) \right) \tag{26}$$

$$\varphi = -\alpha \frac{\sinh\left[\kappa \left(y - \frac{L}{2}\right)\right]}{\varepsilon \kappa^2} \tag{27}$$

$$\mathcal{C}\_{+} = \mathcal{C}\_{M} + \frac{a \sinh\left(\kappa \left(y - \frac{l}{2}\right)\right)}{2zF} \tag{28}$$

$$\mathcal{C}\_{-}=\mathcal{C}\_{M}-\frac{\alpha \sinh\left(\kappa \left(y-\frac{L}{2}\right)\right)}{2zF} \tag{29}$$

The equilibrium solutions of Equations (25)–(29) can be easily verified by the recent results of Reference [46] for the present water solution with *<sup>κ</sup>* <sup>=</sup> 1.027 <sup>×</sup> <sup>10</sup><sup>8</sup> *zC* 1 2 *<sup>M</sup>* (in S.I. (System International)) (from Equation (19)) and *zNa* = *zCl* = 1. For *L* = 0.0015, 0.015, 0.15 m, where tanh *κ L* 2 <sup>≈</sup> <sup>1</sup> and thus sinh *κ L* 2 <sup>≈</sup> cosh *κ L* 2 and for *<sup>C</sup><sup>M</sup>* <sup>=</sup> 8.6, 100, 400 mol/m<sup>3</sup> , and by considering the effect as negligible due to the compact layer, i.e., *λSκ* → 0, Equations (25)–(27) are exactly the same as in Reference [46]. However, the concentration Equations (28) and (29) are different than the ones of Reference [46] that read as

$$\mathcal{C}\_{+} = \mathcal{C}\_{M} e^{-\frac{F\rho(0)}{RT} \frac{\sinh\left[\kappa\left(\frac{L}{2} - y\right)\right]}{\sinh\left(\kappa\frac{L}{2}\right)}}\tag{30}$$

$$\mathcal{C}\_{-}=\mathcal{C}\_{M}e^{+\frac{F\_{\Psi}(0)}{RT}\frac{\sinh\left[\kappa\left(\frac{L}{2}-y\right)\right]}{\sinh\left(\kappa\frac{L}{2}\right)}}\tag{31}$$

The first observation is that Equations (28) and (29) are expected to differ from Equations (30) and (31) of Reference [46] due to the linearization of the exponential term (*e <sup>x</sup>* = 1 + *x* + *<sup>x</sup>* 2 2! + *<sup>x</sup>* 3 3! + · · ·). This happens when *C*<sup>+</sup> + *C*<sup>−</sup> = 2*CM*, which is the same as the linearization of the exponential function when only the first two terms of the Taylor expansion are kept. Thus, Equations (30) and (31) are more accurate and are used for comparison in the next sections where nonlinear terms are considered for the PNP equation. Moreover, the linear approximation is valid along the duct width, Reference [46], for

$$\left| \frac{F}{RT} \varphi \right| < 1 \text{ or } \varphi(0) \le 0.026 \text{ V} \tag{32}$$

which is also applied for the linearized PNP solution. By using Equations (25)–(29) and after some calculus, it can easily be determined that the electric potential and electric field become zero at the fluid bulk, as the width of the duct, the initial concentration and the potential *ϕ*(0) are increased as is also discussed in Reference [46]. As it is illustrated below, the linear approximation gives reasonable results at the fluid's bulk.

Furthermore, the surface charge density (in C/m<sup>2</sup> ) can be found by integrating Equation (25) of the linearized PNP solution as

$$
\sigma\_l = \int\_0^y \rho dy = \frac{\varphi(0)\varepsilon\kappa}{(\lambda\_s \kappa + 1)} [\frac{\cosh\left(\kappa \left(y - \frac{L}{2}\right)\right)}{\cosh\left(\kappa \frac{L}{2}\right)} - 1] \omega
$$

and its absolute value adjacent to each electrode is given for *y* = *<sup>L</sup>* 2 as

$$|\sigma\_l| = \frac{\varrho(0)\varepsilon\kappa}{(\lambda\_s\kappa + 1)} (\frac{1}{\cosh\left(\kappa\frac{L}{2}\right)} - 1) \tag{33}$$

Based on the above approximations, <sup>1</sup> cosh(*κ L* 2 ) → 0, the surface charge density is found to be exactly the same as the electric circuit RC model as

$$|\sigma\_l| = \varphi(0)\frac{\varepsilon}{(\lambda\_S + \kappa^{-1})} = \varphi(0)\mathcal{C} \tag{34}$$

**Figure 2.** Time constant variation 5*τ* as a function of *C<sup>M</sup>* for various duct widths.

**Figure 2.** Time constant variation 5 as a function of ெ for various duct widths.

#### *3.3. Long Time Behavior* **4. Exact Evaluation of the Electric Field and the Surface Charge Density**

As the time increases, then ሺ1 െ ି ഓሻ→1 and ି <sup>ഓ</sup> → 0, and Equations (14)–(21) have their equilibrium forms as ρ ൌ sinhቆκ൬y െ <sup>L</sup> The calculation of the electric field and the surface charge density under steady-state conditions is presented below without the linear approximation. Starting from the electric field definition:

$$E = \frac{-\partial\rho}{\partial y} \rightarrow \frac{\partial E}{\partial y} = -\frac{\partial^2 \rho}{\partial y^2} \rightarrow \frac{\partial E}{\partial \rho}\frac{\partial \rho}{\partial y} = \frac{\rho}{\varepsilon} \rightarrow -2E\frac{\partial E}{\partial \rho} = \frac{2\rho}{\varepsilon} \rightarrow d\left(E^2\right) = -\frac{2\rho}{\varepsilon}d\rho\sqrt{\rho}$$

൰ቇ (25)

<sup>ଵ</sup> ଶൗ (in S.I.

ቁ ) and for

φ ൌ െ sinh ቂκ ቀy െ <sup>L</sup> 2ቁቃ εκଶ (27) and by considering *ϕ y* = *<sup>L</sup>* 2 = 0 and *E<sup>M</sup>* the electric field intensity is in the center of the duct, and by integrating it, it is found that

$$\begin{cases} \int\_{E\_M}^{E} d(E^2) = -\int\_0^{\varphi} \frac{2\rho}{\varepsilon} d\varphi \to \mathbb{E}^2 - E\_M^2 = \frac{4}{\varepsilon} \mathbb{C}\_M F \int\_0^{\varphi} \sinh(\frac{\varepsilon F}{RT} \rho) d\rho \to \\\ E = \sqrt{E\_M^2 + \frac{4 \mathcal{C}\_M \mathcal{R} T}{\varepsilon} \left[ \cosh(\frac{\varepsilon F}{RT} \rho) - 1 \right]} \end{cases}$$

The equilibrium solutions of Equations (25)–(29) can be easily verified by the recent

where

$$\rho = \mathbb{C}\_{+}zF - \mathbb{C}\_{-}zF \rightarrow \rho = -2\mathbb{C}\_{M}F\sinh(\frac{zF}{RT}\rho)$$

ଶ ଶ ଶ ெ ൌ 8.6, 100, 400 mol/mଷ, and by considering the effect as negligible due to the compact layer, i.e., ௌ→0, Equations (25)–(27) are exactly the same as in Reference [46]. However, the concentration Equations (28) and (29) are different than the ones of Reference [46] that read as which is valid only inside the diffuse layer and not in the compact layer, thus, in regions where the Boltzmann distribution is valid. In this region, *E<sup>M</sup>* is almost diminished relative to the values of the second part of the undersquare quantity; thus,

$$E \cong \sqrt{\frac{4\mathcal{C}\_M RT}{\varepsilon} [\cosh(\frac{zF}{RT}\varphi) - 1]} = \sqrt{\frac{8\mathcal{C}\_M RT}{\varepsilon}} \left| \sinh\left(\frac{zF}{2RT}\varphi\right) \right|^2$$

and close to the positive electrode:

$$E = \sqrt{\frac{8C\_MRT}{\varepsilon}}\sinh(\frac{zF}{2RT}\rho) = -\frac{\partial\rho}{\partial y}\text{ for }\lambda\_s \le y \le \frac{L}{2} \tag{35}$$

The integration of the 1D Poisson equation from *y* = *λ<sup>S</sup>* to *y* = *<sup>L</sup>* 2 can give the surface charge density near the walls as <sup>R</sup> *<sup>L</sup>* 2 *λS d* 2*ϕ dy*<sup>2</sup> *dy* = − R *L* 2 *λS ρ ε dy* <sup>→</sup> *<sup>d</sup><sup>ϕ</sup> dy L* 2 − *dϕ dy λS* = − R *<sup>L</sup>*/2 *λS ρdy <sup>ε</sup>* = −*σ ε* , where *σ* is the surface charge density (in C/m<sup>2</sup> ) near the positive electrode. Applying that *E<sup>M</sup>* = *dϕ dy L* 2 = 0, it is − *dϕ dy λS* <sup>=</sup> <sup>−</sup>*<sup>σ</sup> <sup>ε</sup>* → *σ* = −*εEλ<sup>S</sup>* , which due to Equation (35) gives

$$
\sigma = -\sqrt{8\varepsilon \mathcal{C}\_M RT} \sinh(\frac{zF}{2RT}\rho\_s) \tag{36}
$$

where *ϕ<sup>s</sup>* is the potential at the outer Helmholtz plane (OHP) at distance *λS*.

It should be noticed that due to continuity, the electric field magnitude is equal at both sides of the OHP, i.e., of the diffuse layer but also of the compact layer side. Moreover, electric field magnitude is almost constant along the compact layer. Thus, for the positive electrode, it is *ϕ<sup>S</sup>* = *ϕ*(0) + *λ<sup>s</sup> dϕ dy λS* = *σλs <sup>ε</sup>* + *ϕ*(0), and by replacing *ϕ<sup>S</sup>* in Equation (36), we have

$$\sigma = -\sqrt{8\varepsilon \mathcal{C}\_M RT} \sinh[\frac{zF}{2RT}(\frac{\sigma \lambda\_s}{\varepsilon} + \varphi(0))] \tag{37}$$

or after substituting

$$\sigma = -3.76 \cdot 10^{-3} \mathcal{C}\_{M}^{\frac{1}{2}} \sinh[z(19.34\rho(0) + 13.65\sigma)] \text{ (S.I.)}\tag{38}$$

The distribution of the surface charge density as a function of *ϕ*(0) is shown in Figure 3 for various concentrations considering the case of salty water, *zNa* = *zCl* = 1. As may be seen by comparing these results with the those from Reference [46], where no compact layer was considered, the deployment of the compact layer results in a reduced charge. Moreover, the charge is practically independent from the solution salt concentration.

As it may be observed, the relation (38) is different than Equation (34) due to linearized PNP equations or the electric models. However, it collapses to this upon linearization, i.e., for sinh[ *zF* <sup>2</sup>*RT σλs <sup>ε</sup>* + *ϕ*(0) i <sup>≈</sup> *zF* <sup>2</sup>*RT σλs <sup>ε</sup>* + *ϕ*(0) , something that can not be applied in general at the double layer. Now, the potential *ϕ<sup>s</sup>* at the outer Helmholtz plane is

$$
\varphi\_S = \varphi(0) - \lambda\_s \sqrt{\frac{8 \mathcal{C}\_M RT}{\varepsilon}} \sinh(\frac{zF}{2RT} \varphi\_s) \tag{39}
$$

or

$$\varphi\_{\rm S} = \varphi(0) - 2.65 \times 10^{-3} \mathcal{C}\_{M}^{\frac{1}{2}} \sinh(19.34 \cdot z \varphi\_{\rm s}) \text{ (S.I)}\tag{40}$$

The potential at the outer Helmholtz plane is shown in Figure 4 as a function of *ϕ*(0) for various concentrations considering the case of salty water, *zNa* = *zCl* = 1. As it can be observed, *ϕ<sup>S</sup>* is scaled almost linear with *ϕ*(0) for values up to 0.1 V. However, as *ϕ*(0) is increased further, *ϕ<sup>S</sup>* is saturated to values lower than 1 V independently of *ϕ*(0).

Now, the potential distribution along *<sup>y</sup>* in the diffuse layer for *<sup>λ</sup><sup>s</sup>* <sup>≤</sup> *<sup>y</sup>* <sup>≤</sup> *<sup>L</sup>* 2 can be found from Equation (35) as

$$\int\_{\varphi\_s}^{\varphi} \frac{d\varphi}{\sinh(\frac{F}{2RT}\varphi)} = -\sqrt{\frac{8\mathbb{C}\_M \mathbb{R} \mathcal{T}}{\varepsilon}} \int\_{\lambda\_s}^{y} dy \to \varphi = \frac{4\mathbb{R}\mathbb{T}}{zF} \tanh^{-1}\{\tanh(\frac{zF\varphi\_s}{4RT}) \cdot e^{-\kappa(y-\lambda\_s)}\}\tag{41}$$

Thus, the electric field intensity inside the compact layer is given by

$$E\_{\lambda\_S} = -\frac{\sigma}{\varepsilon} \to E\_{\lambda\_S} = \sqrt{\frac{8C\_MRT}{\varepsilon}}\sinh[\frac{zF}{2RT}\left(-\lambda\_S E\_{\lambda\_S} + \varphi(0)\right)]\tag{42}$$

or

$$E\_{\lambda\_S} = 5.31 \times 10^6 \mathcal{C}\_M^{\frac{1}{2}} \sinh[z \left(19.34 \varphi(0) - 9.67 \times 10^{-9} \, E\_{\lambda\_S} \right)] \tag{43}$$

As it is observed from Figure 5 for the case of salty water, *zNa* = *zCl* = 1, the electric field intensity at the compact layer takes huge values independently of the ion concentration, especially as the potential increases.

**Figure 3.** Surface charge density dependance on ሺ0ሻ at the outer Helmholtz plane for various ெ. **Figure 3.** Surface charge density dependance on *ϕ*(0) at the outer Helmholtz plane for various *CM*.

2ൗ can be

concentration, especially as the potential increases. **Figure 4.** Potential ௌ dependance on ሺ0ሻ at the outer Helmholtz plane for various ெ. **Figure 4.** Potential *ϕ<sup>S</sup>* dependance on *ϕ*(0) at the outer Helmholtz plane for various *CM*.

**5. Nonlinear PNP Solution Figure 5.** Electric field intensity inside the compact layer dependance on ሺ0ሻ for various ெ. **Figure 5.** Electric field intensity inside the compact layer dependance on *ϕ*(0) for various *CM*.

zଶFଶ RT <sup>E</sup>

பమ ப୷మ zeμ ப ப௬ డశ

∂φ ∂ Cି

∂ଶφ∂yଶ െ zeμ

zଶFଶ RT <sup>E</sup>

∂ሺCା Cିሻ ∂y

∂φ ∂ Cି

డ௬ (44)

(45)

∂ሺCା Cିሻ ∂y

(46)

(46)

(45)

డ௬ (44)

பమ ப୷మ zeμ ப ப௬ డశ

∂ଶφ ∂yଶ െ zeμ

εRT ρ െ

∂ଶCି ∂yଶ െzeμCି

ப୲ ൌ Dபమେశ

Starting from ൌ ሺା െ ିሻ and setting Equation (3) in the form

ப୷మ zeμCା

Starting from ൌ ሺା െ ିሻ and setting Equation (3) in the form

∂yଶ െ ଶFଶሺା ିሻ

εRT ρ െ

∂yଶ െzeμCି

ப୷మ zeμCା

பେశ ப୲ ൌ Dபమେశ

பେశ

∂Cି ∂t ൌ D

∂Cି ∂t ൌ D

1 D ∂ρ ∂t ൌ ∂ଶρ

1 D ∂ρ ∂t ൌ ∂ଶρ

and by subtracting, we have

and by subtracting, we have

**5. Nonlinear PNP Solution**

## **5. Nonlinear PNP Solution**

Starting from *ρ* = *zF*(*C*<sup>+</sup> − *C*−) and setting Equation (3) in the form

$$\frac{\partial \mathcal{C}\_{+}}{\partial t} = +D \frac{\partial^2 \mathcal{C}\_{+}}{\partial y^2} + ze\mu \mathcal{C}\_{+} \frac{\partial^2 \varrho}{\partial y^2} + ze\mu \frac{\partial \varrho}{\partial y} \frac{\partial \mathcal{C}\_{+}}{\partial y} \tag{44}$$

$$\frac{\partial \mathcal{C}\_{-}}{\partial t} = +D \frac{\partial^2 \mathcal{C}\_{-}}{\partial y^2} - ze\mu \mathcal{C}\_{-} \frac{\partial^2 \varrho}{\partial y^2} - ze\mu \frac{\partial \varrho}{\partial y} \frac{\partial \mathcal{C}\_{-}}{\partial y} \tag{45}$$

and by subtracting, we have

$$\frac{1}{D}\frac{\partial\rho}{\partial t} = \frac{\partial^2\rho}{\partial y^2} - \frac{z^2F^2(\mathbb{C}\_+ + \mathbb{C}\_-)}{\epsilon RT}\rho - \frac{z^2F^2}{RT}E\frac{\partial(\mathbb{C}\_+ + \mathbb{C}\_-)}{\partial y} \tag{46}$$

The above Equation (46) can be transformed into the linear Equation (12) subject to *C*<sup>+</sup> + *C*<sup>−</sup> = *constant*. Since along the y-direction the ion concentration is changing, the decrease in the positive ion needs to correspond to the increase in the negative ion for their summary to be constant and equal to 2*CM*, i.e., *C*<sup>+</sup> + *C*<sup>−</sup> = 2*CM*.

This is valid throughout the duct width (for the case of salty water, *zNa* = *zCl* = 1) when *F RT ϕ* <sup>&</sup>lt; <sup>1</sup> or *<sup>ϕ</sup>* <sup>≤</sup> 0.026 V, i.e., in the case of very weak external electric fields as already proved in Reference [46]. However, this is true only in the solution's bulk as the applied potential increases and faults inside the double layer. Since we proved that the time scale of the problem is of the order of a second, we study the final equilibrium state of the ion drift. Similar to Section 3.3, we can assume that the final concentration profiles upon equilibrium that came after the solution of Equations (44)–(46) are of the form

$$\mathcal{C}\_{+} = \mathcal{C}\_{M} e^{+\frac{Fz\rho(0)}{RT(\lambda\_{\rm s}\kappa + \tanh(\kappa\frac{L}{2}))}} \frac{\sinh[\kappa(y - \frac{L}{2})]}{\cosh(\kappa\frac{L}{2})} \tag{47}$$

and

$$\mathcal{C}\_{-}=\mathcal{C}\_{M}e^{-\frac{Fz\rho(0)}{RT(\lambda\_{\rm s}\kappa+\tanh(\kappa\frac{L}{2}))}}\frac{\sinh[\kappa(y-\frac{L}{2})]}{\cosh(\kappa\frac{L}{2})}\tag{48}$$

The effect of the nonlinearity can be investigated when a third term in the Taylor series of the exponential function, *e <sup>x</sup>* = 1 + *x* + *<sup>x</sup>* 2 2! + *<sup>x</sup>* 3 3! + · · · is considered. Then, for the additional term, it is

$$\mathcal{C}\_{+} + \mathcal{C}\_{-} = \mathcal{C}\_{M}(2 + W)\_{\prime} \text{ and } \frac{\partial(\mathcal{C}\_{+} + \mathcal{C}\_{-})}{\partial y} = \mathcal{C}\_{M} \text{Q} \kappa \sinh[2\kappa \left(y - \frac{L}{2}\right)]\_{\prime} \tag{49}$$

where *<sup>W</sup>* = [ *Fzϕ*(0)sinh[*κ*(*y*<sup>−</sup> *<sup>L</sup>* <sup>2</sup> )] *RT*(*λsκ*+tanh(*κ L* <sup>2</sup> )) cosh(*<sup>κ</sup> L* 2 ) ] 2 , and *Q* = [ *Fzϕ*(0) *RT*(*λsκ*+tanh(*κ L* <sup>2</sup> )) cosh(*<sup>κ</sup> L* 2 ) ] 2 .

Thus, from Equation (46) and by setting *∂ρ <sup>∂</sup><sup>t</sup>* = 0 for the steady-state, it is

$$0 = \frac{\partial^2 \rho}{\partial y^2} - \frac{z^2 F^2}{\varepsilon RT} \mathbb{C}\_M (2 + \mathcal{W}) \rho - Y \mathbb{C} \mathbb{C}\_M \text{ where } \mathcal{Y} = \frac{z^2 F^2}{RT} Q \kappa \sinh[2\kappa \left( y - \frac{L}{2} \right)] \tag{50}$$

The shift of the linearity is due to *W* and *Y*. It will be useful to study their profiles along *y* for various potentials *ϕ*(0) (for the case of salty water, *zNa* = *zCl* = 1). In the range of 0.026, 0.2 and 2 V, concentrations are in the range of 8.6, 100 and 400 mol/m<sup>3</sup> with a duct width of *L* = 0.015 m, while the effect of *L* is insignificant.

As it is observed in Figures 6 and 7, only inside the compact layer can a deviation from linearity be encountered, while the difference is very negligible outside this layer as one would expect. Moreover, the width of the deviation is found to be inversely proportional to *C<sup>M</sup>* as it is also discussed in Reference [46].

**Figure 6.** Effect of the nonlinearity on W at *ϕ*(0) = 0.026 (**a**), 0.2 (**b**) and 2 (**c**) along *y* for various *CM*.

ெ.

**Figure 6.** Effect of the nonlinearity on W at ሺ0ሻ ൌ 0.026 ሺሻ, 0.2 ሺሻ and 2 ሺሻ along for various

**Figure 7.** Effect of the nonlinearity on Y along *y* for various *ϕ*(0) at *C<sup>M</sup>* = 8.6 mol/m<sup>3</sup> (**a**), *C<sup>M</sup>* = 100 mol/m<sup>3</sup> (**b**) and *C<sup>M</sup>* = 400 mol/m<sup>3</sup> (**c**).

#### **6. Comparison with Experimental Results**

The only available and relevant experimental data are from the capacitive deionization method and are obtained for the water desalination case. A typical duct width of the order *L* = 1.5 mm is proposed and constructors indicate a performance of about 65% for low salt concentrations. In a similar approach, the authors of Reference [40] succeeded in reducing an initial ion concentration of *Cbe f* = 1.686 gr/L = 29 mole/m<sup>3</sup> to the final concentration of *Ca f ter* = 0.497 gr/L = 8.6 mole/m<sup>3</sup> with *V* = 0.4 V; thus, the electric charge that it is removed per surface area is about <sup>|</sup>*σ*<sup>|</sup> <sup>=</sup> <sup>2967</sup> Cb/m<sup>2</sup> . This large ion drift can only be achieved by the existence of an electric field in the fluid bulk. This is due to the micropores in the capacitive deionization method. However, this is not possible here due to the non-Faradaic conditions to which our analysis refers. Clearly, the phenomenon of ion drift stops long before we have a noticeable reduction in ion in the bulk of the solution. This is due to the zeroing of the electric field in the bulk due to the creation of the double layer on the side walls of the duct. Thus, we need to continuously absorb the compact layer in order to permit additional ions to move close to the duct sidewalls. This can be conducted in an almost continuous way because as we showed, the time to reach the final equilibrium state is of the order of one second or less. This is equivalent to maintaining a non-zero electric field at the fluid's bulk. This mechanism of ion removal from the salted water is discussed in Reference [37] since in this case, the solution had a constant electric field, which is impossible upon creation of the double layer, and in order for the electric field to exist, we have to continuously eliminate this layer. Thus, electric field existence in the duct is the crucial parameter to ensure that an efficient desalination can be achieved by this method. The absorption of the double layer replaces the porous electrodes and if achieved experimentally may be the solution for continuous desalination without the inevitable interruptions that must occur in capacitive deionization to discharge the electrodes.

#### **7. Conclusions**

This study examined the ion concentration, electric field, electric potential profiles and surface charge density inside a salt solution when it is under the effect of an external electric field as a function of time and position. The compact layer is discussed as it forms an additional mechanism to diminish the electric field in the bulk of the fluid. The conditions and the region where the linear approximation gives satisfactory results are examined.

The energy consumption of conventional methods using membranes is in the order of 0.30–2.11 kWh/m<sup>3</sup> , quite costly but provides desalination of a large scale. With a capacitive deionization method, we have a smaller energy consumption of about 0.2 kWh/m<sup>3</sup> . At present with this method, we can manage smaller amounts of water, and it can be used as a secondary desalination method. Research continues intensively on this method mainly with regard to electrodes in order to achieve greater absorbency. With the EID method, we have the same energy consumption, but we are exempt from the cost of electrodes.

The importance of maintaining a non-zero electric field in the fluid's bulk is shown, which can only be succeeded by the continuous absorption of the double layer.

**Author Contributions:** Conceptualization, V.B.; methodology, V.B.; software, I.E.S.; investigation, V.B.; writing—original draft preparation, V.B.; writing—review and editing, I.E.S.; visualization, I.E.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **Appendix A. General Theory of Ion Movement**

It is supposed that a positive ion of concentration C<sup>+</sup> is dissolved in a water solution that is flowing in the streamwise direction of a duct (Figure 1). Moreover, a uniform electric field of magnitude *E* is applied normal to the flow, and the corresponding force acting on the ions is

$$f\_{(1\text{ ion})}^E = \text{zeE} \tag{A1}$$

Then, the result is that ions start to drift. Thus, a spatial concentration distribution along the *y*-direction is formed, which causes the concentration gradient *<sup>∂</sup>C*<sup>+</sup> *∂y* . The pseudoforce per ion due to this is given by the relation:

$$f\_{(1\text{ ion})}^{\nabla \mathbb{C}} = -\frac{RT}{N\_A \mathbb{C}\_+} \frac{\partial \mathbb{C}\_+}{\partial y} \tag{A2}$$

The negative sign indicates that the force direction is opposite to the concentration gradient direction, along the *y*-direction. Moreover, during the ion movement, a friction force also exists between the moving ion and the solvent. This force is opposite to the ion velocity <sup>→</sup> *υ <sup>y</sup>*, which is given approximately by the relation (we consider as usual the ion shape as a small sphere with radius *r* that moves with a small velocity inside a fluid of dynamic viscosity *ν*):

$$f\_{(1\text{ion})}^{fr} = -6\pi\nu r \cdot \upsilon\_y \tag{A3}$$

Quickly, equilibrium is established; thus,

$$
\Sigma f\_{(1ion)}^y = 0 \tag{A4}
$$

By using Equations (A1)–(A3), the velocity *υ<sup>y</sup>* can be found as

$$\upsilon\_{y} = \frac{zeE}{6\pi\upsilon r} - \frac{RT}{\mathbb{C}\_{+}N\_{A}6\pi\upsilon r} \frac{\partial \mathbb{C}\_{+}}{\partial y} \tag{A5}$$

#### **Appendix B. The RC Models—Time Scale**

In Electrochemistry, there is an old theoretical approximation that when an electrolyte is under the effect of an external electric field, it is modeled as an effective electric circuit constituted by one capacitor and one resistance connected inline [52]. In the following, the most important models are reviewed.

## *Appendix B.1. Gouy-Chapman Model*

The very first model used was attributed to Gouy–Chapman [52] in which only the diffuse layer that is formed from the ions' confinement nearby the electrodes' surface is considered. This thin layer can be simulated by a capacitor with a plate separation width: *λ<sup>D</sup>* = *κ* <sup>−</sup><sup>1</sup> = q *εRT* 2*F* 2*z* <sup>2</sup>*C<sup>M</sup>* . Thus, the capacity per unit area is equal to *C<sup>D</sup>* = *<sup>ε</sup> κ*−<sup>1</sup> . Since an additional capacitor is formed at the other wall of the duct, the total capacity per unit area is then

$$\mathbb{C}\_D^{tot} = \frac{\varepsilon}{\mathfrak{Z}(\kappa^{-1})} \tag{A6}$$

Assuming that the bulk of the solution is in an almost uniform initial ion concentration, the electric force due to the external electric field on the ions is fast balanced by the friction force, and their asymptotic velocity can be found from (A4) when *<sup>∂</sup>C*<sup>+</sup> *<sup>∂</sup><sup>y</sup>* = 0 for the positive ions as in *<sup>υ</sup><sup>y</sup>* <sup>=</sup> <sup>−</sup> *ze* 6*πνr ∂ϕ ∂y* . Their ionic flux is also given by *J*<sup>+</sup> = *CMυ<sup>y</sup>* = −*C<sup>M</sup> ze* 6*πνr ∂ϕ ∂y* , while *J*<sup>−</sup> = *C<sup>M</sup> ze* 6*πνr ∂ϕ ∂y* is similarly given for the negative one. The electric current per unit surface is given by *<sup>i</sup>* <sup>=</sup> *zeNA*(*J*<sup>+</sup> <sup>−</sup> *<sup>J</sup>*−) or substituting *<sup>i</sup>* <sup>=</sup> <sup>−</sup> *<sup>ε</sup><sup>D</sup>* (*κ*−<sup>1</sup>) 2 *∂ϕ ∂y* . Comparing this relationship with the Ohm's law *i* = −*σ ∂ϕ ∂y* , the special conductivity is given by the relation *σ* = *<sup>ε</sup><sup>D</sup>* (*k*−<sup>1</sup>) 2 . The resistance per unit surface is given by *R* = *<sup>L</sup> σ* and, thus,

$$R = \frac{L\left(\kappa^{-1}\right)^2}{D\varepsilon} \tag{A7}$$

The time constant of the *RC* circuit determines that characteristic time of the phenomenon as

$$
\pi = \mathcal{R} \mathcal{C}\_D^{tot} = \frac{L\kappa^{-1}}{2D} \tag{A8}
$$

while 5*τ* is its practical duration according to the RC circuit theory. In the above Gouy– Chapman theory, only the diffuse layer and no compact layer is considered.

#### *Appendix B.2. Stern's Model*

Assuming that ions have finite dimensions and, thus, that they can approach in an ionic radius distance to the walls of the duct that is of the order 1 . *A* (further increased by the solvent molecules that surround the ion), the Stern model needs to be used. According to this model, the double layer is formed by two areas:

a. The compact part that is simulated by a Helmholtz capacitor of effective width of the order of molecule that is almost constant. The term effective is due to the possibility that the solvent (i.e., the water here) may not be considered in this region as a continuous media, and the permittivity ε is considered to be similar to the fluid's bulk. Thus, the width of the compact part may be considered the size of the ionic radius, i.e., *λ<sup>s</sup>* ∼ 1 − 10 . *A* [51], while in the present calculations, it is assumed to be equal to *λ<sup>s</sup>* ≈ 5 . *A*. Thus, the capacity of the compact part per unit area is given by *C<sup>H</sup>* = *<sup>ε</sup> λS* .

b. The diffuse layer that is simulated by a Gouy–Chapman-type capacitor with effective width of the diffuse part *λ<sup>D</sup>* = *κ* −1 , and the known capacity per area unit is given by *C<sup>D</sup>* = *<sup>ε</sup> κ*−<sup>1</sup> .

The above two capacitors are both inline connected and connected with the two capacitors existed at the other wall side (near the second electrode). The two electrodes constitute the Stern model [52]. Thus, the total capacity per surface unit is then

$$\mathcal{C}^{tot} = \frac{\varepsilon}{2(\lambda\_S + \kappa^{-1})} \tag{A9}$$

The time constant of the RC circuit of the Stern model is then

$$\pi = \frac{L}{2D\kappa(\lambda\_\text{s}\kappa + 1)}\tag{A10}$$

#### **Appendix C. Solution of the Linear PNP Equation**

Using the Laplace transformation in Equation (12) and assuming *ρ*(*t* = 0) = 0 for the initial current density, we have

$$\frac{\partial^2}{\partial y^2}\overline{\rho} - M^2 \overline{\rho} = 0 \tag{A11}$$

where

$$M^2 = \kappa^2 + \frac{S}{D} \tag{A12}$$

and *ρ* = R <sup>∞</sup> 0 *dte*−*Stρdt* is the |Laplace transformation of *ρ*. The solution of Equation (A11) is of the form

$$\overline{\rho} = B \cosh(My') + A \sinh(My') \tag{A13}$$

where *y* <sup>0</sup> <sup>=</sup> *<sup>y</sup>* <sup>−</sup> *<sup>L</sup>* 2 , and *A*, *B* are constants. Due to antisymmetric conditions around the duct center at *y* = *<sup>L</sup>* 2 , it is *ρ*(*y* 0 ) = −*ρ*(−*y* 0 ), and, thus, *B* = 0, i.e.,

$$\overline{\rho} = A \sinh \left( M \left( y - \frac{L}{2} \right) \right) \tag{A14}$$

where

$$M^2 = \kappa^2 + \frac{S}{D} \tag{A15}$$

In order to determine the *A* constant, the Laplace transformed Poisson Equation (5) is used, *<sup>∂</sup>* 2*ϕ ∂y* <sup>2</sup> = − *ρ ε* , and integrated once:

$$\frac{\partial \overline{\boldsymbol{\varphi}}}{\partial \boldsymbol{y}} = -\frac{A}{\varepsilon \boldsymbol{M}} \cosh \left( \boldsymbol{M} \left( \boldsymbol{y} - \frac{\boldsymbol{L}}{2} \right) \right) + \boldsymbol{\Gamma} \tag{A16}$$

where Γ is a constant to be determined by considering the current density, *i* = *zeNA*(*J*<sup>+</sup> − *J*−), and it vanishes at the duct walls. Using Equation (12), the current density can be found as

$$i = -\kappa^2 \varepsilon D \frac{\partial \rho}{\partial y} - D \frac{\partial \rho}{\partial y} \tag{A17}$$

While the Laplace transformed Equation (A17) reads as

$$\overline{\dot{a}} = -\kappa^2 \varepsilon D \frac{\partial \overline{\varphi}}{\partial y} - D \frac{\partial \overline{\varphi}}{\partial y} \tag{A18}$$

and *i*(*y* = 0) = 0. Substituting *ρ* and *∂ϕ ∂y* from Equations (A14) and (A16), respectively, the constant Γ can be found as

$$
\Gamma = \frac{AM}{\varepsilon} \left( M^{-2} - \kappa^{-2} \right) \cosh \left( M \frac{L}{2} \right) \tag{A19}
$$

Moreover, integrating Equation (A16) and using the antisymmetric condition, *ϕ*(*y*) = *<sup>ϕ</sup>*(*<sup>L</sup>* <sup>−</sup> *<sup>y</sup>*), we find *<sup>ϕ</sup>* <sup>=</sup> <sup>−</sup> *<sup>A</sup> εM*<sup>2</sup> sinh<sup>h</sup> *M <sup>y</sup>* <sup>−</sup> *<sup>L</sup>* 2 i <sup>+</sup> <sup>Γ</sup>(*<sup>y</sup>* <sup>−</sup> *<sup>L</sup>*), and from Equation (A19):

$$\overline{\varphi} = -\frac{A \cosh\left(M\frac{L}{2}\right)}{\varepsilon M^2} \left\{ \frac{\sinh\left[M\left(y - \frac{L}{2}\right)\right]}{\cosh\left[M\frac{L}{2}\right]} + \frac{MS\left(y - \frac{L}{2}\right)}{D\kappa^2} \right\} \tag{A20}$$

The constant *A* can be obtained by the boundary condition at *y* = 0 from the transformed Equation (6), *ϕ* = +*ϕ*(0)*s* <sup>−</sup><sup>1</sup> + *λ<sup>s</sup> ∂ϕ ∂y* , when *ϕ* and *∂ϕ ∂y* are substituted by Equations (A20) and (A16) for *y* = 0. Thus, *A* is found as

$$A = \frac{\varrho(0)s^{-1}\varepsilon M^2}{\cosh\left(M\frac{L}{2}\right)} \frac{1}{\frac{MSL}{2Dk^2}\left(1 + \frac{2\lambda\_s}{L}\right) + \lambda\_s M + \tanh\left(M\frac{L}{2}\right)}\tag{A21}$$

The Debye time is defined as *τ<sup>D</sup>* = *λ* 2 *D <sup>D</sup>* = *κ* <sup>2</sup>*D* <sup>−</sup><sup>1</sup> and is of the order *<sup>τ</sup><sup>D</sup>* <sup>∼</sup> <sup>10</sup>−<sup>7</sup> *s*, while we are interested in system responses at higher time scales. From the Laplace transformation, *f*(*s*) = R <sup>∞</sup> 0 *dτ e* −*sτ f*(*τ*), and it may be observed that for the function not to be zero as time increases, *S* should be less than a threshold value. Given that *τ* −1 *<sup>D</sup>* = *κ* <sup>2</sup>*D*, it is found that *S κ* <sup>2</sup>*D*, and from Equation (A12), it appears that *<sup>M</sup>* <sup>≈</sup> *<sup>κ</sup>*. Then,

$$A = \mathfrak{a} \frac{\mathbb{S}^{-1}}{1 + \mathfrak{r}\mathbb{S}} \tag{A22}$$

$$\text{where } \kappa = \frac{\varrho(0)\varepsilon\kappa^2}{\cosh\left(\kappa\frac{L}{2}\right)\left(\lambda\_s\kappa + \tanh\left(\kappa\frac{L}{2}\right)\right)}\tag{A23}$$

$$\text{rand } \tau = \frac{L\left(1 + \frac{2\lambda\_s}{L}\right)}{2D\kappa \left(\lambda\_s \kappa + \tanh\left(\kappa \frac{L}{2}\right)\right)}\tag{A24}$$

Thus,

$$\overline{\rho} = \kappa \sinh\left(\kappa \left(y - \frac{L}{2}\right)\right) \frac{\mathcal{S}^{-1}}{1 + \tau \mathcal{S}} \tag{A25}$$

$$\frac{\partial \overline{\rho}}{\partial y} = -\frac{a}{\varepsilon \kappa} \cosh \left( \kappa \left( y - \frac{L}{2} \right) \right) \frac{\mathcal{S}^{-1}}{1 + \tau \mathcal{S}} - a \frac{1}{1 + \tau \mathcal{S}} \frac{\cosh \left( \kappa \frac{L}{2} \right)}{\varepsilon D \kappa^3} \tag{A26}$$

$$\overline{\varphi} = -a \frac{\cosh\left(\kappa \frac{L}{2}\right)}{\varepsilon \kappa^2} \left\{ \frac{\sinh\left[\kappa \left(y - \frac{L}{2}\right)\right]}{\cosh\left[\kappa \frac{L}{2}\right]} \frac{\mathcal{S}^{-1}}{1 + \tau \mathcal{S}} + \frac{1}{1 + \tau \mathcal{S}} \frac{\left(y - \frac{L}{2}\right)}{D\kappa} \right\} \tag{A27}$$

The inverse Laplace transformation of functions is *<sup>S</sup>* −1 1+*τS* and <sup>1</sup> 1+*τS* is [*L*] −1 h *S* −1 1+*τS* i = 1 − *e* − *t <sup>τ</sup>* and [*L*] −1 h 1 1+*τS* i = <sup>1</sup> *τ e* − *t <sup>τ</sup>* , respectively, resulting then in Equations (14)–(21).
