Mathematical Modeling of the Influence of Equilibrium Coefficient Variation on the Steady-State Transport of a Binary Electrolyte in the Cross-Section of a Desalination Channel

1. Introduction

Currently, there is an acute shortage of clean water, more than 40% of the world's population does not have access to clean water. Even more important is the problem of water scarcity in technological processes, as well as in agriculture [1,2,3]. The most effective method of solving this problem is electromembrane methods of water purification [4,5,6,7,8,9]. One of the promising directions is the use of electrodialysis machines, which are alternating concentration and desalination chambers separated by anion exchange and cation exchange membranes [10,11,12,13,14,15,16,17].

Let's consider the steady-state transport of binary salt ions, such as $NaCl$. The mathematical models in the diffusion layer and the one-dimensional cross-section of a desalination channel formed by an anion-exchange and cation-exchange membrane differ significantly, and the diffusion layer model cannot be considered a submodel of the channel cross-section, even though physically there are two diffusion layers in the channel cross-section, namely, near the anion-exchange and cation-exchange membranes. Indeed, let's consider the diffusion layer near one of the membranes, for example, a cation-exchange membrane. The equations describing the mathematical models of ion transport in the diffusion layer and the channel cross-section are identical, as are the boundary conditions at the imaginary solution/cation-exchange membrane boundary at $x=H$. The boundary conditions of these models differ at $x=0$, since the physical processes at this point differ. In the channel cross-section, the anion-exchange membrane is located at point $x=0$, while for the diffusion layer, point $x=0$ represents the depth of the solution where the concentrations are known and the electroneutrality condition is satisfied. Diffusion layer models assume that the solution is perfectly mixed at point $x=0$, so the concentration of cations and anions is the same and corresponds to the initial solution concentration. Thus, transport in the diffusion layer retains some features of the two-dimensional model. When considering a channel cross-section, the fundamental assumption is that the solution is immobile. This leads to fundamental differences in the mathematical models of ion transport in the channel cross-section and in the diffusion layer. For example, under the natural assumption of ideal selectivity of ion-exchange membranes, the anion flux through the cation-exchange membrane is zero, and all current is carried by cations. The magnitude of this current can be arbitrary, depending on the potential drop at points $x=0$ and $x=H$( Figure 1 ) In channel cross-section models, in addition to the anion current being zero, the cation current is also zero, since the anion-exchange membrane is also perfectly selective. Thus, the salt ion current in channel cross-section models is zero. Each model has its own advantages and disadvantages. For example, channel cross-section models allow one to consider the recombination/dissociation process in the extreme steady-state case, when the influence of salt ions is minimal, since the salt ion current is zero. Furthermore, the mathematical model of binary electrolyte transport in a channel cross-section provides a more detailed description of ion flows and interactions, allowing for complex geometric and physical conditions to be taken into account. This leads to increased prediction accuracy and process optimization compared to the diffusion layer model.

Water dissociation at interfaces and changes in the equilibrium coefficient are important for understanding the processes occurring during substance transfer through selective ion-exchange membranes. Changes in the rate coefficient of water molecule dissociation at interfaces (membrane/solution or membrane/membrane) can have several causes, as indicated in [18]. This paper examines the effect of an electric field on the equilibrium coefficient of the dissociation/recombination reaction of water molecules using a new mathematical model of steady-state salt ion transport in a channel cross-section, taking into account the dissociation/recombination reaction of water molecules.

In the works of Onsager [18], Timashev S.F., and Sheldeshov N.V. [19,20], the change in the equilibrium coefficient is associated with the magnitude of the electric field strength. In some studies [22,23], in contrast to the above-mentioned studies, the change in the equilibrium coefficient is associated not simply with the magnitude of the electric field strength, but with the magnitude of the space charge.

A change in the equilibrium coefficient, according to the proposed mathematical model, leads to a significant increase in the fluxes of $H^{+}$ and $O{H}^{-}$ ions. This leads to a local increase in the equilibrium coefficient, for example, in the quasi-equilibrium boundary layer region by a factor of 10–30, which corresponds in order of magnitude to experimental data [21].

2. A Mathematical Model of the Influence of the Dependence of the Dissociation Coefficient on the Magnitude of the Space Charge on the Stationary Transport of Salt Ions in the Channel Cross-Section

