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

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


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][2][3][4][5][6][7][8][9][10][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 , 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. 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 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.

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 ⃗. Thus, ions are moving due to the combined motion of the water flow along the streamwise direction and along the yaxis due to the electric force with velocity ⃗ . The y-axis ionic flux of the positive ions is given by the relation = , where is the concentration of the positive ions. From the analysis that is presented in Appendix A and Equation (A5), we have where z is the number of overflow protons; = 1.6 × 10 Cb; ≡ (see Appendix 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. In the above Equation (1), the electric field is substituted by = − . Following the same procedure, the ionic flux of the negative ions is read as where z is the number of overflow electrons. Considering that both positive and negative ions have the same z, the same mobility = = and, thus, = = , the ionic fluxes can be read more simply as where the sign of the ionic fluxes indicates the direction of motion and z is considered positive.
The concentration conservation equation may read as where the concentration ± is in (mol/m ) and the ionic fluxes in mol/(m • s). The charge density due to the above ion concentration reads as and F = 96485.34 C/mol is the Faraday constant. By applying the 1D Poisson equation in the y-direction of the duct, it is found that The above Equations (3) and (5) constitute the so-called Poisson-Nernst-Planck (PNP) equations system, which is solved in the next Section.
It should be noticed that the present study is in accordance with the Debye-Huckel theory [47][48][49][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 = , where ≈ 80 is the relative permittivity of the water and = 8.85 × 10 F/m.

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 = , the electric potential is kept equal to zero, ( = ) = 0 , while (0) and − (0) are the potential at y = 0 and = , 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 where 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., ± = ∓ ± − ± = 0 for = 0, Finally, the initial conditions of the model are such that a uniform ion distribution is applied at t < 0, together with a zero potential.

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][55][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 Considering that C + C = 2 , where 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: it is found that where κ = Thus, the width κ 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: and κ = Thus, from the charge density, it is found that the concentrations of the positive and negative ions are varied as Validation of the Electric Circuit Models for Long Ducts In case of a duct of width of mm order or wider, ≥ 10 m, and since ~1 − 10 ̇, it is found for the rate ≤ 10 and, thus, may be considered almost equal to zero. Moreover, for such length of the duct and since ~10 m , it can also be considered that tanh κ ≈ 1. Thus, after these simplifications: 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 Equation (A10)). Since the fluid bulk is in quasi-neutral condition, its resistance per unit surface is equal to = ( ) (Equation (A7)), and for the time constant effective capacitor = (due to the charge confinement at the edges), we have 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 2 for the case of ions, for the salted water solution of ( = 1.3 • 10 m s ), as a function of the initial concentration 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.

Long Time Behavior
As the time increases, then (1 − ) → 1 and → 0, and Equations (14)- (21) have their equilibrium forms as The equilibrium solutions of Equations (25)- (29) can be easily verified by the recent results of Reference [46] for the present water solution with = 1.027 • 10 (in S.I.) (from Equation (19)) and = = 1. For = 0.0015, 0.015, 0.15 m, where tanh κ ≈ 1 and thus ( ℎ κ ≈ ℎ κ ) and for = 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 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 . This happens when C + C = 2 , 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 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 ) can be found by integrating Equation (25) of the linearized PNP solution as and its absolute value adjacent to each electrode is given for = as Based on the above approximations, → 0, the surface charge density is found to be exactly the same as the electric circuit RC model as

Exact Evaluation of the Electric Field and the Surface Charge Density
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: and by considering ( = ) = 0 and the electric field intensity is in the center of the duct, and by integrating it, it is found that 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, is almost diminished relative to the values of the second part of the undersquare quantity; thus, where is the potential at the outer Helmholtz plane (OHP) at distance . 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.
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, = = 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[ + (0) ≈ + (0) , something that can not be applied in general at the double layer. Now, the potential at the outer Helmholtz plane is  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, = = 1. As it can be observed, is scaled almost linear with (0) for values up to 0.1 V. However, as (0) is increased further, is saturated to values lower than 1 V independently of (0). Now, the potential distribution along y in the diffuse layer for ≤ ≤ 2 can be found from Equation (35) as Thus, the electric field intensity inside the compact layer is given by As it is observed from Figure 5 for the case of salty water, = = 1, the electric field intensity at the compact layer takes huge values independently of the ion concentration, especially as the potential increases.

Nonlinear PNP Solution
Starting from = ( − ) and setting Equation ( The above Equation (46) can be transformed into the linear Equation (12) subject to C + C = . 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 , i.e., C + C = 2 . This is valid throughout the duct width (for the case of salty water, = = 1) when φ < 1 or φ ≤ 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) The effect of the nonlinearity can be investigated when a third term in the Taylor series of the exponential function, = 1 + + ! + ! + ⋯ is considered. Then, for the additional term, it is where = [ ( ) ] , and = [ ( ) ] .
Thus, from Equation (46) and by setting = 0 for the steady-state, it is 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, = = 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 3 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 as it is also discussed in Reference [46].

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 = 1.686 gr/L = 29 mole/m to the final concentration of = 0.497 gr/L = 8.6 mole/m with V = 0.4 V; thus, the electric charge that it is removed per surface area is about | | = 2967 Cb/m . 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.

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 3 , 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 3 . 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.

Appendix A. General Theory of Ion Movement
It is supposed that a positive ion of concentration C+ 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 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 . The pseudoforce per ion due to this is given by the relation: 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 ⃗ , 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 ν): Quickly, equilibrium is established; thus, By using Equations (Α1)-(Α3), the velocity can be found as

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: Thus, the capacity per unit area is equal to = .
Since an additional capacitor is formed at the other wall of the duct, the total capacity per unit area is then 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 (Α4) when . The resistance per unit surface is given by = and, thus, The time constant of the RC circuit determines that characteristic time of the phenomenon as 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 ̇ (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., ~1 − 10 ̇ [51], while in the present calculations, it is assumed to be equal to ≈ 5 ̇. Thus, the capacity of the compact part per unit area is given by = .
b. The diffuse layer that is simulated by a Gouy-Chapman-type capacitor with effective width of the diffuse part = , and the known capacity per area unit is given by = . 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 The time constant of the RC circuit of the Stern model is then

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 In order to determine the A constant, the Laplace transformed Poisson Equation (5) is used, = − , and integrated once: where Γ is a constant to be determined by considering the current density, i = zeN (J − J ), and it vanishes at the duct walls. Using Equation (12), the current density can be found as