Exploring the Impact of Steric Effects on Ion Removal of Water Solutions under the Influence of an Electric Field

Abstract

In this study, we examine the movement of ions that are in a water solution which flows along a duct, due to the existence of an electric field, taking into account the size of the ions, a phenomenon known as the steric effect. We compare the results from the above model with the classical one (the one that uses the Boltzmann distribution where ions are considered dimensionless) for various parameters such as surface charge density, electric field and differential capacitance. It is shown that for dilute water solutions (10¹⁹–10²⁴ ions/m³ final concentration at the center of the duct), with ions of valence $z=1$ (let us say saline water), steric effects become important for potentials greater than 1 V, and the phenomenon is more pronounced at higher concentrations. Furthermore, the steric effect model is applied to the calculation of the percentage of reduction in ion concentration in the main volume of the solution as a function of duct width for various electrode potentials and initial ion concentrations. Removal times are also calculated using Modified PNP equations which take into account steric effects. It is found that with a potential of 2.6 V, a 96% reduction in ions is achieved in the main volume of the solution for duct width 0.1 mm for 10²¹ ions/m³ final concentration at the center of the duct within approximately 1.6 s, while the percentage drops to 80% for duct width 1 mm. For smaller potentials, no noticeable decrease in concentration is observed, while for higher potentials, there are more impressive results, but we must be very careful because there is the case of other electrochemical phenomena taking place. The results are better when reducing the width of the duct, but relatively large widths are considered for the method to be practically applicable. With the increase in the concentration of the ions, their reduction percentage in the main volume of the solution decreases but remains significant up to 10²³ ions/m³ final concentration at the center of the duct. In addition, the completion time is shown to be proportional to the duct width. Therefore, for example, with the other parameters the same (2.6 V, 10²¹ ions/m³) but with $L~1\mathrm{mm}$, the completion time can be estimated to be approximately 16 s. This observation enables us to estimate the completion time for different duct widths, eliminating the need for repeated numerical computation of the MPNP equations.

1. Introduction

The recent rapid progression of climate change is causing significant disturbances to our planet, not only disrupting its natural balance but also foreboding more extensive suffering for humanity in the future. These disturbances include major natural disasters like floods, frequent heatwaves, extended wildfires, desertification of large areas and insufficient water for irrigation and essential human needs. Given these significant changes, the need for ensuring the availability of safe drinking water is more pressing than ever. Achieving this goal necessitates the removal of ions from various sources, including seawater, brackish water and industrial water contaminated with heavy metals. Numerous physical and chemical techniques have been developed for this purpose. The use of an electric field, which employs electrodes in either faradaic or non-faradaic processes, plays a pivotal role in these methods. A similar approach involving nanochannels and using molecular dynamics demonstrated that an electric field can successfully drift salt ions [1,2,3]. However, practical experience has shown that the application of an electric field alone yields limited results, especially when dealing with ducts of relatively large width and solutions containing high ion concentrations, as demonstrated by Bartzis et al. [4,5,6,7,8,9]. In such cases, it becomes essential to combine this method with the use of membranes or porous electrodes, a process known as capacitive deionization. This necessity arises from the formation of an Electric Double Layer (EDL).

The EDL is formed in front of a charged electrode by ions carrying a charge opposite to that of the electrode. These ions polarize the electrode, thereby inhibiting its ability to influence other ions in the solution. Until recently, the dominant theory for studying the EDL was based on the Stern model, which considers the presence of two distinct sub-layers within the EDL [10,11]. The first sub-layer, referred to as the Stern layer or compact layer, arises because ions, constrained by their effective radius, cannot approach the electrode any closer. So, taking into consideration this particular layer, it is a first effort to account for the finite dimensions of ions. In this layer, the potential decreases linearly (so-called “Graham ansatz”), which is valid only in the absence of specific adsorption, a point that has raised experimental questions [10]. In the second layer, known as the diffuse layer, originally proposed by Gouy and Chapman [10] for a dilute solution state where ion dimensions are considered negligible, ions are widely dispersed, and they follow a Boltzmann distribution. However, this is also an approximation, as real solutions may not always be dilute, and ions are separated by a minimum distance they cannot breach, as previously noted. Nevertheless, by combining these two layers that constitute the Electric Double Layer (EDL), an attempt is made to depict the solution’s behavior. When dealing with substantial electrode polarizations and solutions that are not very dilute, the compact Stern layer predominantly dictates the solution’s behavior. Otherwise, the diffuse layer is dominant.

It became evident rather swiftly that the diffuse layer model (GCh model) would prove inadequate, even with the inclusion of the compact layer, when the extent of the diffuse layer (Goey length) approaches the order of ions’ diameter [10,11]. In ordinary diluted solutions, electrochemical reactions occur not only between the ions but also between the ions and the solvent as the potential increases. These reactions take place before steric effects become significant, which is why steric effects were initially overlooked. However, the significance of these considerations gained prominence with the advent of research on room-temperature ionic liquids (RTILs) and microfluidics. The discussion around this issue has been reignited by the work of numerous authors who, acknowledging the steric effects of ions, have proposed a distribution that differs from Boltzmann’s and uniquely characterizes the diffuse layer [10,11,12,13,14,15,16,17,18].

In our previous research, we studied the behavior of water solutions containing ion salts, heavy metal ions, toxins, viruses and anthocyanins, and we relied on the Stern model [4,5,6,7,8,9], where the solutions were flowing through ducts exposed to electric field forces generated by an external capacitor. We also calculated the specific potential and concentrations at which this model produced reliable quantitative results.

In this current study, we study the drift of ions in an aqueous solution which, as in previous research, flows along a duct and is exposed to an electric field, taking into account the ion size (steric effects), leading to a new distribution of ions (Modified Boltzmann distribution) [16,17,18,19,20,21,22,23,24,25,26,27,28,29]. In Section 4, we focus on the final state of the solution by calculating various parameters such surface charge density, electric field and differential capacitance, as well as ion concentration. First of all, we compare the values of the above parameters with those obtained when we consider the ions to be dimensionless (Boltzmann distribution), examining the compatibility of the two analyses and at which values of concentration and electrode potential the steric effects become significant. Although in the analysis, the theory developed is general, the application is for dilute water solutions of ions with valence $z=1$ and, in particular, for salt ions (e.g., saline water). In the same section, we also examine the percentage reduction in ion concentration in the main volume of the solution for potentials and initial concentrations where steric effects become significant, as a function of duct width for various values of potentials and concentrations. In Section 5, using the Modified PNP (Poisson Nernst Planck) equations, we study the temporal evolution of ion drift, giving an estimation for the completion times of the ion drift as a function of the initial concentration and the electrode potential. Finally, in Section 6, we comment on the results.