According to the concepts developed by Onsager, Timashev S.F., Sheldeshov N.V. et al. [18,19,20] for bipolar membranes, the current density of $H^{+}$($O{H}^{-}$) ions generated in the membrane system is determined by an empirical formula, where the equilibrium coefficient exponentially depends on the strength of the external electric field.

In electrodialysis systems with monopolar membranes, few studies have been devoted to the study of this dependence [22,23].

There are considerations from which it follows that, in contrast to the above-mentioned works, in systems with monopolar membranes, the change in the equilibrium coefficient can be associated not simply with the magnitude of the electric field strength, but with the magnitude of the space charge, that is, $k_w(\rho )=k_{w_0}+b\left|\rho \right|$, where $k_{w_0}={\displaystyle \frac{k_d}{k_r}}(0)C_{H_{2}O}={10}^{-8}$, $\rho =F(C_1(x)-C_2(x)+C_3(x)-C_4(x))$ the value of the space charge at point $x$, $b$ is an adjustable parameter, $k_d$ dissociation coefficient, $k_r$recombination coefficient, $C_i(x)$as the concentrations, of the i-th ion, $C_{H_{2}O}$ - water concentration, $F$is Faraday’s constant. The dependence of ion flows on the magnitude of space charge was first mentioned in the work of V.V. Nikonenko et al. [21]. To explain the influence of the magnitude of space charge on the equilibrium coefficient, let us consider the process of increasing dissociation of water molecules near the anion-exchange membrane. In this region, the space charge is determined by the predominance of the concentration of $C{l}^{-}$ ions over the concentrations of all other ions. As is known, a water molecule is a dipole, at one end of which are located $H^{+}$ and $O{H}^{-}$.$C{l}^{-}$ions surround the water molecule from the $H^{+}$side, weakening (shielding) and breaking its bond with $O{H}^{-}$. As a result, the resulting free $O{H}^{-}$ ions are transferred unhindered through the anion-exchange membrane and do not participate in further processes, while $H^{+}$ions move into the solution, and the higher the concentration of $C{l}^{-}$ions, the greater the flow of $H^{+}$ions. Thus, the greater the magnitude of the space charge, the more intense this process, which is in fact the electrolytic dissociation of water molecules, occurs; that is, the equilibrium coefficient kw depends on the magnitude of the space charge. This corresponds to the point of view of Onsager L., that the increase in the equilibrium coefficient is due to the gradual screening of the ion's electric field by the surrounding space charge [2]. In the region of the space charge in the vicinity of the cation exchange membrane (CEM), similar processes occur, only the role of $C{l}^{-}$ions is played by $N{a}^{+}$ or $K^{+}$ ions, $O{H}^{-}$ions are injected into the solution, and $H^{+}$ions are transferred through the CEM and do not participate in further transfer processes.

2.1. System of Equations

Under the above conditions, the steady-state transfer of salt ions for a 1:1 electrolyte into the channel cross-section, taking into account the space charge and the dissociation/recombination reaction, taking into account the dependence of the dissociation coefficient on the magnitude of the space charge, is described by the extended system of Nernst-Planck-Poisson equations [23].

Let us denote $\varphi (x,\epsilon )$ as the potential,$j_i(x,\epsilon )$, ${С}_i(x,\epsilon )$, $D_i$ as the flux, concentrations, and diffusion coefficients of the i-th ion, respectively. In the steady-state case $j_1=j_2=0$ [23]. Therefore, the current carried by the salt ions is equal to zero, and, consequently, the total current (formula) is equal to the current carried by the $H^{+}$ and $O{H}^{-}$ ions: $I_c=j_3-j_4$. Thus, the steady-state transport in the one-dimensional cross-section of the desalination channel can only be equilibrium with respect to the salt ions.

Thus, the system of equations has the form:

$$-{\displaystyle \frac{d{j}_i}{dx}}+R_i=0,i=1,\mathrm{...},4$$

(1)

$$j_i=-z_i{\displaystyle \frac{F}{RT}}{D}_i{C}_i{\displaystyle \frac{d\varphi }{dx}}-D_i{\displaystyle \frac{d{C}_i}{dx}},i=1,\mathrm{...},4$$

