Theoretical Analysis of Electroconvection in the Electrodialysis Desalination Channel under the Action of Direct Current

The development of electroconvection in electromembrane systems is a factor that increases the efficiency of the electrolyte solution desalination process. The desalination of the solution, manifested by a change in the distribution of the ion concentration, can affect the mechanisms of development of electroconvection. The purpose of this work is to study the electroconvective flow developing in the desalination channel under various desalination scenarios. The study was carried out on the basis of a mathematical model of the transfer of binary electrolyte ions in the desalination channel formed between the anion and cation exchange membranes under the action of DC current. An analytical estimation of the threshold current density reflecting the conditions of the system transition into a quasi-stationary state has been obtained. The chronopotentiograms of the desalination channel and the thickness of the electroconvective mixing layer are calculated for both pre-threshold and supra-threshold values of the current density.


Introduction
Electrodialysis is actively used for the production of drinking water in the food and pharmaceutical industries [1,2]. The functioning of electrodialysis systems is determined by the selective properties of ion exchange membranes, which mainly pass ions of one sign (counterions) and prevent the movement of ions of another sign (coions) [3].
One of the methods of electrodialysis system research is chronopotentiometry, that is, studying the change in the potential jump in the system over time under the action of DC current. Based on the analysis of the results of chronopotentiometry, the mechanisms of mass transfer processes and the properties of ion exchange membranes are investigated [4][5][6][7][8]. A typical chronopotentiogram (ChP) of a single ion exchange membrane consists of the following sequence of characteristic sections [4][5][6][7][8]: 1.
An initial, almost vertical section, defined by the initial ohmic resistance of the solution; 2.
A region of slow growth of the potential jump, which is associated with a decrease in the concentration in the depleted solution near the membrane surface as a result of electrodiffusion processes; 3.
At overlimiting current densities, an inflection point of the ChP curve is noted, which is associated with the appearance of another ion transport mechanism [4]; thus, the third transition region of the ChP corresponds to the development of additional ion transport mechanisms; 4.
A section where the system is in a stationary or quasi-stationary state (when the potential jump oscillates about a certain fixed value).
Experimental studies [4,5] have shown that for dilute solutions of electrolytes, electroconvection is an additional mechanism responsible for the increase in the intensity of mass transfer and the appearance of an inflection point on the ChP. Electroconvection is understood as the movement of an electrolyte solution under the action of an electric field on a space charge [9][10][11][12]. The difference in the behavior of ions of opposite signs in an ion exchange membrane leads to the formation of a space charge region in a depleted solution near membrane surface, which extends with an increase in the flowing current density [13]. Therefore, at overlimiting current densities, electroconvection significantly affects the transfer of ions in membrane systems with dilute electrolyte solutions [11,[14][15][16][17][18][19].
The theoretical calculation of the ChP for a layer of a dilute electrolyte solution near the surface of a homogeneous ion exchange membrane, taking into account the development of electroconvection, which qualitatively coincides with the experimental data, was performed in [20]. On the calculated ChP, all the characteristic areas described above are distinguished (indicated by numbers in Figure 1). It is shown that the appearance and increase in the electroconvective mixing layer (the plot of this layer thickness in Figure 1 is denoted by d ec ) in the region near the membrane surface corresponds to the transition region of the ChP. Slipping of vortices and vortex complexes along the channel under the action of forced flow causes oscillations in the potential jump. In a quasi-stationary state (in which the potential jump oscillates about certain fixed value), the thickness of the electroconvective mixing layer also becomes quasi-stationary [20]. It was numerically shown that the development of electroconvection significantly reduces the potential jump in the quasi-stationary state of the electrolyte layer near the ion exchange membrane.  4. A section where the system is in a stationary or quasi-stationary state (when the potential jump oscillates about a certain fixed value).
Experimental studies [4,5] have shown that for dilute solutions of electrolytes, electroconvection is an additional mechanism responsible for the increase in the intensity of mass transfer and the appearance of an inflection point on the ChP. Electroconvection is understood as the movement of an electrolyte solution under the action of an electric field on a space charge [9][10][11][12]. The difference in the behavior of ions of opposite signs in an ion exchange membrane leads to the formation of a space charge region in a depleted solution near membrane surface, which extends with an increase in the flowing current density [13]. Therefore, at overlimiting current densities, electroconvection significantly affects the transfer of ions in membrane systems with dilute electrolyte solutions [11,[14][15][16][17][18][19].
The theoretical calculation of the ChP for a layer of a dilute electrolyte solution near the surface of a homogeneous ion exchange membrane, taking into account the development of electroconvection, which qualitatively coincides with the experimental data, was performed in [20]. On the calculated СhP, all the characteristic areas described above are distinguished (indicated by numbers in Figure 1). It is shown that the appearance and increase in the electroconvective mixing layer (the plot of this layer thickness in Figure 1 is denoted by ec d ) in the region near the membrane surface corresponds to the transition region of the ChP. Slipping of vortices and vortex complexes along the channel under the action of forced flow causes oscillations in the potential jump. In a quasi-stationary state (in which the potential jump oscillates about certain fixed value), the thickness of the electroconvective mixing layer also becomes quasi-stationary [20]. It was numerically shown that the development of electroconvection significantly reduces the potential jump in the quasi-stationary state of the electrolyte layer near the ion exchange membrane.  Since the ChP was calculated for the electrolyte layer near the surface of a single ion exchange membrane, the outer boundary of this layer is taken as the core of the flow with constant uniform distributions of ion concentrations, flow velocity, and electric potential [20]. This assumption causes the onset of a stationary or quasi-stationary state at any value of the current density, without taking into account the possibility of reducing the concentration in the core of the flow. In this work, numerical simulation of the mass transfer process in the desalination channel formed between the anion exchange (AEM) and cation Since the ChP was calculated for the electrolyte layer near the surface of a single ion exchange membrane, the outer boundary of this layer is taken as the core of the flow with constant uniform distributions of ion concentrations, flow velocity, and electric potential [20]. This assumption causes the onset of a stationary or quasi-stationary state at any value of the current density, without taking into account the possibility of reducing the concentration in the core of the flow. In this work, numerical simulation of the mass transfer process in the desalination channel formed between the anion exchange (AEM) and cation exchange (CEM) membranes is carried out. The conditions for the establishment of a quasi-stationary state in the desalination channel under the action of DC current are considered, and an analytical estimate of the threshold current density for the transition of the system to the quasi-stationary state is obtained. The ChPs of the desalination channel and the thickness of the electroconvective mixing layer are calculated for both pre-threshold and, for the first time, supra-threshold values of the current density.

Mathematical Model
Consider a binary electrolyte solution in an electrodialysis desalting channel formed between two ion exchange membranes (AEM and CEM), which is pumped at the average velocity V 0 ( Figure 2). Let x and y be the normal and tangent coordinates to the membrane surfaces, respectively; x = 0 and x = H correspond to the interfaces of the electrolyte solution with AEM and CEM; and y = 0 and y = L correspond to the inlet and outlet from the channel ( Figure 2).
Membranes 2022, 12, x FOR PEER REVIEW 3 of 15 exchange (CEM) membranes is carried out. The conditions for the establishment of a quasi-stationary state in the desalination channel under the action of DC current are considered, and an analytical estimate of the threshold current density for the transition of the system to the quasi-stationary state is obtained. The ChPs of the desalination channel and the thickness of the electroconvective mixing layer are calculated for both pre-threshold and, for the first time, supra-threshold values of the current density.

Mathematical Model
Consider a binary electrolyte solution in an electrodialysis desalting channel formed between two ion exchange membranes (AEM and CEM), which is pumped at the average velocity 0 V ( Figure 2). Let x and y be the normal and tangent coordinates to the membrane surfaces, respectively;  The non-stationary process of transfer of binary electrolyte ions in membrane systems in the absence of chemical reactions taking into account electroconvection is described by the following equations: The non-stationary process of transfer of binary electrolyte ions in membrane systems in the absence of chemical reactions taking into account electroconvection is described by the following equations: → j n = −z n D n c n ∇ϕ − D n ∇c n + Pe c n Equations (1)-(4) are given in a dimensionless form: time, t, is scaled by the value H/V 0 ; spatial coordinates, x and y, by the channel width H; velocity, → V, by the average velocity of the forced flow V 0 ; pressure, p, by the velocity head ρV 2 0 ; the concentration of the n-th ion, c n , by the input concentration of the electrolyte c 0 ; electric potential, ϕ, by the value RT/F; the diffusion coefficient of the n-th ion, D n , by the electrolyte diffusion ; and the flux of ions of the n-th kind, → j n , by the value Dc 0 /H. Here, are the dimensionless parameters [21]; Z n is the charge number of the n-th ion; F is the Faraday constant; R is the gas constant; T is the absolute temperature; ε 0 is the electrical constant; ε r is the relative permittivity of the solution; ρ 0 is the density of the solution; ν is the kinematic viscosity; and ε r , ρ 0 and ν are assumed to be constant.  [22]: When modeling the DC current mode, the conditions of constancy of the average current density at the solution/membrane interfaces, i av , must be satisfied [23]: The vector field of the total current density is solenoidal, div → i = 0, so the electric current function, η, can be defined as follows [23]: On the basis of Equations (2), (5) and (7), Equation (8) was obtained, which makes it possible to determine the current density distribution in the DC current mode [20]: The numerical solution of the boundary value problem of the model requires the definition of the boundary conditions for Equations (1), (3), (4) and (8). Let us set boundary conditions similar to the conditions in [20] with the difference that the entire desalination channel from AEM to CEM is considered.
At the channel inlet, x ∈ [0, 1], y = 0, the parabolic Poiseuille velocity profile is taken, Equation (9); for ion concentrations, the Dankwerts conditions are taken, which determine the fact that the rate of entry of ions into the channel is equal to the rate at which ions cross the plane y = 0 through a combination of electromigration, diffusion, and convection, Equation (10); and the condition for the electric potential obtained from Equations (2) and (6) taking into account the zero density of the tangential current, i y (x, 0, t) = 0 (the tangential component of the displacement current is negligibly small), Equation (11): −z n D n c n ∂ϕ ∂y − D n ∂c n ∂y At the channel outlet, x ∈ [0, 1], y = L, the "soft" boundary conditions are accepted [24], Equation (12); tangential derivatives of ion concentrations and potential are set equal to zero, Equations (13) and (14): At the solution/CEM interface, x = 1, y ∈ [0, L], the no-slip condition is applied, Equation (15); the counterion concentration, c 1 , is set as a constant value that is N c times greater than the input concentration of solution [13], Equation (16); continuous flux of coions, where the transfer numbers of ions in membranes, T nA and T nC , determine the portion of the current carried by ions of this type (moreover, T 1A + T 2A = 1 and T 1C + T 2C = 1) [25], Equation (17); the normal derivative of the electric field potential is given as a function of the electric current density (i x = ∂η/∂y), Equation (18) [26]: The similar conditions are accepted for the velocity and concentrations of ions at the solution/AEM interface, x = 0, y ∈ [0, L]: In the system of Equations (1)-(4), (8) the electric potential enters only in the form of derivatives on spatial coordinates; therefore, only the potential jump is significant. For the convenience of calculations, we set the zero potential on the left boundary, Equation (22): Boundary conditions for the electrical current function [23]: At the initial time, t = 0, the flow is parabolic, Equation (24); the electrical neutrality condition is satisfied at all points of the desalination channel, and the ion concentrations are equal to the initial electrolyte concentration c 0 , Equation (25); the potential is zero everywhere, Equation (26): ϕ(x, y, 0) = 0.

The Condition for the Transition of the System to a Quasi-Stationary State
To determine the condition for the transition of the electrolyte solution in the desalination channel to a quasi-stationary state, we perform a sequence of transformations of the equation for the cation concentration, Equation (3): We multiply Equation (27) by the charge number of the cation, z 1 , and perform the integration over the considered region of the desalination channel: The first term on the right side of Equation (28) is determined taking into account the boundary conditions (17) and (21), and the fact that T 1C + T 2C = 1: In the last transformation (29), the conditions for the constancy of the average current density at the solution/membrane interfaces, Equation (6), are taken into account.
The second term on the right side of Equation (27) is determined taking into account the boundary conditions (9)-(14): Then, Equation (27) can be written as: The terms on the right side of Equation (31) are the values of the partial current of cations through the boundaries of the system under consideration. The condition for the transition of the system to a stationary (or quasi-stationary) state is the equality to zero (or smallness) of the right side of Equation (31): where i av is the current density and c 1 (x, L) is the concentration distribution at the outlet of the channel in the stationary state (or their time-averaged values in the quasi-stationary The zero flux of cations at the outlet of the channel determines the balance state of the processes of desalination of the electrolyte solution and its entry into the channel through forced flow. The current density for these conditions is determined by Equation (34): or in dimensional form (here L d is the dimensional length of the channel): Therefore, Equation (34) defines the threshold current density for the transition of the system to the quasi-stationary state, i tr : at value of current density less than i tr the system goes into the stationary (or quasi-stationary at overlimiting currents) state, and when the value of current density is exceeded i tr , the electrolyte is completely depleted, and the potential jump increases rapidly.
Note that Equation (32) makes it possible to estimate the current efficiency (CE) and the energy per ion removal (EPIR) in the DC current mode. CE shows how efficiently the current is used to remove ions [27]: EPIR is determined by Equation (37) (scaled by RT) [27]: For a quasi-stationary state, taking into account Equation (32), it can be written: where ∆ϕ is the potential jump in the stationary state or the time-averaged value in the quasi-stationary state. EPIR in dimensional form (J/mol): Thus, in the quasi-stationary state of desalination by DC current CE < 1 due to imperfect membrane selectivity, EPIR is determined by the potential jump and ion transport numbers in membranes.
When analyzing the results of studying processes in membrane systems, it is customary to normalize the current density by the value of the diffusion limiting current density, i lim . For the solution under consideration (NaCl) we use the limiting current density of the CEM defined by Equation (41) [29]: where t 1 is the transport number of the cation in solution, t 1 = 0.395 [30]. The limiting current density, i lim , determines the state of almost complete depletion of the ion concentration at the solution/CEM interface. When the limiting current density i lim is exceeded, an extended space charge region is formed and an electroconvective flow develops in the region near the CEM surface [20,28]. When the current density is higher Experimental studies of electrodialysis systems, as a rule, are carried out under conditions when the threshold current density of the quasi-stationary state, i tr , is greater than the values of the limiting current density for CEM and AEM, that is For example, the system settings described above at channel length L = 8 correspond to these conditions. Then, the threshold current density calculated by Equation (34) will be equal to i tr ≈ 1.76i lim . Figure 3 shows ChPs calculated for the following values of the current density i av /i lim : 1, 1.1, 1.5, 1.6, 1.7, 1.8, 1.9.      At the limiting current density, i av = i lim , at the ChP of the system under consideration, there is a successive change in the initial ohmic region by a region of slow growth of the potential jump (associated with a decrease in concentration in the regions near the surfaces of both membranes) and a transition to a stationary state without the development of electroconvection (Figure 4).
At values of the current density above the limiting value, i lim , electroconvection develops in the region near the CEM, and an electroconvective mixing layer is formed (see Figure 4). In order to quantitatively characterize the intensity of electroconvection, the averaged over the channel length thickness of the electroconvective mixing layer, d ec , is calculated. The boundary of the electroconvective mixing layer was determined as a point at which the difference in the root-mean-square value of the velocity in calculations with and without taking into account electroconvection exceeds 15% of the average forced flow velocity, V 0 , (by analogy with [31]). In the considered case of a desalination channel with two membranes, the thickness of the electroconvective mixing layer is calculated for the channel halves at the CEM (d ec CEM ) and at the AEM (d ec AEM ), Figure 3c,d.
On the ChPs for the current density equal to 1.1i lim and 1.5i lim , at the moment of development of the electroconvective mixing layer near the surface of CEM, a decrease in the growth rate of the potential jump and a transition to a quasi-stationary state, accompanied by an increase in d ec CEM , are observed. The shapes of these curves differ significantly in the transition region, since as the current density approaches the value of the limiting state of the solution at AEM, i lim (T 1C − t 1 )/(T 2A − t 2 ), the effect of desalination in this region increases. The ChPs for the current density, which exceeds the limiting value for the CEM, but are less than the limiting current value for the AEM, differs from the ChPs for the underlimiting and limiting densities by the appearance of a transition section, which corresponds to an increase in the thickness of the electroconvective mixing layer before the transition to the quasi-stationary state. Thus, the ChPs for the current density in the range i lim < i av < i lim (T 1C − t 1 )/(T 2A − t 2 ) contain all the regions characteristic of a single membrane, observed experimentally [4][5][6][7][8]. At values of the current density above the limiting value, lim i , electroconvection develops in the region near the CEM, and an electroconvective mixing layer is formed (see Figure 4). In order to quantitatively characterize the intensity of electroconvection, the averaged over the channel length thickness of the electroconvective mixing layer, dec, is calculated. The boundary of the electroconvective mixing layer was determined as a point at which the difference in the root-mean-square value of the velocity in calculations with and without taking into account electroconvection exceeds 15% of the average forced flow velocity, V0, (by analogy with [31]). In the considered case of a desalination channel with two membranes, the thickness of the electroconvective mixing layer is calculated for the channel halves at the CEM (dec CEM) and at the AEM (dec AEM), Figure 3c , the effect of desalination in this region increases. The ChPs for the current density, which exceeds the limiting value for the CEM, but are less than the limiting current value for the AEM, differs from the On the curve of the thickness d ec CEM for the current density of 1.1i lim , there are periodic oscillations of small amplitude (≈0.001 in the quasi-stationary state), which are associated with the movement of single electroconvective vortices (Figure 4) along the membrane surface from the inlet to the outlet of the desalination channel under the action of forced electrolyte flow. Thickness d ec CEM fluctuations for the current density of 1.5i lim have a more complex shape and a larger amplitude (≈0.012 in the quasi-stationary state), which is due to the fact that instead of single vortices, larger vortex complexes are observed (Figure 4).
At current densities of 1.6i lim , 1.7i lim , 1.8i lim , and 1.9i lim (which exceed the value of the limiting state density at AEM, i lim (T 1C − t 1 )/(T 2A − t 2 )), an electroconvective flow also develops in the region near the surface of the AEM (Figures 3d and 4), while at the ChPs, there is a second slowdown in the growth of the potential jump. Thus, ChPs for the current density, which exceeds the limiting value for the AEM, i lim (T 1C − t 1 )/(T 2A − t 2 ), but is less than the threshold value, i tr , differ from the one described above in that the transition state in which the electroconvective flow develops is divided into two fragments, one of which is associated with the appearance of electroconvective vortices in the region near the CEM, and the second one, in the region near the surface of the AEM. With an increase in current density, an increase in the size of electroconvective vortex complexes is fixed, accompanied by an increase in the amplitude of electroconvective mixing layers' thickness fluctuations at the CEM (d ec CEM ) and at the AEM (d ec AEM ).
At a current density equal to 1.8i lim and 1.9i lim (which exceed the threshold value, i tr ), all the characteristic sections of the ChPs described above are observed, except for the quasi-stationary state. Desalination of the solution occurs until the concentration is completely depleted (Figure 3b), while an increase in the potential jump is observed.
Changes in the structure of the electroconvective flow in the desalination channel for the current density 1.9i lim can be seen in Figure 5. The appearance of electroconvection near the CEM occurs almost simultaneously along the entire membrane surface, which is associated with the relative homogeneity (in the longitudinal direction) of solution depletion in this region (Figure 5a). The appearance of electroconvection near the AEM surface occurs later at a higher degree of solution desalination (Figure 3b). The ion concentration field at this point in time (Figure 5c) is characterized by large concentration field gradients in the region at the channel inlet compared to the region at the channel outlet. Therefore, for the considered parameters of the system, electroconvection near AEM is developed earlier in the region near the channel inlet. This explains the rapid growth of d ec CEM and the more gentle growth of d ec AEM immediately after the development of electroconvection. The depletion of the electrolyte solution in the region near the channel outlet is accompanied by the destruction of the space charge regions and, consequently, by a decrease in electroconvection (Figure 5f) [32,33]. On all ChPs, for a current density above the limiting value, potential jump oscillations are observed, which appear almost simultaneously with the development of electroconvection (Figure 3a,c,d). This is due to the non-stationary dynamics of vortices (slip of On all ChPs, for a current density above the limiting value, potential jump oscillations are observed, which appear almost simultaneously with the development of electroconvection (Figure 3a,c,d). This is due to the non-stationary dynamics of vortices (slip of vortices in the tangential direction or changes in the structure of the vortex complex). The non-stationary dynamics of vortices cause fluctuations in the ion concentration, and, consequently, the space charge density in the region of the electroconvective mixing layer, which determines the fluctuations in the potential jump. Similar fluctuations in the potential jump and the size of electroconvective vortices were recorded in a direct experimental study of electroconvection near the surface of an ion exchange membrane under the action of DC current in [34]. Figure 6 shows the results of calculations of the EPIR (by Equation (39)), performed with and without taking into account the development of electroconvection for the following values of current density i av /i lim : 0, 0.1, . . . , 1.7. Figure 6 shows that EPIR increases with increasing current density. Exceeding the limiting density in the calculation without taking into account electroconvection is accompanied by a sharp increase in EPIR. The appearance of electroconvection first at the CEM (at i lim ) and later at the AEM (see the value at 1.6i lim ) significantly reduces EPIR. Therefore, the development of electroconvection significantly reduces the EPIR. At the supra-threshold current density, the sharp increase in the potential jump causes large EPIR.

Conclusions
The study of the desalination process of the electrolyte solution in the electrodialysis desalination channel under the action of DC current by means of numerical simulation is carried out. The model makes it possible to describe the development of electroconvection caused by the action of an electric field on the space charge formed near the surface of ion exchange membranes, taking into account the process of desalting the electrolyte solution. The analytical estimate of the threshold current density of the transition to a quasi-stationary state is derived. ChPs calculated for pre-threshold current densities are characterized by a transition to a quasi-stationary state over time. The thickness of the electroconvective mixing layer under these conditions is also a quasi-stationary quantity. The ChP for suprathreshold current densities reflects an unlimited growth of the potential jump with a gradual depletion of the ion concentration. In supra-threshold current modes, electroconvection covers the entire channel over time. The depletion of the electrolyte concentration causes the destruction of the space charge regions and the attenuation of electroconvection. It is shown that in pre-threshold current regimes, the development of electroconvection significantly reduces the energy spent on ions' removing.

Conclusions
The study of the desalination process of the electrolyte solution in the electrodialysis desalination channel under the action of DC current by means of numerical simulation is carried out. The model makes it possible to describe the development of electroconvection caused by the action of an electric field on the space charge formed near the surface of ion exchange membranes, taking into account the process of desalting the electrolyte solution. The analytical estimate of the threshold current density of the transition to a quasi-stationary state is derived. ChPs calculated for pre-threshold current densities are characterized by a transition to a quasi-stationary state over time. The thickness of the electroconvective mixing layer under these conditions is also a quasi-stationary quantity. The ChP for supra-threshold current densities reflects an unlimited growth of the potential jump with a gradual depletion of the ion concentration. In supra-threshold current modes, electroconvection covers the entire channel over time. The depletion of the electrolyte concentration causes the destruction of the space charge regions and the attenuation of electroconvection. It is shown that in pre-threshold current regimes, the development of electroconvection significantly reduces the energy spent on ions' removing. Funding: The reported study was funded by RFBR and DFG according to the research project № 20-58-12018 NNIO (RFBR).