The mean field approximation for electrostatic interactions is utilized, which considers the field as collectively determined by the solvent and ions. The question then arises: can we confidently adopt this approach? This is a matter of ongoing debate, and the final decision will likely be based on experimental data. However, this consideration is closely related to the concept of the Bjerrum length, which is defined as the distance at which the Coulomb energy between two elementary charges in a medium with a macroscopic dielectric constant ε equals the thermal energy $kT$.

$$L_B=\frac{e^2}{4\pi \epsilon kT}$$

In an aquatic environment, this distance is typically around 7 Å [16]. When the distance exceeds this, it can be inferred that the electrical interaction energy is less than the thermal energy. This allows us to assume that the mean field approximation is valid.

2. Theoretical Description of the Proposed Experiment

The proposed experimental device is based on our previous work [4,7,9] and consists of the following:

(a) A capacitor with its two electrodes at voltage V, as is shown in Figure 1 and Figure 2, which create between them an almost spatially constant electric field E, directed from the positive to the negative electrode. The distance between the electrodes is L.

(b) The main duct into which the water fluid flows, which is placed along the electrodes, as shown in Figure 1 and Figure 2. Thus, the velocity of flow $\overrightarrow{\upsilon }$ is perpendicular to the intensity of the electric field E, which is along y-axis (Figure 2) and, at first approximation, is considered constant over the entire range of the fluid. Due to the electrical force, as will be discussed extensively below, the charged ions inside the water solution will move towards the oppositely charged electrodes, creating an accumulation of ions near the side walls of the duct.

(c) Two smaller ducts near the side walls (Figure 2) in order to collect excess ions. The gap between the duct and the electrodes is kept as small as possible to ensure the spatial stability of the electric field. Figure 3 depicts the accumulation of the ions near the walls of the duct.

Current research aligns with the Debye–Huckel theory [30,31,32,33], which regards water as a continuous dielectric medium. Consequently, water influences the solution’s electrical permittivity and changes it to $\mathsf{\epsilon }={\mathsf{\epsilon }}_{\mathrm{r}}{\mathsf{\epsilon }}_0$, where ${\mathsf{\epsilon }}_{\mathrm{r}}\approx 80$ is the relative permittivity of the water, and ${\mathsf{\epsilon }}_0=8.85\times {10}^{-12}\mathrm{F}/\mathrm{m}$.

3. Creation of the Modified Stern Model

3.1. Description of the Theoretical Model

Considering that ions are not point-like but have dimensions in reality, we can treat them as spheres with a constant effective diameter $\alpha$ during their movement. Consequently, when they accumulate near the duct walls, they can approach a maximum distance ${\lambda }_s=\alpha /2$, equal to their effective radius, thereby forming the Stern layer. Section 3.2 will provide a more detailed discussion about ${\lambda }_s$. In summary, a double layer is formed, which, according to Stern’s theory, consists of two consecutive sub-layers, as depicted in Figure 4.

(a) The compact layer is in contact with the duct walls, which is also known as the Stern layer. This layer is represented by an effective capacitor of width ${\lambda }_S$ with capacitance per unit area:

$$c_H=\frac{\epsilon }{{\lambda }_S}$$

(1)

(b) The diffuse layer next to the compact layer is represented by an effective capacitor with capacitance per unit area ${\mathrm{c}}_{\mathrm{D}}$:

$$c_D=\frac{\epsilon }{\lambda }$$

(2)

where λ is its effective width of the diffuse layer (Figure 4).

The differentiation from the Gouy–Chapman theory occurs within this layer. We do not consider the solution to be dilute enough, nor do we consider the ions as points when studying the distribution of charges based on the Boltzmann distribution. In other words, we take into account the dimensions of the ions, as previously mentioned, and therefore, the concentration can reach a maximum value. In order to define the boundary conditions, the potential is considered zero at the center of the conductor $\phi (y=L/2)=0$, and therefore, at positions $y=0$ and $y=L$, the potentials are $\phi (0)$ and $-\phi (0)$, respectively. Within the compact layer, the potential is assumed to vary linearly with distance, so

$${\pm \phi }_S=\pm \phi \left(0\right)\pm {\lambda }_s\frac{\partial \phi }{\partial y}\mathrm{for}y=0,L$$

(3)

where ${\phi }_S$ is the potential at the Outer Helmhotz Plane (OHP) at distance ${\lambda }_s$ (Figure 4).

3.2. Final State Parameters

As previously mentioned, we aim to develop a theory that considers the size of ions. Therefore, initially, we need to define certain quantities that are related to the size of the ions and will be useful later on. Following the analysis of Kilic et al. [16,17] and initially examining the case of the symmetric electrolyte $z:z$ with ions of the same size, we will begin by discussing some characteristic lengths that need to be defined. The first is Debye length:

$${\lambda }_D=\sqrt{\frac{\epsilon kT}{2z^2{C}_M{e}^2}}$$

(4)

where $k$ is the Boltzmann constant, $T$ is the absolute temperature (it is considered $T=300\mathrm{K}$ throughout this study), $z$ is the number of overflow protons or electrons, a quantity that is considered to be positive throughout this paper, and $C_M\left(\frac{\mathrm{i}\mathrm{o}\mathrm{n}}{{\mathrm{m}}^3}\right)$ is the positive and negative separately ion concentration at the center of the duct, where we consider the potential to be zero. The elementary charge is $e=1.6\times {10}^{-19}\mathrm{C}$.

The second characteristic length is the average distance $L_{bef}$ of ions in the solution before applying the electric field, which can be calculated by assuming complete ionic decay:

$$2C_{bef}=\frac{ions}{V}=\frac{1}{{\displaystyle L_{bef}^3}}\to L_{bef}={\left(2C_{bef}\right)}^{-\frac{1}{3}}$$

(5)

where $C_{bef}$ is the positive and negative separately constant ion concentration before the application of the electric field.

Similarly, we can calculate the average distance $L_{{\rm M}}$ of ions in the middle of the conductor:

$$2C_{{\rm M}}=\frac{ions}{V}=\frac{1}{{\displaystyle L_{{\rm M}}^3}}\to L_{{\rm M}}={\left(2C_{{\rm M}}\right)}^{-\frac{1}{3}}$$

(6)

Next, we will establish a connection between the dimensions of the ion and the maximum density achievable. If we consider $\alpha$ to be the effective diameter of the ion, we have the following:

$$C_{max}=\frac{1}{a^3}$$

(7)