(2)

$${\displaystyle \frac{d^2\varphi }{d{x}^2}}=-{\displaystyle \frac{\rho }{{\epsilon }_r}}$$

(3)

$$R_3=R_4=k_d(\rho )C_{H_{2}O}-k_r{C}_3{C}_4=k_r(k_w(\rho )-C_3{C}_4)$$

(4)

$$z_1=1,z_2=-1,z_3=1,z_4=-1$$

Here, (1) represents the material balance equations; (2) are the Nernst–Planck equations for the fluxes of potassium or sodium ions ($i=1$), chloride ions ($C{l}^{-},i=2$), hydrogen ions ($H^{+},i=3$), and hydroxide ions ($O{H}^{-},i=4$); (3) is the Poisson equation for the electric potential; (4) are the formulas describing the dissociation/recombination reaction of water molecules; (5) is the current flow equation, which states that the current density $I$ through the cross-section of the desalination channel is determined by the ionic fluxes, i.e., it corresponds to the Faraday current density.

Here, ${\epsilon }_r$ is the dielectric permittivity of the solution, $R$ is the universal gas constant, $\varphi$ is the potential, $z_i$charge numbers, T - absolute temperature of the solution, $R_i$ - reaction speed, $E=-{\displaystyle \frac{d\varphi }{dx}}$ is the electric field strength, $k_r(0)=1,1\cdot {10}^8$ ${\displaystyle \frac{m^3}{s*mol}}$is the recombination rate constant, $k_d(0)=2\cdot {10}^{-5}$ is the dissociation rate constant of water molecules, and, consequently, $k_{w_0}={\displaystyle \frac{k_d}{k_r}}(0)C_{H_{2}O}={10}^{-8}$ ${\displaystyle \frac{mo{l}^2}{m^6}}$ is the equilibrium coefficient (ionic product of water).

2.2. Boundary Conditions

At $x=0$:

$$C_2(0)=C_{2a},$$

(5)

$${\left.\left(-{\displaystyle \frac{F}{RT}}{C}_1{D}_1{\displaystyle \frac{\partial \varphi }{\partial x}}-D_1{\displaystyle \frac{\partial C_1}{\partial x}}\right)\right|}_{х=0}=0,$$

(6)

$${\left.\left(-{\displaystyle \frac{F}{RT}}{C}_3{D}_3{\displaystyle \frac{\partial \varphi }{\partial x}}-D_3{\displaystyle \frac{\partial C_3}{\partial x}}\right)\right|}_{x=0}=j_{3a},$$

(7)

$${\displaystyle \frac{d{C}_4(0)}{dx}}=0,$$

(8)

$$\varphi (0)=d.$$

(9)

At $x=H$:

$$C_1(H)=C_{1k},$$

(10)

$${\left.\left({\displaystyle \frac{F}{RT}}{C}_2{D}_2{\displaystyle \frac{\partial \varphi }{\partial x}}-D_2{\displaystyle \frac{\partial C_2}{\partial x}}\right)\right|}_{x=H}=0,$$

(11)

$${\displaystyle \frac{d{C}_3(H)}{dx}}=0,$$

(12)

$${\left.\left({\displaystyle \frac{F}{RT}}{C}_4{D}_4{\displaystyle \frac{\partial \varphi }{\partial x}}-D_4{\displaystyle \frac{\partial C_4}{\partial x}}\right)\right|}_{x=H}=j_{4k},$$

(13)

$$\varphi (H)=0.$$

(14)

The concentrations of anions (5) and cations (10) at the membrane surfaces are determined by their exchange capacity, with the membranes assumed to be ideally selective (6), (11). It is assumed that in the surface layer of the membranes the generation of $H^{+}$ and $O{H}^{-}$ions occurs, and these ions freely pass through the cation-exchange and anion-exchange membranes (8), (12), no longer participating in the subsequent transport processes. At the same time, from the surface of the anion-exchange membrane a flux of $H^{+}$ ions is injected into the solution, determined by condition (7), while from the surface of the cation-exchange membrane a flux of $O{H}^{-}$ ions is injected (13). The membrane surfaces are considered equipotential, and the potential jump is defined by relations (9) (14).