Here, we should emphasize that the diameter should not be completely equated with the effective diameter. In most ions, the diameter is of the order of 1 Å; there may be a solvation shell, which can extend it to 3 Å. If we also consider the Bjerrum length, which, as we mentioned, is at least in the order of 7 Å, it increases even further. Let us consider it to be as large as the Bjerrum length, i.e., $\alpha =7\mathrm{\AA }$. Especially in the Stern layer, due to the large fields, the water dipoles are highly aligned. As a result, the relative susceptibility decreases from, let us say, $80{\epsilon }_0$ to about $8{\epsilon }_0$, changing the Bjerrum length to 7 nm. Therefore, based on what we have already mentioned, the effective width of the Stern (compact) layer is ${\lambda }_S=3.5\mathrm{n}\mathrm{m}$. The width of the conductor L is considered to be much larger than the Debye length, so ${\lambda }_D\ll L$.

From these characteristic quantities, we can define three very useful dimensionless quantities. One is the $\alpha /{\lambda }_D$, which is small for ordinary electrolytes and relatively dilute solutions. The other is the fraction of the volume occupied by ions in the solution medium:

$$\nu =\frac{V_{ion}}{V_{L_M}}=\frac{a^3}{{\displaystyle L_M^3}}\to \nu =2C_M{a}^3$$

(8)

In a similar way, we can define the fraction of the volume of ions in the initial solution before the application of the field:

$${\nu }^{\prime }=\frac{V_{ion}}{V_{bef}}=\frac{a^3}{{\displaystyle L_{bef}^3}}\to {\nu }^{\prime }=2C_{bef}{a}^3$$

(9)

4. Final State Study

The phenomenological equation of free energy can be considered to be of the following form:

$$F=U-TS$$

(10)

where

$$U=\int dr(-\frac{\epsilon }{2}{\left|\nabla \phi \right|}^2+ze{C}_{+}\phi -ze{C}_{-}\phi )$$

(11)

and

$$-TS=\frac{kT}{a^3}\int dr[C_{+}{a}^3\ln \left(C_{+}{a}^3\right)+C_{-}{a}^3\ln \left(C_{-}{a}^3\right)+\left(1-C_{+}{a}^3-C_{-}{a}^3\right)ln\left(1-C_{+}{a}^3-C_{-}{a}^3\right)]$$

(12)

In Equation (11), the first term represents the self-energy of the electric field created by a potential difference applied to the boundaries of the output. The next two terms represent the electrostatic energy of the ions. The entropic contribution $-TS$ in Equation (12) is due to steric effects. The first and second terms represent the entropies of positive and negative ions, respectively, which we have assumed to be of the same size $a$. The last term accounts for the possibility of large concentrations, as it represents the entropy of solvent molecules.

Following the method of Borukhov et al. [14,15] we obtain the following:

$$C_{+}=\frac{C_M{e}^{-\frac{ze\phi }{kT}}}{1+2v{sinh}^2\left(\frac{ze\phi }{2kT}\right)}$$

(13)

$$C_{-}=\frac{C_M{e}^{+\frac{ze\phi }{kT}}}{1+2v{sinh}^2(\frac{ze\phi }{2kT})}$$

(14)

where $C_{\pm }\left(\mathrm{i}\mathrm{o}\mathrm{n}·{\mathrm{m}}^{-3}\right)$ are the concentrations of the positive and negative ions, and $\nu$ is given by Equation (8) (Modified Boltzmann distribution).

Therefore, the charge density corresponding to the above ion concentrations is given by the following relation:

$$\rho =C_{+}ze-C_{-}ze=-2ze{C}_M\frac{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\frac{ze\phi }{kT})}{1+2\nu {sinh}^2(\frac{ze\phi }{2kT})}$$

(15)

where the Poisson equation takes the following form:

$$\frac{{\partial }^2\phi }{\partial y^2}=-\frac{\rho }{\epsilon }$$

$$\frac{{\partial }^2\phi }{\partial y^2}=\frac{2ze{C}_M}{\epsilon }\frac{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\frac{ze\phi }{kT})}{1+2\nu {sinh}^2(\frac{ze\phi }{2kT})}$$

(16)

where $C_M$ is the concentration of positive and negative ions separately at the center of the duct.

The calculation of the electric field in the final equilibrium state can be performed as follows:

$$\mathrm{E}=\frac{-\partial \mathsf{\phi }}{\partial \mathrm{y}}\to \frac{\vartheta E}{\vartheta y}=-\frac{{\vartheta }^2\phi }{{\vartheta y}^2}\to \frac{\partial \mathrm{E}}{\partial \mathsf{\phi }}\frac{\partial \mathsf{\phi }}{\partial \mathrm{y}}=\frac{\mathsf{\rho }}{\mathsf{\epsilon }}\to -2E\frac{\vartheta E}{\vartheta \phi }=\frac{2\rho }{\epsilon }$$

$$d(E^2)=-\frac{2\rho }{\epsilon }d\phi$$

Assuming $\phi (y=\frac{L}{2})=0$ and $E_M$ the electric field at the center of the duct (it is set approximately equal to zero since it has a very small value compared to that near the walls of the duct), we have the following:

$${\displaystyle {\int }_{{{\rm E}}_{{\rm M}}}^{{\rm E}}}d\left(E^2\right)=-{\displaystyle {\int }_0^{\phi }}\frac{2\rho }{\epsilon }d\phi \to$$

$$E=-\frac{\partial \phi }{\partial y}\approx \sqrt{\frac{4C_{M}kT}{\epsilon v}\mathit{ln}[1+2\nu {\mathit{sinh}}^2(\frac{ze}{2kT}\phi )]}$$

(17)

The previous equation is only valid within the area of the diffuse layer, assuming that $E_M$ is considered zero. Near the positive electrode, we observe the following:

$$E=-\frac{\partial \phi }{\partial y}\approx \sqrt{\frac{4C_{M}kT}{\epsilon v}\ln [1+2\nu {\mathit{sinh}}^2(\frac{ze}{2kT}\phi )]}\phi >0$$

(18)

for

$${\lambda }_s\le y\le \raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$$

Integrating the Poisson Equation (16) we can calculate the surface charge density $\sigma (\frac{\mathrm{C}}{{\mathrm{m}}^2})$ from ${y=\lambda }_S$ to $y=\frac{L}{2}$ as follows:

$${\displaystyle {\int }_{{\lambda }_S}^{\frac{L}{2}}}\frac{d^2\phi }{d{y}^2}dy=-{\displaystyle {\int }_{{\lambda }_S}^{\frac{L}{2}}}\frac{\rho }{\epsilon }dy\to {\left.\frac{d\phi }{dy}\right|}_{\frac{L}{2}}-{\left.\frac{d\phi }{dy}\right|}_{{\lambda }_S}=-\frac{{\displaystyle {\int }_{{\lambda }_S}^{L/2}}\rho dy}{\epsilon }=-\frac{\sigma }{\epsilon }$$