3. Analysis of the Numerical Solution

The numerical solution was obtained using the finite element method. To analyze the results, calculations were performed with a wide range of initial data, the fitting parameter $b$, and so on. The most representative results are presented below: in each figure, the results for $b=0$ are shown on the left, while those for $b=0.0001$ are generally shown on the right.

To analyze the results of the numerical solution, we will use the following functions:

1. the charge distribution density normalized by Faraday’s constant ${\rho }_0={\displaystyle \frac{\rho }{F}}=C_1(x)-C_2(x)+C_3(x)-C_4(x)$,

2. the equilibrium function $P_w=k_{w_0}+bF\left|C_1-C_2+C_3-C_4\right|-C_3{C}_4$, or $P_w=k_{w_0}+b_0\left|C_1-C_2+C_3-C_4\right|-C_3{C}_4$, where $b_0=bF$,

3. the equilibrium coefficient $k_w(\rho )=k_{w_0}+b\left|\rho \right|$.

3.1. Concentration Distribution

As shown in Figure 2 (a, b), the qualitative behavior of the concentrations of all four ion species for $b=0$ and $b=0.0001$ is similar. Specifically, near the ion-exchange membranes, boundary layers are observed where the concentrations of chloride ions $C{l}^{-}$ and hydrogen ions $H^{+}$ at the anion-exchange membrane ($x=0$), as well as sodium ions $N{a}^{+}$ and hydroxide ions $O{H}^{-}$ at the cation-exchange membrane ($x=H$), vary at an exponential rate. Beyond these layers, there are relatively large regions adjacent to the boundaries where the concentrations of the above-mentioned ions change linearly or are nearly zero. Moreover, the concentrations of $C{l}^{-}$ and $H^{+}$ match with high accuracy, as do those of $N{a}^{+}$ and $O{H}^{-}$. In the central part of the channel, there exists a region with nearly constant concentration distributions, the values of which are practically independent of $b$.

Despite this qualitative similarity, there are several quantitative differences between the concentration distributions at $b=0$ and $b=0.0001$. The size and location of the central region differ at $b=0.0001$, it shifts to the right and becomes smaller. In addition, although the central region is narrower at $b=0.0001$, the linear concentration distributions are steeper, since the concentration value at the boundary of the adjacent layer is significantly higher.

From the analysis of the concentration distributions, it can be concluded that the cross-section of the desalination channel consists of at least seven regions: two boundary layers, two regions of linear concentration variation, and one region of nearly constant concentrations, as well as two intermediate layers between the regions of linear variation and the constant region. To fully understand the significance of these regions, it is necessary to consider additional characteristics of the transport processes.

Furthermore, unlike the diffusion layer (for example, near the cation-exchange membrane), where the steady state may be nonequilibrium — that is, the salt-ion current may not equal zero, and therefore the cation concentration may decrease almost to zero — in the channel cross-section, the steady state can only be equilibrium. In other words, the salt-ion current is zero, which means that salt ions cannot be completely removed. Consequently, during the transport of $H^{+}$ and $O{H}^{-}$ ions, the redistribution of salt ions occurs, but they are not carried beyond the channel boundaries.

3.2. Spatial Charge Distribution

In Figure 3, the spatial charge distribution shows that the channel cross-section can also be divided into seven intervals. The boundary layers represent regions of spatial charge associated with the properties of the interfacial boundary, that is, part of the conventional electric double layer [24], with the anion-exchange membrane carrying a negative charge and the cation-exchange membrane a positive charge. In the regions of linear concentration distribution, the condition of electroneutrality is satisfied with high accuracy. In the central part of the channel, however, an additional region of spatial charge appears (internal boundary layers), which is not related to the interfacial boundaries. To explain the emergence of these internal boundary layers, it is necessary to consider the function $P{}_W$ (the rate of water molecule dissociation/recombination reaction).

3.3. Graphs of the Equilibrium Function $P_w$