Considering that $E_M=-{\left.\frac{d\phi }{dy}\right|}_{\frac{L}{2}}=0$ we have $-{\left.\frac{d\phi }{dy}\right|}_{{\lambda }_S}=-\frac{\sigma }{\epsilon }\to \sigma =-\epsilon {{\rm E}}_{{\lambda }_S}$ which due to Equation (18) becomes the following:

$$\sigma =-\sqrt{\frac{4kT\epsilon C_M}{\nu }ln[1+2\nu {sinh}^2(\frac{ze{\phi }_s}{2kT})]}$$

(19)

At this point, we should emphasize that due to the principle of continuity, the electric field is equal on both sides of the OHP, i.e., both on the side of the compact layer and on the side of the diffuse layer. Furthermore, considering its stability within the compact layer, we derive for the area near the positive electrode from Equation (3) the following:

$${\phi }_S=\phi \left(0\right)+{\lambda }_s{\left.\frac{d\phi }{dy}\right|}_{{\lambda }_S}=\frac{\sigma {\lambda }_s}{\epsilon }+\phi \left(0\right)$$

Or, using it in Equation (19),

$$\sigma =-\sqrt{\frac{4kT\epsilon C_M}{\nu }ln[1+2\nu {sinh}^2(\frac{ze}{2kT}\left(\frac{\sigma {\lambda }_s}{\epsilon }+\phi \left(0\right)\right))]}$$

(20)

and

$${\phi }_S=\phi \left(0\right)-{\lambda }_s\sqrt{\frac{4kT{C}_M}{\epsilon \nu }ln[1+2\nu {sinh}^2(\frac{ze{\phi }_s}{2kT})]}$$

(21)

Figure 5 shows the relationship of ${\phi }_S$ versus $\phi (0)$ for various concentrations, as determined by Equation (21), which takes into account steric effects (Modified Poisson Boltzmann PB model), in comparison with the results derived from the corresponding equation that utilizes the Poisson Boltzmann (PB) model, as elaborated in references [4,7]:

$${\phi }_S=\phi (0)-{\lambda }_s\sqrt{\frac{8C_{M}kT}{\epsilon }}\sinh (\frac{ze}{2kT}{\phi }_s)$$

(22)

for ${\phi }_S\le {\displaystyle {\phi }_S^{max}}=\frac{\kappa T}{ze}\mathrm{l}\mathrm{n}(\frac{1}{{{\alpha }^{3}C}_M})$.

We observe the complete match of the two analyses for potentials $\mathsf{\phi }\left(0\right)\le 1\mathrm{V}$, while the exact point where the deviation begins depends on the concentration $C_M$. The maximum ${\displaystyle {\phi }_S^{max}}$ that the PB distribution holds is given by Equation (22) [7] (${\displaystyle {\phi }_S^{max}}=0.5\mathrm{V}$ for $C_M={10}^{19}\frac{\mathrm{i}\mathrm{o}\mathrm{n}}{{\mathrm{m}}^3}$ and ${\displaystyle {\phi }_S^{max}}=0.2\mathrm{V}$ for $C_M={10}^{24}\frac{\mathrm{i}\mathrm{o}\mathrm{n}}{{\mathrm{m}}^3}$), which is again in full agreement with our results.

Another observation that can be made is that steric effects reduce the maximum concentration of ions near the electrode and, therefore, increase the value of ${\phi }_s$ for the same $\phi (0)$, for values of $\phi (0)$ greater than 1 V.

The same conclusion follows from the diagram of surface charge density $\sigma$ (Equation (20)) as a function of $C_M$, depicted in Figure 6. In the same figure, we have incorporated the prediction into the linear regime that holds when $z\phi (0)\le 0.026\mathrm{V}$ for which the following relation holds [4]:

$$\left|{\sigma }_l\right|=\phi \left(0\right)\frac{\epsilon }{({\lambda }_S+{\lambda }_D)}$$

(23)

We also observe a complete alignment of the Modified PB model with the PB model [7]:

$$\sigma =-\sqrt{8{\epsilon C}_{M}kT}\sinh [\frac{ze}{2kT}(\frac{\sigma {\lambda }_s}{\epsilon }+\phi \left(0\right))]$$

(24)

again for $\phi \left(0\right)\le 1\mathrm{V}$, while for higher values, a deviation is observed. More specifically, for $\phi \left(0\right)=0.026\mathrm{V}$, there is a complete match of the three models (linear, PB Distribution, Modified PB Distribution). For $\phi \left(0\right)=0.26\mathrm{V}$, there is complete alignment between both the PB and Modified PB distributions. However, for $\phi \left(0\right)=2.6\mathrm{V}$, the divergence between the two models gradually begins. It appears that the inclusion of steric effects reduces the prediction of charge concentration, and as $\phi \left(0\right)$ increases, the deviation also increases.

The potential distribution as a function of y can be calculated for the region near the positive electrode as below:

$$\frac{d\phi }{dy}=-\sqrt{\frac{4C_{M}kT}{\epsilon v}\mathit{ln}[1+2\nu {\mathit{sinh}}^2(\frac{ze}{2kT}\phi )]}\phi >0\mathrm{for}{\lambda }_s\le y\le \raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$$

So,

$${\displaystyle {\int }_{{\phi }_s}^{\phi }}\frac{d\phi }{\sqrt{\mathit{ln}[1+2\nu {\mathit{sinh}}^2(\frac{ze}{2kT}\phi )]}}=-\sqrt{\frac{4C_{M}kT}{\epsilon v}}{\displaystyle {\int }_{{\lambda }_s}^y}dy\to$$

$$y={\lambda }_s+\sqrt{\frac{\epsilon v}{4C_{M}kT}}{\displaystyle {\int }_{\phi }^{{\phi }_s}}\frac{d\phi }{\sqrt{\mathit{ln}[1+2\nu {\mathit{sinh}}^2(\frac{ze}{2kT}\phi )]}}\mathrm{for}{\lambda }_s\le y\le \raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$$

(25)

while for the region near the negative electrode, we have

$$\frac{d\phi }{dy}=-\sqrt{\frac{4C_{M}kT}{\epsilon v}\mathit{ln}[1+2\nu {\mathit{sinh}}^2(\frac{ze}{2kT}\phi )]}\phi <0\mathrm{for}\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\le y\le L-{\lambda }_s$$

(26)

So,