Figure 4 shows that in the regions of spatial charge, within the boundary layers near the ion-exchange membranes, intense water molecule dissociation occurs, as the electric field strength in these regions reaches very high values. Moreover, for$b=0.0001$, water dissociation is more intense than for $b=0$. In $[0,x{}_1)$, the ions $O{H}^{-}$ are completely transported through the anion-exchange membrane (boundary condition 19a), so their concentration is practically zero. At the same time, the concentration of ions $H^{+}$rapidly increases. The concentration of ions $C{l}^{-}$ decreases exponentially from $K^{+}$, while the concentration of (formula) remains practically zero.

In $({х}_4,H]$, the situation occurs in the opposite manner: the ions $H^{+}$ are completely transported through the cation-exchange membrane (boundary condition 12), so their concentration is practically zero. At the same time, the concentration of ions $O{H}^{-}$ increases rapidly. The concentration of ions $K^{+}$ grows exponentially toward $C{}_{1k}$, while the concentration of $C{l}^{-}$ remains practically zero.

The spatial charge concentration distribution in the central part of the channel, as shown above, significantly depends on the potential jump. Let us first consider the case when the potential is below the limiting value, for example, 0.3 V. In this case, a region $(x{}_2,x{}_4)$ forms in the central part of the channel where recombination dominates over dissociation; we will refer to this as the recombination region. The structure of the recombination region is almost symmetric with respect to the channel center for a $KCl$ solution, with a slight rightward shift, since the diffusion coefficients of potassium and chloride ions are nearly equal. However, for a $NaCl$ solution, a significant asymmetry is observed because the diffusion coefficient of sodium ions is smaller than that of chloride ions.

A characteristic feature of this region Is the presence of a deep, single local minimum of the equilibrium function, caused by the dominance of recombination over dissociation. The minimum point is defined by the location where the hydrogen ion concentration equals the hydroxide ion concentration, which is why recombination is maximal at this point. To the left of the minimum point, the hydrogen ion concentration exceeds that of hydroxide ions. Since in the recombination region the concentrations of potassium and chloride ions are much lower than those of hydrogen and hydroxide ions, an excess of hydrogen ions occurs to the left of the minimum point, creating a region of positive spatial charge. Similarly, to the right of the minimum point, a region of negative spatial charge forms due to an excess of hydroxide ions.

From this, it can be concluded that the formation of the internal boundary layer of spatial charge is associated with the finite recombination rate. Specifically, $H^{+}$ ions approach from the left and $O{H}^{-}$ ions from the right, and they recombine. Since the recombination rate is finite, albeit very high, not all ions recombine immediately, so they accumulate in a “queue” for recombination. This results in a positively charged spatial region on the left due to $H^{+}$ ions and a negatively charged spatial region on the right due to $O{H}^{-}$ions.

At a potential jump of ${\Delta }_r\varphi =1.5В$, a characteristic feature is the presence of two local minima in the equilibrium function (Figure 4), which is due to the very high electric field strength between these points. This causes dissociation to be so enhanced that near-equilibrium conditions are almost reached between the two local minima of the equilibrium function.

3.4. Graphs of the Fluxes $j_3$($H^{+}$) and $j{}_4$($O{H}^{-}$)

Figure 5 shows that accounting for the dependence of the equilibrium coefficient on the spatial charge leads to an increase in the fluxes: for $b=0.0001$, the fluxes rise by an order of magnitude in the boundary layers near the membranes and approximately double in the central part of the channel. By selecting an appropriate value of$b$, the flux can be increased by 10–30 times, as observed in experiments shown in Figure 5 and Figure 6 of [21]. Thus, choosing the dependence of the equilibrium coefficient on the spatial charge, along with a proper fitting parameter $b$, allows for both qualitatively and quantitatively accurate description of the enhancement of the fluxes of the products of water molecule dissociation.

In the electroneutral region, one of the fluxes is zero while the other proportional to the current I. In the equilibrium region, where the concentrations are practically constant, both fluxes take constant values, which are different from zero. This leads to the formation of two boundary layers for each flux in the vicinity of points $x_1$ and $x_4$.

3.5. Graphs of the Electric Field Intensity $E$