$$y={L-\lambda }_s-\sqrt{\frac{\epsilon v}{4C_{M}kT}}{\displaystyle {\int }_{{-\phi }_s}^{\phi }}\frac{d\phi }{\sqrt{\mathit{ln}[1+2\nu {\mathit{sinh}}^2(\frac{ze}{2kT}\phi )]}}\phi <0\mathrm{for}\raisebox{1ex}{$\mathrm{L}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\le \mathrm{y}\le \mathrm{L}-{\mathsf{\lambda }}_{\mathrm{s}}$$

(27)

Next, we will calculate the total differential capacitance which is defined as follows:

$$c_{tot}=\frac{\left|d\sigma \right|}{d\phi (0)}$$

(28)

which takes the following form:

$$c_{tot}=\frac{d\sigma }{d\phi (0)}=\frac{d\sigma }{{d\phi }_S}\frac{{d\phi }_S}{d\phi (0)}$$

Calculating the first fraction, we obtain

$$A=\frac{d\sigma }{{d\phi }_S}=-\sqrt{\frac{\epsilon C_M\nu z^2{e}^2}{kT}}\frac{\sinh (\frac{\mathrm{z}\mathrm{e}{\mathsf{\phi }}_{\mathrm{s}}}{\mathrm{k}\mathrm{T}})}{\left(1+2\nu {sinh}^2\left(\frac{ze{\phi }_s}{2kT}\right)\right)\sqrt{\ln \left[1+2\nu {sinh}^2\left(\frac{ze{\phi }_s}{2kT}\right)\right]}}$$

and because of Equation (3), we derive for the second fraction:

$$\frac{{d\phi }_S}{d\phi \left(0\right)}=\frac{{\lambda }_s}{\epsilon }\frac{d\sigma }{d\phi \left(0\right)}+1$$

So by combining the two previous equations, we have:

$$c_{tot}=\frac{\left|d\sigma \right|}{d\phi (0)}=\frac{\left|A\right|}{1+\frac{{\lambda }_s\left|A\right|}{\epsilon }}\equiv \frac{1}{c_H}+\frac{1}{c_D}$$

(29)

where $c_H$, $c_D$ are also defined by Equations (1) and (2) (see also Figure 4). So, we have

$$c_H=\frac{\epsilon }{{\lambda }_S}$$

(30)

$$c_D=\left|A\right|=\sqrt{\frac{\epsilon C_M\nu z^2{e}^2}{kT}}\frac{\mathit{sinh}(\frac{ze{\phi }_s}{kT})}{(1+2\nu {sinh}^2\left(\frac{ze{\phi }_s}{2kT}\right))\sqrt{ln[1+2\nu {sinh}^2(\frac{ze{\phi }_s}{2kT})]}}$$

(31)

Expressing the $c_D$ based on ${\lambda }_D$, we obtain the following:

$$c_D=\frac{\epsilon }{{\lambda }_D}\mathit{cosh}(\frac{ze{\phi }_s}{2kT})\frac{1}{(1+2\nu {sinh}^2\left(\frac{ze{\phi }_s}{2kT}\right))}\sqrt{\frac{2\nu {sinh}^2\left(\frac{ze{\phi }_s}{2kT}\right)}{ln[1+2\nu {sinh}^2(\frac{ze{\phi }_s}{2kT})]}}$$

(32)

Using Equation (2), $c_D=\frac{\epsilon }{\lambda }$, where λ is the effective thickness of the effective capacitor of the diffuse layer, we have

$$\frac{\lambda }{{\lambda }_D}={[\mathit{cosh}(\frac{ze{\phi }_s}{2kT})\frac{1}{\left(1+2\nu {sinh}^2\left(\frac{ze{\phi }_s}{2kT}\right)\right)}\sqrt{\frac{2\nu {sinh}^2\left(\frac{ze{\phi }_s}{2kT}\right)}{\ln \left[1+2\nu {sinh}^2\left(\frac{ze{\phi }_s}{2kT}\right)\right]}}]}^{-1}$$

(33)

Figure 7, which is based on Equation (33) in combination with Equation (21), represents the dependence of $\frac{\lambda }{{\lambda }_D}$ on $\phi (0)$. To understand the figure, the terms of Equation (33) have to be analyzed. The first term defines Gouy–Chapman behavior, while the other two are derived from the Modified PB Model. When $\nu <1/3$, as in our case of dilute solutions ($7\times {10}^{-9}<\nu <7\times {10}^{-5}$), we have inverse camel-shape curves. The term camel-shape was introduced by Kornyshev [10,11] to describe the shape of the graph $\raisebox{1ex}{$C$}\!\left/ \!\raisebox{-1ex}{$C_D$}\right.=f(\phi \left(0\right))$ when $\nu <1/3$. Of course, in our graphs, the curves are inverse because they represent $\frac{\lambda }{{\lambda }_D}=f(\phi \left(0\right))$. For small values of $v$ ($\nu <1/3$) and for low voltages, the system adheres to the Gouy–Chapman distribution, also known as the PB distribution. However, when the voltage exceeds $\phi \left(0\right)>1\mathrm{V}$, steric effects start to influence the system, causing the width of the diffuse layer to expand, which subsequently leads to a reduction in capacity. It is worth noting that these findings are entirely consistent with the results presented in [10,11,16].

The area of the duct that is of primary interest for practical applications is the one near its center, where the ion concentration, denoted as $C_M$, reaches its minimum value. Therefore, if the objective is to obtain a liquid with the lowest possible ion concentrations, the liquid should be extracted from this central area of the duct. It is important to note that the concentration remains practically constant in a broad region around the center of the duct. Conversely, higher ion concentrations are found near the duct walls.

However, since the concentration $C_M$ also depends on the ion concentration in the solution prior to the application of the electric field, the ratio $\frac{C_M}{C_{bef}}$ has to be calculated. The calculation is based on negative ions, which are concentrated near the positive electrode, bearing in mind that the situation is completely symmetrical for positive ions, which are concentrated near the negative electrode. As previously mentioned, the concentration of negative ions outside the compact layer as a function of the potential $\phi$ is given by Equation (14):

$$C_{-}=\frac{C_M{e}^{+\frac{ze\phi }{kT}}}{1+2v{sinh}^2(\frac{ze\phi }{2kT})}$$

where $\phi$ is given by Equations (25) and (27).

Using the distribution above, we can calculate $C_{bef}$ as follows:

$$C_{bef}=\frac{C_M{\displaystyle {\int }_{{\lambda }_s}^{L-{\lambda }_s}}\frac{e^{+\frac{ze\phi }{kT}}}{1+2v{sinh}^2(\frac{ze\phi }{2kT})}dy}{L}$$

(34)

So $\frac{C_M}{C_{bef}}$ has the following form:

$$\frac{C_M}{C_{bef}}=\frac{L}{{\displaystyle {\int }_{{\lambda }_s}^{L-{\lambda }_s}}\frac{e^{+\frac{ze\phi }{kT}}}{1+2v{sinh}^2(\frac{ze\phi }{2kT})}dy}$$

(35)

Figure 8 shows the $\frac{C_M}{C_{bef}}$ ratio as a function of the duct width $L$ for various $C_M\left(\frac{\mathrm{i}\mathrm{o}\mathrm{n}}{{\mathrm{m}}^3}\right)$ and for various potentials $\phi \left(0\right)$, $z=1$, with Figure 8a–c utilizing a different target concentration $C_M$. A general observation is that with increasing $C_M$, the $\frac{C_M}{C_{bef}}$ ratio increases, i.e., the ion removal efficiency decreases. The second observation is that as the duct width L increases, the ion removal efficiency decreases. Both results are logical, since it is expected that by increasing $\phi \left(0\right)$, we can achieve better purification ($\frac{C_M}{C_{bef}}$ decreases).

For practical applications to yield workable results, the duct should be somewhat large, i.e., of the order ${10}^{-4}\mathrm{m}\le L\le {10}^{-3}\mathrm{m}$. From Figure 8a, where $C_M={10}^{21}\raisebox{1ex}{$\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$, it is calculated that $\frac{C_M}{C_{bef}}=0.0386$ for potential $\phi \left(0\right)=2.6\mathrm{V}$ and $L={10}^{-4}\mathrm{m}$, so 96% of the charge has been removed from the central area of the duct. For $L={10}^{-3}\mathrm{m}$, $\frac{C_M}{C_{bef}}=0.2$, i.e., 80% of the charge has been removed from the central area of the duct. Respectively, for $C_M={10}^{22}\raisebox{1ex}{$\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$ (Figure 8b) and for potential $\phi \left(0\right)=2.6\mathrm{V}$, we take that $\frac{C_M}{C_{bef}}=0.2$ for $L={10}^{-4}\mathrm{m}$ (80% reduction) and $\frac{C_M}{C_{bef}}=0.8$ for $L={10}^{-3}\mathrm{m}$ (20% reduction). Figure 8c shows that for $C_M={10}^{23}\raisebox{1ex}{$\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$ and potential $\phi \left(0\right)=2.6\mathrm{V}$, the ratio $\frac{C_M}{C_{bef}}=0.7$ for $L={10}^{-4}\mathrm{m}$ (30% reduction), and negligible reduction is observed for $L={10}^{-3}\mathrm{m}$.

Smaller potentials for similar concentrations $C_M$ give a negligible decrease in concentration. For larger potentials (e.g., $\phi \left(0\right)=26\mathrm{V}$), we obtain impressive results; for example, when $L={10}^{-3}\mathrm{m}$ and for $C_M={10}^{21}\raisebox{1ex}{$\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$, we obtain a value of $\frac{C_M}{C_{bef}}\cong 0.07$ (93% reduction), for $C_M={10}^{22}\raisebox{1ex}{$\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.,\frac{C_M}{C_{bef}}\cong 0.3$ (70% reduction), for $C_M={10}^{23}\raisebox{1ex}{$\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.,\frac{C_M}{C_{bef}}\cong 0.8$ (20% reduction). However, in the case of such elevated potentials, other electrochemical phenomena can affect the process and should be adequately considered in the model [10].

5. Modified PNP Equations

In this section, we will extract the Modified Nernst Planck equations in order to study the temporal evolution of ion drift. Starting from the phenomenological equation of free energy, as it has been already presented in Section 4 (Equations (10)–(12)), the electrochemical potential ${\mu }_{+}$ of positive ions is calculated as follows [14,15,17]:

$${\mu }_{+}=\frac{\partial F}{\partial C_{+}}=ze\phi +kT[\ln \left(C_{+}{\alpha }^3\right)-\ln \left(1-C_{+}{\alpha }^3-C_{-}{\alpha }^3\right)]$$

(36)

The y axis ion flux of the positive ions is given by the following relation:

$$J_{+}=-\frac{D_{+}}{kT}ze{C}_{+}\frac{\vartheta \phi }{\vartheta y}-D_{+}\frac{\partial C_{+}}{\partial y}-\frac{D_{+}{\alpha }^3{C}_{+}}{1-\left(C_{+}+C_{-}\right){\alpha }^3}\frac{\vartheta }{\vartheta y}(C_{+}+C_{-})$$

(37)

where $z$ represents the number of excess protons. The elementary charge is $e=1.6\times {10}^{-19}\mathrm{C}$. The mobility of positive ions, denoted by $b_{+}\equiv \frac{1}{6\pi {\nu }_{v}r}$, depends on the dynamic viscosity of water ${\nu }_v$ and the effective radius of ions $r;$ $\phi$ signifies the electric potential within the fluid’s bulk. The diffusion coefficient of the positive ion is calculated by $D_{+}=b_{+}kT$, where k stands for the Boltzmann constant ($k\equiv \frac{R}{N_A}$, $R=8.314\mathrm{J}/(\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{K})$, ${{\rm N}}_A$ is the Avogadro constant), and T represents the absolute temperature (considered as T = 300 K throughout this study).

In the same way, we have for the electrochemical potential and the flux of the negative ions,

$${\mu }_{-}=\frac{\partial F}{\partial C_{-}}=-ze\phi +kT[\ln \left(C_{-}{\alpha }^3\right)-\ln \left(1-C_{+}{\alpha }^3-C_{-}{\alpha }^3\right)]$$

(38)

$$J_{-}=\frac{D_{-}}{kT}ze{C}_{-}\frac{\vartheta \phi }{\vartheta y}-D_{-}\frac{\partial C_{-}}{\partial y}-\frac{D_{-}{\alpha }^3{C}_{-}}{1-\left(C_{+}+C_{-}\right){\alpha }^3}\frac{\vartheta }{\vartheta y}(C_{+}+C_{-})$$

(39)

In Equations (37) and (39) above, the electric field is replaced by ${\rm E}=-\frac{\partial \phi }{\partial y}$. Assuming both positive and negative ions possess identical $z$ values and sizes, their respective mobilities are also equal, $b_{+}=b_{-}=b$, consequently making diffusion coefficients equivalent as well: $D_{+}=D_{-}=D$.

So, the Modified Nernst Planch equations are as follows:

$$\frac{\partial C_{\pm }}{\partial t}=-\frac{\partial }{\partial y}\left(J_{\pm }\right)$$

(40)

where $J_{\pm }$ are given by Equations (37) and (39).

The charge density, resulting from the concentration of the ions, is as follows:

$$\rho =C_{+}ze-C_{-}\mathrm{z}e=ze\left(C_{+}-C_{-}\right)$$

By applying the 1D Modified Poisson equation in the y-direction of the duct, it is found that

$$\frac{{\partial }^2\phi }{\partial y^2}=-\frac{\rho }{\epsilon }\to \frac{{\partial }^2\phi }{\partial y^2}=-\frac{ze(C_{+}-C_{-})}{\epsilon }$$

(41)

Equations (40) and (41) form the Modified Poisson Nernst Planck equations (MPNP equations).

5.1. Boundary Conditions

Prior to solving the MPNP equations, it is essential to establish the suitable boundary conditions. At the duct’s midpoint $\frac{L}{2}$, the electric potential remains zero $\phi (y=\frac{L}{2})=0$. The potentials at y = 0 and $y=L$ are denoted as $\phi \left(0\right)$ and $-\phi \left(0\right)$, respectively. Therefore, assuming a linear relationship between distance and potential in the compact part of the double layer adjacent to the electrodes, the equation becomes

$$\phi =\pm \phi (0)\pm {\lambda }_s\frac{\partial \phi }{\partial y}\mathrm{for}y={\lambda }_S,L-{\lambda }_S$$

(42)

As always, ${\lambda }_s$ represents the effective thickness of the compact part of the double layer [19,20]. Additionally, given that this process is non-Faradaic [16,17,20], indicating no transfer of charge through the electrodes, both ion fluxes and the current density must equate to zero at the boundaries, that is,

$$J_{\pm }=D\frac{\partial C_{\pm }}{\partial y}\pm \frac{D}{kT}ze{C}_{\pm }\frac{\vartheta \phi }{\vartheta y}+\frac{D{\alpha }^3{C}_{\pm }}{1-\left(C_{+}+C_{-}\right){\alpha }^3}\frac{\vartheta }{\vartheta y}\left(C_{+}+C_{-}\right)=0\mathrm{for}y=0,L$$

(43)

$$i=ze{N}_A\left(J_{+}-J_{-}\right)=0\mathrm{for}y=0,L$$

(44)

5.2. Initial Conditions

The model’s initial conditions involve implementing a uniform ion distribution for $t=0$, accompanied by a zero potential.

$$C_{+}\left(y,0\right)={\displaystyle C_{bef}^{+}}=C_{bef}{C}_{-}\left(y,0\right)={\displaystyle C_{bef}^{-}}=C_{bef}$$

(45)

At $t=0$, the potential has the distribution

$$\phi =\phi \left(0\right)\left(1-\frac{2}{y}L\right)$$

(46)

5.3. Dimensionless Formulation and Numerical Solution

In order to solve the system of Equations (40) and (41) numerically, taking into account the boundary conditions (42)–(44) and initial conditions (45) and (46), these equations are transformed into a dimensionless form by applying the following size definitions:

$$\tilde{y}=\frac{2y}{L},\tilde{t}=\frac{2tD}{{\displaystyle {\lambda }_D^{\prime }}L},{\displaystyle {\lambda }_D^{\prime }}=\sqrt{\frac{\epsilon kT}{2z^2{C}_{bef}{e}^2}},\tilde{C}=\frac{C_{+}+C_{-}}{2C_{bef}},\tilde{\rho }=\frac{C_{+}-C_{-}}{2C_{bef}},$$

$$\tilde{\phi }=\frac{ze\phi }{kT},\delta =\frac{{\lambda }_S}{{\displaystyle {\lambda }_D^{\prime }}},ϵ=\frac{2{\displaystyle {\lambda }_D^{\prime }}}{L},{\nu }^{\prime }=2C_{bef}{a}^3,\tilde{\phi }(0)=\frac{ze\phi (0)}{kT},$$

So, the Equations (40) and (41) take the form:

$$\frac{\partial \tilde{C}}{\partial \tilde{t}}=ϵ\frac{\partial }{\partial \tilde{y}}(\frac{\partial \tilde{C}}{\partial \tilde{y}}+\tilde{\rho }\frac{\partial \tilde{\phi }}{\partial \tilde{y}}+\frac{{\nu }^{\prime }\tilde{C}}{1-{\nu }^{\prime }\tilde{C}}\frac{\partial \tilde{C}}{\partial \tilde{y}})$$

(47)

$$\frac{\partial \tilde{\rho }}{\partial \tilde{t}}=ϵ\frac{\partial }{\partial \tilde{y}}(\frac{\partial \tilde{\rho }}{\partial \tilde{y}}+\tilde{C}\frac{\partial \tilde{\phi }}{\partial \tilde{y}}+\frac{{\nu }^{\prime }\tilde{\rho }}{1-{\nu }^{\prime }\tilde{C}}\frac{\partial \tilde{C}}{\partial \tilde{y}})$$

(48)

$$-{ϵ}^2\frac{{\partial }^2\tilde{\phi }}{\partial {\tilde{y}}^2}=\tilde{\rho }$$

(49)

The boundary conditions at $\tilde{y}=0.2$ become

$$\frac{\partial \tilde{C}}{\partial \tilde{y}}+\tilde{\rho }\frac{\partial \tilde{\phi }}{\partial \tilde{y}}+\frac{{\nu }^{\prime }\tilde{C}}{1-{\nu }^{\prime }\tilde{C}}\frac{\partial \tilde{C}}{\partial \tilde{y}}=0 \frac{\partial \tilde{\rho }}{\partial \tilde{y}}+\tilde{C}\frac{\partial \tilde{\phi }}{\partial \tilde{y}}+\frac{{\nu }^{\prime }\tilde{\rho }}{1-{\nu }^{\prime }\tilde{C}}\frac{\partial \tilde{C}}{\partial \tilde{y}}=0 \tilde{\phi }=\pm \tilde{\phi }(0)\pm \delta ϵ\frac{\partial \tilde{\phi }}{\partial \tilde{y}}$$

(50)

and the initial conditions for $t=0$ are as follows:

$$\tilde{C}(\tilde{y},0)=1,\tilde{\rho }(\tilde{y},0)=0 \tilde{\phi }=\tilde{\phi }(0)(1-\tilde{y})$$

(51)

From Equations (47)–(49) together with boundary (50) and initial conditions (51), we determine $\tilde{C}$, $\tilde{\rho }$ and $\tilde{\phi }$ as functions of time. The equations are solved numerically for different values of $C_{bef}$ and $\phi (0)$ in order to estimate the time required to reach the steady state. The reliability of the entire analysis is checked by reproducing the results of section III of Ref. [17] for the same parameter values. However, the numerical calculations are extremely time consuming for high concentrations, high potentials and duct widths, with singularities emerging in some cases. For this reason and in order to evaluate the model, we have selected two potentials, namely $\phi \left(0\right)=0.26\mathrm{V}$ and $\phi \left(0\right)=2.6\mathrm{V}$, and we have numerically solved the equations for a duct width of $L={10}^{-4}\mathrm{m}$.

To integrate the results with those from the preceding Section 4, we combine the definitions of the dimensionless quantities $\tilde{C}$ and $\tilde{\rho }$. This yields the following:

$$\frac{C_{+}}{C_{bef}}=\tilde{C}+\tilde{\rho }\hspace{1em}\hspace{1em}\hspace{1em}\frac{C_{-}}{C_{bef}}=\tilde{C}-\tilde{\rho }$$

(52)

In the case of $\phi \left(0\right)=0.26\mathrm{V}$, the concentration $C_{bef}$ is derived from Figure 8a, where $C_M={10}^{21}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$. Given that in Figure 8a, the final state $\frac{C_M}{C_{bef}}=0.72$ for $L={10}^{-4}\mathrm{m}$, it follows that $C_{bef}=1.385\times {10}^{21}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$. Figure 9a depicts the concentration $\frac{C_{-}}{C_{bef}}$ for $C_{bef}=1.385\times {10}^{21}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$, $\phi \left(0\right)=0.26\mathrm{V}$ as a function of the dimensionless duct width $\tilde{y}$ for 10 different times. The value of $\frac{C_M}{C_{bef}}$ that we are mainly interested in ($\frac{C_{-}}{C_{bef}}$ at the center of the duct where $\tilde{y}=1$) reaches the value 0.72 in dimensionless time $\tilde{t}=87$, which corresponds to $t=0.6\mathrm{s}$ (for our calculations, we considered that we have salt ions with $D=1.3\times {10}^{-9}{\mathrm{m}}^2/\mathrm{s}$, $z=1$ and duct width $L={10}^{-4}\mathrm{m}$).

Increasing the potential to $\phi \left(0\right)=2.6\mathrm{V}$ leads to a lower value of the final state $\frac{C_M}{C_{bef}}=0.0386$, as it is again derived from Figure 8a for $L={10}^{-4}\mathrm{m}$. Given that $C_M={10}^{21}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$, it follows that $C_{bef}=2.59\times {10}^{22}\raisebox{1ex}{$\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^3$}\right.$. Figure 10a depicts the concentration $\frac{C_{-}}{C_{bef}}$ as a function of the dimensionless duct width $\tilde{y}$ for 10 different times for $C_{bef}=2.59\times {10}^{22}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$. The ratio $\frac{C_M}{C_{bef}}$ ($\frac{C_{-}}{C_{bef}}$ at the center of the duct where $\tilde{y}=1$) reaches the value of 0.0386 in dimensionless time $\tilde{t}\cong 870$. This means that a reduction of ≈96% is possible for the ions at the central region of the duct for times in the order of $t=1.6\mathrm{s}$ (again, for our calculations, we considered that we have salt ions with $D=1.3\times {10}^{-9}{\mathrm{m}}^2/\mathrm{s}$, $z=1$ and duct width $L={10}^{-4}\mathrm{m}$). We note that when $\phi \left(0\right)$ is multiplied by ten, $\tilde{t}$ also increases tenfold, given that the amplitude L and concentration $C_M$ remain constant.

Figure 9b and Figure 10b show the value of $\frac{C_{-}}{C_{bef}}$ versus $\tilde{t}$ and $\tilde{y}$ in a 3-dimensional plot for $C_{bef}=1.385\times {10}^{21}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$, $\phi \left(0\right)=0.26\mathrm{V}$ and $C_{bef}=2.59\times {10}^{22}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$, $\phi \left(0\right)=2.6\mathrm{V}$, respectively. It is evident that at the higher potential (Figure 10b), the concentration rapidly increases at $\tilde{y}=0$, but it drops rapidly to zero in a region having a width that becomes wider with time (the zero-concentration width is approximately $0<\tilde{y}\le 0.2$ for small times). For the low potential, the concentration remains nonzero, as is clearly shown in Figure 9b. The explanation for this is that higher potential $\phi \left(0\right)$ exerts great force on neighboring charges almost instantaneously, resulting in emptying the space near the electrodes, which, of course, then receives charges from the main volume of the duct. Despite the mathematical problems, the method developed in Section 5 provides additional information to that obtained from Section 4 for the final state, not only because it determines time but also for two additional reasons. Firstly, a prerequisite for solving MPNP equations is the prior knowledge of the concentration $C_{bef}$ of the solution before the field is applied. In the method outlined in Section 4, it is necessary to predefine the final target concentration $C_M$ at the center of the duct. The second advantage is that the quantities used in MPNP equations are dimensionless. This also constrains the parameter of the duct width $L$, which is incorporated into the dimensionless parameters $\tilde{t}$ and $\tilde{y}$. As can be observed from Figure 9 and Figure 10, these are independent of $L$. From the definitions of $\tilde{t}$ and $\tilde{y}$, we see that when the duct width $L$ is doubled, the time t for the solution to reach the same concentration fraction also doubles, i.e., it follows the same curve across the width of the duct. Consequently, for $L={10}^{-3}\mathrm{m}$ (tenfold increase), the time t should be approximately ten times that calculated for $L={10}^{-4}\mathrm{m}$, i.e., it becomes 6 s and 16 s in the case of $\phi \left(0\right)=0.26\mathrm{V}$, $C_{bef}=1.385\times {10}^{21}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$ and ${\phi \left(0\right)=2.6\mathrm{V},C}_{bef}=2.59\times {10}^{22}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$, respectively.

6. Conclusions

This paper explores the model of ion movement under the influence of an electric field, taking into account steric effects due to ion size and compares the results with those of the Boltzmann distribution, where the ions are considered dimensionless. The general conclusion is that for relatively dilute water solutions (10¹⁹–10²⁴ ions/m³ final concentration at the center of the duct) with ions of valence $z=1$ (let us say saline water), the differences are minor for potentials below 1 V. However, for potentials above 1 V, the differences become significant, necessitating the use of the steric effect model for studying ion movement. Using the modified steric effect model, we demonstrated that for the initial concentrations, specifically $C_{bef}=2.59\times {10}^{22}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$ ($C_M={10}^{21}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$), we can achieve a 96% reduction in ion concentration using a 2.6 V potential for duct width 0.1 mm, which is achievable within a completion time of 1.6 s, while the percentage drops to 80% for duct width 1 mm ($C_M={10}^{21}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$). For smaller potentials, no noticeable decrease in concentration is observed, while for higher potentials, there is the case of other electrochemical phenomena taking place, so we must be very careful (at least theoretically we obtain impressive results). With the increase in the concentration of the ions, their reduction percentage in the main volume of the solution decreases but remains significant up to $C_M={10}^{23}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}/{\mathrm{m}}^3$ (final concentration at the center of the duct).

Another observation is that when we increased the potential $\phi (0)$ tenfold keeping all other parameters constant (such as the duct width L and the final concentration in the middle of the duct $C_M$), the dimensionless time $\tilde{t}$ also increased tenfold. This can be very useful, freeing us from repeated calculations.

Furthermore, the completion time has been found to be proportional to the duct width. Thus, for a duct width $L~1\mathrm{m}\mathrm{m}$, the completion time can be estimated to be around 16 s, given the same parameters. This insight allows us to determine the completion time for various duct widths, also removing the necessity for repeated numerical computations of the MPNP equations.