Figure 7 shows that the electric field intensity decreases numerically in the boundary layers with increasing parameter $b$. This is because, as noted above, the release of $H^{+}$ ions at the anion-exchange membrane and $O{H}^{-}$ ions near the cation-exchange membrane is enhanced. Since the spatial charges near the anion- and cation-exchange membranes are negative and positive, respectively, the $H^{+}$ and $O{H}^{-}$ ions reduce the magnitude of these spatial charges, leading to a corresponding decrease in the electric field intensity in these regions. To explain the increase in the electric field intensity in the middle of the channel, it should be noted that the concentrations in this region are constant and do not depend on $b$, while the fluxes increase. After a series of transformations, we obtain ${\displaystyle \frac{d(C_1-C_2+C_3-C_4)}{dx}}=(C_1+C_2+C_3+C_4)E-{\displaystyle \frac{1}{D_3}}{j}_3+{\displaystyle \frac{1}{D_4}}{j}_4$.

Since the concentrations are practically constant, ${\displaystyle \frac{d(C_1-C_2+C_3-C_4)}{dx}}\approx 0$, therefore

$$(C_1+C_2+C_3+C_4)E-{\displaystyle \frac{1}{D_3}}{j}_3+{\displaystyle \frac{1}{D_4}}{j}_4\approx 0$$

$$E\approx {\displaystyle \frac{1}{D_4{D}_3}}{\displaystyle \frac{D_4{j}_3-D_3{j}_4}{C_1+C_2+C_3+C_4}}$$

Thus, with an increase in b, the fluxes increase, the numerator grows while the denominator remains constant, and consequently, the electric field intensity increases proportionally with the fluxes.

3.6. Equilibrium Constant $K_W$

For $b=0$, the equilibrium constant $K{}_W$ is constant, whereas for $b=0.0001$it becomes variable and behaves similarly to the spatial charge, except that $K_W$ is positive and its maximum value depends on $b$, increasing as $b$ increases. Thus, for $b=0.0001$, the equilibrium constant increased by approximately 57 times. It can be shown that by selecting $b$, the maximum value of the equilibrium constant can be made arbitrarily large, which leads to an increase in the fluxes of $H^{+}$ and $O{H}^{-}$ ions. The increase in $H^{+}$ and $O{H}^{-}$ fluxes reduces the magnitude of the spatial charge near the anion- and cation-exchange membranes, since the signs of the spatial charge near the membranes and those of the $H^{+}$ and $O{H}^{-}$ ions are opposite. This, in turn, diminishes the effect of the spatial charge on enhancing the fluxes of $H^{+}$ and $O{H}^{-}$ ions, leading to a certain equilibrium. As a result, neither the electric field intensity nor the fluxes can increase indefinitely.

4. Structure of the Desalination Channel Cross-Section

From the analysis presented above, it is evident that the channel cross-section can be divided into seven intervals:

• Boundary layers near the anion- and cation-exchange membranes, which are quasi-equilibrium regions of spatial charge for $b=0$and their analogs for $b=0.0001$. In addition, there are two internal boundary layers of spatial charge, whose formation is associated with the recombination reaction of water molecules.
• Between the two internal boundary layers, there is an equilibrium region where the concentrations and electric field intensity are constant. The electric field intensity in this region is relatively high, which is associated with the increased fluxes of $H^{+}$and $O{H}^{-}$ions.
• In the remaining two regions, the conditions of equilibrium and local electroneutrality are satisfied, with the concentration distribution being nearly linear.

Using this structure of the desalination channel cross-section, it is possible to obtain an approximate analytical, for example, asymptotic solution to the boundary-value problem.

5. Conclusions

The study established that in systems with monopolar membranes, the change in the equilibrium constant can be linked not merely to the electric field strength, but to the magnitude of the spatial charge. According to the proposed mathematical model, variations in the equilibrium constant lead to a significant increase in the fluxes of $H^{+}$ and $O{H}^{-}$ions. In particular, the equilibrium constant locally increases, for example, in the quasi-equilibrium boundary layer region, by a factor of 10–30, which is in the same order of magnitude as the experimental data [21].

The work demonstrated that in the regions of spatial charge, within the boundary layers near the ion-exchange membranes, intense water molecule dissociation occurs, with higher values of the parameter $b$ leading to stronger dissociation. It was also shown that an internal boundary layer (the recombination region) forms deeper in the solution, associated with the recombination of $H^{+}$ and $O{H}^{-}$ ions. Furthermore, that with increasing b the flows increase and the electric field strength increases proportionally.