Mathematical Modeling of Monovalent Permselectivity of a Bilayer Ion-Exchange Membrane as a Function of Current Density

Modification of an ion-exchange membrane with a thin layer, the charge of which is opposite to the charge of the substrate membrane, has proven to be an effective approach to obtaining a composite membrane with permselectivity towards monovalent ions. However, the mechanism of permselectivity is not clear enough. We report a 1D model based on the Nernst–Planck–Poisson equation system. Unlike other similar models, we introduce activity coefficients, which change when passing from one layer of the membrane to another. This makes it possible to accurately take into account the fact that the substrate membranes usually selectively sorb multiply charged counterions. We show that the main cause for the change in the permselectivity coefficient, P1/2, with increasing current density, j, is the change in the membrane/solution layer, which controls the fluxes of the competing mono- and divalent ions. At low current densities, counterion fluxes are controlled by transfer through the substrate membrane, which causes selective divalent ion transfer. When the current increases, the kinetic control goes first to the modification layer (which leads to the predominant transfer of monovalent ions) and then, at currents close to the limiting current, to the depleted diffusion layer (which results in a complete loss of the permselectivity). Thus, the dependence P1/2 − j passes through a maximum. An analytical solution is obtained for approximate assessment of the maximum value of P1/2 and the corresponding fluxes of the competing ions. The maximum P1/2 values, plotted as a function of the Na+ ion current density at which this maximum is reached, gives the theoretical trade-off curve between the membrane permselectivity and permeability of the bilayer monovalent selective ion-exchange membrane under consideration.


Introduction
To date, membrane methods of purification, separation and concentration are the most environmentally friendly and economically promising. One of these methods is electrodialysis (ED), the long-term use of which on an industrial scale confirms its effectiveness and expediency [1][2][3].
Conventional electrodialysis desalination, which involves removing any ions and thereby reducing the total salinity of a solution [1], is in demand, e.g., for wastewater treatment, the production of deionized water [4,5] and water for the food industry [6]. However, there are also applications for electrodialysis in which the removal of certain kinds of ions is important. For example, NaCl and KCl are removed from whey in the dairy industry [7], while calcium and organic ions are valuable components and their removal is undesirable. In upgrading groundwater for irrigation use, monovalent Na + and Cl − ions are removed, and multivalent hardness cations and sulfate anions, which are necessary for optimal plant growth, are retained [8]. Selectrodialysis is a novel kind of electrodialysis where c 0 i and C 0 i = |z i |c 0 i are the molar (mol·m −3 ) and equivalent (eq·m −3 ) ion concentrations in the bulk solution, respectively.
The problem of the dependence of P 1/2 on the electric current density, j, was vigorously discussed in the 1960s and 1970s in connection with the development of a technology for obtaining table salt from sea water by electrodialysis. A great contribution to the study of the P 1/2 − j dependence was made by the authors of refs. [24][25][26][27][28]. At present, this issue is once again relevant due to the rapid development of technologies for the production and application of IEMs that are selective to singly charged ions [18,19,[29][30][31][32][33]. The general principle for the manufacture of such membranes is the formation of a thin active surface layer, which serves as a barrier to the transfer of ions that are counterions to the substrate membrane. This barrier creates only insignificant resistance for singly charged ions but is a serious obstacle for two-and especially for three-charged ions. The ideology of creating this barrier goes back to the works of Sata [34]. Such a barrier can be formed by applying a hydrophobic film [35], as well as by applying a single layer [17,22,30,[36][37][38] having fixed ions with a charge opposite to that of the substrate membrane, or by several surface layers (layer-by-layer) [17,[39][40][41] with alternating charge signs of fixed ions. The best results are obtained when applying several binary layers: the permselectivity coefficient increases with an increase in the number of such layers, n, but when n > 10, the increase in P 1/2 is no longer observed [18,42].
For monopolar IEMs, there is a continuous loss of selectivity transfer of the preferentially transferred ion with increasing current density [27,28,36,43,44]. This effect is due to a change in the control of the counterion transfer kinetics that occurs with increasing current density. At low currents, the fluxes of competing counterions are controlled by the membrane: usually divalent ions preferentially pass through the membrane due to their greater sorption by the membrane phase [45]. As the current increases, the membrane still controls the counterion-coion permselectivity, but the kinetic control of the permselectivity between two competing counterions is increasingly shifted to the depleted diffusion layer, in which there is no selective transport of competing ions [28,46,47]. However, in the case of an IEM modified with a surface layer with a charge that is opposite to that of the substrate membrane, the P 1/2 value firstly increases with j, then the curve passes through a maximum, after which a further increase in current density leads to a drop in specific permselectivity [22,36].
Roghmans et al. [36] investigated the possible reasons for the decrease in the specific permselectivity of a commercial CMX cation-exchange membrane modified with a thin anion-exchange layer by increasing current density when treating a ternary electrolyte solution. Based on the results of numerical simulations, they found that the presence of defects in the modification layer (domains not coated with selective material) cannot be the cause of the experimentally observed reduced permselectivity at high currents. The authors [36] suggested that structural changes in the modification layer, such as swelling, can cause a decrease in the fixed ions' charge density in the layer, thereby reducing ion repulsion and permselectivity.
Golubenko et al. [22] studied the permselectivity of an anion-exchange membrane (AEM) modified with a thin cation-exchange layer. They agree with the hypothesis made by Roghmans et al. [36] and offered a possible specific mechanism for reducing the concentration of fixed groups: when water splitting begins, some of the fixed SO 3 − ions in the cation-exchange modification layer can transform into the protonated uncharged form, which reduces the density of the fixed charge. Golubenko et al. [22] also proposed another mechanism that can take place in the case of the transfer of SO 4 2− ions, which at a high concentration of H + ions in the cation-exchange layer can transform into the singly charged form HSO 4 − . In our study, we developed a new one-dimensional, stationary Nernst-Planck-Poisson model that disregards water splitting, and on the basis of which a different (compared to Refs. [22,36]) explanation of the extremal form of the P 1/2 − j dependence for a bilayer membrane is given. We show that, at low current densities, counterion fluxes are controlled by transfer through the substrate membrane, as in the case of monopolar single-layer membranes. As the current increases, the kinetic control passes first to the modification layer (providing the selective transfer of monovalent ions), and then to the depleted diffusion layer. The latter, as in the case of monopolar membranes, leads to a complete loss of permselectivity. It follows from the simulation that the P 1/2 − j dependence also passes through a maximum, even when there is neither a decrease in the concentration of fixed ions in the modification layer, nor any water splitting.

Mathematical Model
The system under study consists of four layers: substrate CEM, thin anion-exchange modification layer, and two diffusion layers (DLs) adjacent to its surfaces. Two identical bulk solutions containing two kinds of cations (Na + and Ca 2+ , for definitiveness) and one kind of anion (Cl − ) flank the system under study (Figure 1).
The main assumptions in the formulation of this problem were similar to those in reference [44]:

•
The substrate and modification layers are considered as a homogeneous medium in which fixed charged groups are uniformly distributed; • It is assumed that the solvent flux through the membrane is negligible, and therefore the phenomena of osmosis and electroosmosis are not taken into account; • The impact of convection transport in solution is taken into account implicitly through the diffusion layer thickness, which is considered independent of the voltage applied; • The gradients of temperature, pressure and solution density are ignored; • Water splitting and electroconvection are not taken into account.
The ion-exchange material of both layers can be represented as a charged sponge impregnated with a charged solution of the opposite sign. The membrane system can also be imagined as a bundle of nanometer pores with charged walls. One such pore is shown in Figure 1b: the wall charge on one edge of the pore is positive (anion-exchange modification layer) and on the other side it is negative (cation-exchange layer). Our model considers 1D ion transport along an axis normal to the membrane surface. It is possible to arrange such a distribution of the charged sites on the pore walls such that the average concentration of them in the cross section normal to the transport axis is the same as in the charged sponge representation. Thus, both images of the bilayer membrane, based on a charged sponge or a bundle of nanopores, can serve to visualize the 1D mathematical model developed in this paper. The main assumptions in the formulation of this problem were similar to those in reference [44]:

•
The substrate and modification layers are considered as a homogeneous medium in which fixed charged groups are uniformly distributed; • It is assumed that the solvent flux through the membrane is negligible, and therefore the phenomena of osmosis and electroosmosis are not taken into account; • The impact of convection transport in solution is taken into account implicitly through the diffusion layer thickness, which is considered independent of the voltage applied; • The gradients of temperature, pressure and solution density are ignored; • Water splitting and electroconvection are not taken into account.
The ion-exchange material of both layers can be represented as a charged sponge impregnated with a charged solution of the opposite sign. The membrane system can also be imagined as a bundle of nanometer pores with charged walls. One such pore is shown The stationary transport of ions in the system under study is described by the Nernst-Planck-Poisson equations and the material balance equation: here J i , c i , D i , z i , and γ i are the flux density, molar concentration, diffusion coefficient, charge number and activity coefficient of ion i, where i = 1 and i = 2 corresponds monovalent (Na + ) и divalent (Ca 2+ ) counterions, respectively, and i = 3 corresponds to Cl − coion; R is the gas constant; T is the temperature; F is the Faraday constant; E = − dϕ dx is the electric field strength; ϕ is the electric potential; ε 0 is the vacuum permittivity; ε is the solution relative permittivity; z and Q are the charge number and concentration of fixed ion groups, respectively.
Equations (2)-(4) are valid for the DLs, substrate membrane, and modification layer. However, parameters z and Q depend on the coordinate: z and Q is set to zero in the DLs; z = z ML , Q = Q ML in the modification layer; z = z m , Q = Q m in the substrate membrane. The ion diffusion coefficients, D i , change in a similar way: D i = D s i in the DLs, D i = D ML i in the modification layer, D i = D m i in the substrate membrane. D m i are found from the experimentally determined value of electrical conductivity; since the electrical conductivity of the modification layer cannot be measured, we set the D ML i value τ times less than in solution, τ = 3 for all ions and cases considered. This parameter may be interpreted as a tortuosity factor. The γ i values are taken equal to unity in the DLs and modification layers, while in the substrate membrane γ i = γ m i , where γ m i can differ from 1. For a smoother change in these parameters at the boundaries of the respective layers, the weight boxcar function (rectangular wave) is used, while the thickness of all three interfacial transition regions is chosen to be 1 nm, which is close to the dense part of the electrical double layer [48]. A similar function was used by Evdochenko et al. [39] when modeling the charge density distribution of a polyelectrolyte multilayer membrane.
The current density, j, in the system includes the charges transported by all ionic fluxes: The thicknesses of the depleted and enriched DLs are δ I and δ II respectively; all calculations were performed for δ I = δ II = δ. The origin of coordinates is set at the substrate membrane/modification layer interface (Figure 1a). It is assumed that the concentrations of the components in the bulk solution are known constants; the electric potential is set to zero at the left boundary of the system, and it is set to the value ϕ 0 at the right boundary (potentiostatic electric mode): where d and d ML are the thicknesses of the membrane substrate and modification layer, respectively. The functions γ i c i (x) and ϕ(x) are continuous throughout the four-layer system (between x = −δ I − d ML and x = d + δ II ), including the solution/modification layer interface and the modification layer/membrane interface.
Introduction to the consideration of the activity coefficients in Equation (2) and taking into account their continuous distribution in the layers of the system makes it possible to describe the selective sorption of individual kinds of ions by the membrane-substrate. This is the main difference between this model and those previously developed [22,36].
The system of Equations (2)-(4) with boundary conditions (6)-(8) is a boundary value problem for ordinary differential equations, the numerical solution of which was obtained using the commercially available software package COMSOL Multiphysics 5.6. The input parameters are listed in Table 1. The electric mode is set by the potential difference applied between the external edges of the left-hand and right-hand DLs (Figure 1a). The solution of the problem gives the value of the current density and partial current densities (fluxes) of both kinds of counterions and coions, as well as the distribution of the concentrations of all ions, the potential and the field strength. It is also possible to calculate the space charge density ρ e = F 3 ∑ i=1 z i c i + zQ and the specific permselectivity coefficient, P 1/2 , Equation (1). The calculations are carried out for the case where the equivalent concentrations of both cations in the bulk solutions are the same and equal to 0.02 eq/L. The substrate membrane parameters are close to those of the homogeneous Neosepta CMX cation-exchange membrane studied in ref. [49]. The thicknesses of the depleted and the enriched DLs were selected based on the parameters of the experimental ED cell, which is typical in laboratory studies of electrodialysis ion separation [16].

Input Parameters
All the calculations were carried out using the input parameters presented in Table 1.
Comments and explanations on the methods applied to determine the input parameters are given in.

Limiting Current Density
The developed model describes the transfer of ions through a bilayer membrane bathed in a ternary electrolyte, two kinds of cations (Na + and Ca 2+ ) and one kind of anion (Cl − ). The limiting fluxes of counterions, J i lim , are reached when the concentrations of all three kinds of ions become very small compared to their concentrations in the bulk solution. The fluxes of cations 1 and 2 can be determined by relation (9) [52], which is a generalization of the Peers equation [53]: where J 3 lim is the coion flux density through the membrane, which is directed opposite to the counterion fluxes (back electrodiffusion [54]), the value is taken at the limiting current density; note that z 3 < 0. The term z i J 3 lim δ reflects the fact that chloride ions, which leak back through the membrane and carry a negative charge, cause an additional transport of counterions from the depleted solution to the membrane surface (known as the exaltation effect [52,55]). The summation of the fluxes expressed by Equation (9) (taking into account the charge of the ions) gives the limiting current density for the system under study, j lim : where is the limiting current density in the case of a membrane ideally impermeable to coions [46]; the second term on the right-hand side describes the effect of the exaltation of the limiting current of counterions, where j 3 lim is the limiting partial current density of the coion; is a parameter resembling the coion transport number in solution. Figure 2 shows the specific permselectivity coefficient of a CEM modified with a thin anion-exchange layer as a function of the electric current density at different concentrations of fixed groups in the modification layer, Q ML , and its different thicknesses, d ML . It can be seen that at relatively low current densities, the P 1/2 value increases, reaches a maximum value, and then gradually decreases at higher currents. Similar dependencies were obtained experimentally and theoretically by Golubenko et al. [22].

Specific Permselectivity and CVC
As can be seen from Figure 2, the value of P 1/2 increases with increasing values of d ML and Q ML , as expected: the increase of both parameters leads to a higher resistance of the modification layer towards the divalent cations. In particular, an increase in Q ML results in the increasing Donnan exclusion of cations from the anion-exchange modification layer; since the Donnan exclusion is stronger for divalent coions, the increase in the barrier for the transfer of the divalent cations is essentially greater than for the monovalent cations.
The simulated total and partial current-voltage characteristics (CVCs) of the system under study are shown in Figure 3. It can be seen that at low voltages (up to 0.12 V corresponding to j = 0.4 j 0 lim , where j 0 lim is defined in the previous section), the partial current of Na + ion, j 1 , is lower than that of Ca 2+ , j 2 , that is P 1/2 < 1. However, at j > 0.4 j 0 lim , j 1 > j 2 , and moreover, the rate of growth of j 1 is higher than j 2 up to 0.36 V (j = 0.7 j 0 lim ). The latter is expressed in increasing value of P 1/2 . At j > 0.7 j 0 lim , the situation changes: j 1 increases at a lower rate than j 2 and P 1/2 decreases after passing through its maximum at j = 0.7 j 0 lim . Figure 2 shows the specific permselectivity coefficient of a CEM modified with a thin anion-exchange layer as a function of the electric current density at different concentrations of fixed groups in the modification layer, QML, and its different thicknesses, dML. It can be seen that at relatively low current densities, the P1/2 value increases, reaches a maximum value, and then gradually decreases at higher currents. Similar dependencies were obtained experimentally and theoretically by Golubenko et al. [22]. As can be seen from Figure 2, the value of P1/2 increases with increasing values of dML and QML, as expected: the increase of both parameters leads to a higher resistance of the modification layer towards the divalent cations. In particular, an increase in QML results in the increasing Donnan exclusion of cations from the anion-exchange modification layer; since the Donnan exclusion is stronger for divalent coions, the increase in the barrier for the transfer of the divalent cations is essentially greater than for the monovalent cations.
The simulated total and partial current-voltage characteristics (CVCs) of the system under study are shown in Figure 3. It can be seen that at low voltages (up to 0.12 V corresponding to j = 0.4 j 0 lim, where j 0 lim is defined in the previous section), the partial current of Na + ion, j1, is lower than that of Ca 2+ , j2, that is P1/2 < 1. However, at j > 0.4 j 0 lim, j1 > j2, and moreover, the rate of growth of j1 is higher than j2 up to 0.36 V (j = 0.7 j 0 lim). The latter is  Point 1 on the total CVC corresponds to the typical current density (j = 0.1 j 0 lim), at which mass transfer is controlled by the substrate membrane. The total and partial CVCs in Figure 3 show two inclined plateau regions (the beginning of these plateaus is marked by points 3 and 4, respectively), where a great increase in voltage is needed to obtain a small increase in fluxes of both sodium and calcium ions. However, the increase in partial current of Na + ions between points 3 and 4 is lower than that of Ca 2+ , which indicates a Point 1 on the total CVC corresponds to the typical current density (j = 0.1 j 0 lim ), at which mass transfer is controlled by the substrate membrane. The total and partial CVCs in Figure 3 show two inclined plateau regions (the beginning of these plateaus is marked by points 3 and 4, respectively), where a great increase in voltage is needed to obtain a small increase in fluxes of both sodium and calcium ions. However, the increase in partial current of Na + ions between points 3 and 4 is lower than that of Ca 2+ , which indicates a loss in P 1/2 in this range of currents ( Figure 2). Point 3 on the CVC corresponds to the electric current density (j/j 0 lim ≈ 0.7), at which the maximum permselectivity is observed for the selected input parameters of the model (blue curve in Figure 2a). This point can also be related to the (first) limiting state, which is usually associated in electrochemistry (particularly, in electrodialysis) with the onset of a plateau on the CVC. As we can see in Section 3.3, this state corresponds to the depletion of ions at the modification layer/substrate interface. Point 4 corresponds to the depletion of ions at the solution/modification layer interface. It should also be noted that after the onset of the second limiting state (at j/j 0 lim > 1), a further increase in current density ( Figure 3) is due mainly to the transfer of coions through the substrate membrane and the related exaltation effect (Section 3.1). These effects cause only a slight decrease in P 1/2 at j/j 0 lim > 1 ( Figure 2). Together, the coion transfer and the exaltation effect give a 5% increase in the total current value between ∆ϕ = 17 V (point 4) and ∆ϕ =30 V (not shown in Figure 3).

Concentration Profiles
To explain the behavior of the dependence of the permselectivity coefficient on the current density (Figure 2), the concentration profiles for all three ions have been calculated at electric current densities corresponding to points (1)(2)(3)(4) in Figure 3, in the case where Q ML = 0.5 M, d ML = 20 nm (Figure 4). For better comparison of the ion concentrations in the modification layer, the concentration profiles in it are presented on a logarithmic scale (Figure 4b). The animation of the changes in the ion concentration profiles in the depleted DL with an increase in the potential drop is presented in the Supplementary Materials (Animation S1). Figure 5 shows the dependencies of the differential resistance of individual layers (depleted DL, modification layer and substrate membrane) as a function of the ratio j/j 0 lim . Both the total resistance of the layer k (defined as R tot k = d(∆ϕ tot ) dj ), and the resistance of the layer k with respect to Na + and Ca 2+ ions (R Na k = d(∆ϕ tot ) dj Na + and R Ca k = d(∆ϕ tot ) shown; k is the layer number: k =1 (depleted DL), 2 (modification layer), 3 (substrate layer). Figure 4a shows that in the depleted DL, at current densities of 0.1 j 0 lim and 0.4 j 0 lim , the concentration of Ca 2+ ions at the surface of the modification layer decreases more than the concentration of Na + ions, because the flux of Ca 2+ ions exceeds the flux of Na + ions ( Figure 3). When passing from j = 0.4 j 0 lim to j = 0.7 j 0 lim , the concentration of Na + at the modification layer surface decreases even more, while that of Ca 2+ ions increases (Figure 4a). This is due to the fact that the resistance of the modification layer, especially with respect to Ca 2+ ions, increases sharply with increasing current density in the indicated range. However, at j > 0.7 j 0 lim the resistance of the depleted DL begins to increase, and it increases to a greater extent with respect to Na + (because the concentration of these ions has almost reached its minimum value near the modification layer surface) compared to Ca 2+ . As a result, in this range of current densities, a more significant increase in the calcium flux is observed compared to the sodium flux ( Figure 3). With that, the P 1/2 value decreases sharply and approaches its limiting value, determined by the following equation: Equation (12) is obtained by substituting Equation (9), which give the fluxes of the competing ions at j = j 0 lim , into Equation (1), which defines the P 1/2 quantity. ratio j/j 0 lim. Both the total resistance of the layer k (defined as tot tot k R dj = ), and the resistance of the layer k with respect to Na + and Ca 2+ ions ( ( ) ) are shown; k is the layer number: k =1 (depleted DL), 2 (modification layer), 3 (substrate layer).    As Figure 3 shows, CVCs have two plateau regions that can be associated with the onset of two different limiting states. It follows from the comparison of Figures 3 and 4, that the first limiting state (at j ≈ 0.7 j 0 lim for the parameters presented in the legend to these figures) is reached when the concentration of Na + and Ca 2+ ions in the modification layer/substrate membrane interface (in the bipolar region) reaches critically low values (Figure 4b), causing a significant increase in the potential drop in the modification layer, in the membrane and in the membrane system as a whole when a small increase in current density occurs. A similar phenomenon is observed in systems with bipolar membranes: the limiting current occurs when the ion concentrations at the bipolar boundary approach zero [56,57]. At current densities j > 0.7 j lim , only a slight increase in the Na + flux is observed when there is a significant increase in the flux of doubly charged Ca 2+ ions ( Figure 3). Therefore, the maximum value of P 1/2 is reached at this state (j ≈ 0.7 j 0 lim ); the P 1/2 value decreases, when j > 0.7 j lim (Figure 2a). The Ca 2+ concentration in the electroneutral region of the modification layer is low at all current densities. This region in Figure 4b can be approximately identified as the domain where the reduced concentration of Cl − ions is <1. At j = 0.7 j lim , its value at the right-hand boundary of this layer (always within the electroneutral region) becomes much lower than that at the left-hand boundary (Figure 4b). With increasing j above 0.7 j lim , the concentration of all ions in the bipolar region decreases, which causes an increase in the field strength ( Figure 6a) and an expansion of the space charge region (Figure 6b). Note that in the modification layer, not only the concentrations of coions (Na + and Ca 2+ ), but also the concentration of counterion (Cl − ) decreases nearly to zero (Figure 4b). When the concentrations of all ions approach zero, the ability of the modification layer to be selectively permeable for Na + becomes significantly reduced. With a further increase in j, the concentration of both Na + and Ca 2+ ions in the solution at the surface of the modification layer sharply decreases, the resistance of the depleted DL grows and the control over cation fluxes gradually passes to this DL. As j approaches j lim , the second limiting state occurs when sufficiently low concentrations of all ions are reached at the depleted DL/modification layer interface (Figure 4a). At j = j lim , the P 1/2 value is determined by the parameters of the diffusion layer and, only slightly, by the back electrodiffusion of Cl − ions through the membrane, Equation (12).  The value of P lim 1/2 , found from Equation (12) using the numerically found J 3lim at j = j 0 lim , is equal to 1.10, which is very close to the value of 1.11, found by direct numerical simulation (DNS) at ∆ϕ = 17 V. P lim 1/2 = 1.11 is reached in DNS at ∆ϕ = 40 V. Note that the J 3lim value has a small effect on P lim 1/2 : the latter is equal to 1.12, if we put J 3lim = 0. It should also be noted that the values P lim 1/2 for the modified and unmodified membranes are equal. This qualitatively corresponds to the experimental results reported in references [30,37], as well as to the results of the numerical simulation in reference [22].

Analytical Assessement of the Maximum Value of P 1/2
As our numerical simulation shows, there is a sharp maximum in the dependence of P 1/2 on the current density. This maximum is achieved when the divalent cation concentration vanishes at the right-hand boundary of the modification layer adjacent to the substrate membrane. This state was identified above as the first limiting state of the membrane system. The flux of the monovalent cation in this state reaches a value J 1maxP , which is close to its limiting value attained (in condition of J 3 = 0) at j = j 0 lim , Equation (11). Therefore, taking into account Equation (9), we can write the following equation to evaluate the flux density of the monovalent cations at the point of maximum P 1/2 value: To obtain an assessment of the flux of the divalent cation at this point, J 2maxP , we use Equations (A6) and (A9), deduced in the Appendix A, for the fluxes of this ion in the depleted DL and in the modification layer, respectively, in conditions of the first limiting state. Both fluxes are presented as functions of the divalent cation concentration at the depleted DL/modification layer interface, c s 2 . When equating these fluxes, we obtain an equation, which contains only one unknown quantity c s 2 . In our estimates, we set K z i D = 1, since the activity coefficients of all ions in the DL and the modification layer are equal to 1. Solving this equation (which is a quadratic equation with respect to c s 2 ) yields:

Substituting the solution expressed by Equation (A6) gives
As it follows from Equation (15), to have a maximum possible flux of divalent cations, c s 2 should be as large as possible, approaching c 0 2 . Equation (14) shows that for this, the value of a should be as small as possible. In other words, to obtain a high value of P 1/2 , the modification layer should be thick with a high concentration of fixed charges; the diffusion coefficient of the divalent cation should be small.
Knowing the flux densities of ions 1 [Equation (13)] and 2 [Equation (15)], and applying the definition of P 1/2 [Equation (1)], we can find the maximum value P max 1/2 , using c s 2 calculated from Equation (14): The comparison of the analytical calculations of P max 1/2 , as described above, with the results of the direct numerical simulation (DNS) at different parameters of the modification layer are shown in Figure 7.

Change in the Kinetic Control with Increasing Current Density
As follows from the foregoing, the value of the permselectivity coefficient depends on which layer of the membrane system controls the counterion fluxes. Figure 8 illustrates this relationship. At low current densities, it is the substrate membrane which controls the fluxes since its resistance is the highest (Figure 5). In this current range, the divalent cation is preferably transferred through the bilayer membrane and P1/2 < 1. When the monovalent cation concentration at the modification layer becomes close to zero, and the Donnan exclusion of the divalent cations from this layer is strong enough, the resistance of this layer increases and the kinetic control passes to it. P1/2 becomes greater than 1, and rapidly increases with increasing current density. However, with that, the depletion of ions in the DL growth and its resistance increases correspondingly. The kinetic control passes to this depleted solution layer, which cannot deliver any selectivity to the transport of one of the competing ions. Consequently, the permselectivity is lost when the voltage is sufficiently great.   Figure 7a shows that with a change in the thickness of the modification layer, the analytical calculation has slight discrepancies with the results of the numerical calculation; however, the error does not exceed 10%, which is a very good result, taking into account all the assumptions made during the analytical calculation. A similar situation is observed when the concentration of fixed ionic groups of the modification layer is varied (Figure 7b) but only at Q ML ≤ 0.5 M. At high values of Q ML , there is a great deviation from the analytical calculation from the DNS. There are two main sources of this deviation. When deriving Equations (14)-(16), we assume, first, that the electrically neutral region occupies the main volume of the modification layer, and the ion concentrations at the outer edges of the charged interface are in local equilibrium. Second, a linear concentration distribution is assumed for the competing ions in the depleted DL. However, as Figure 4a shows, the latter is not fulfilled for the divalent cation (Ca 2+ ) at j = 0.7 j 0 lim . The first assumption does not hold when the space charge interfacial regions are relatively large, which occurs at high Q ML values. See also the discussion in Section 3.7.

Change in the Kinetic Control with Increasing Current Density
As follows from the foregoing, the value of the permselectivity coefficient depends on which layer of the membrane system controls the counterion fluxes. Figure 8 illustrates this relationship. At low current densities, it is the substrate membrane which controls the fluxes since its resistance is the highest (Figure 5). In this current range, the divalent cation is preferably transferred through the bilayer membrane and P 1/2 < 1. When the monovalent cation concentration at the modification layer becomes close to zero, and the Donnan exclusion of the divalent cations from this layer is strong enough, the resistance of this layer increases and the kinetic control passes to it. P 1/2 becomes greater than 1, and rapidly increases with increasing current density. However, with that, the depletion of ions in the DL growth and its resistance increases correspondingly. The kinetic control passes to this depleted solution layer, which cannot deliver any selectivity to the transport of one of the competing ions. Consequently, the permselectivity is lost when the voltage is sufficiently great.
clusion of the divalent cations from this layer is strong enough, the resistance of this layer increases and the kinetic control passes to it. P1/2 becomes greater than 1, and rapidly increases with increasing current density. However, with that, the depletion of ions in the DL growth and its resistance increases correspondingly. The kinetic control passes to this depleted solution layer, which cannot deliver any selectivity to the transport of one of the competing ions. Consequently, the permselectivity is lost when the voltage is sufficiently great.   Figure 2 shows that with an increase in the thickness of the anion-exchange modification layer, d ML , the maximum value of P 1/2 of the modified CEM increases; however, this maximum is observed at lower current densities. This dependence is explained by the fact that the presence of a layer on the membrane surface creates a barrier for the transfer of both counterions, but this barrier is higher for polyvalent ions due to their greater charge (the Donnan exclusion effect).

Trade-Off Curve between Membrane Permselectivity and Permeability. The Effect of the Concentration of Fixed Groups in the Modification Layer and Its Thickness
Similarly, the value of P 1/2 is affected by an increase in the concentration of fixed groups in the modification layer, Q ML . Figure 9 allows comparing the P 1/2 − j/j 0 lim dependencies at different values of Q ML . An increase in Q ML (and, as a consequence, the charge density of the fixed ions) leads to an increase in the forces of electrostatic expulsion of cations from the modification layer: the higher the charge number of these ions is, the more significant is the force of expulsion. That is why, with an increase in Q ML , the permselectivity of the modified membrane with respect to singly charged ions increases.  Figure 2 shows that with an increase in the thickness of the anion-exchange modification layer, dML, the maximum value of P1/2 of the modified CEM increases; however, this maximum is observed at lower current densities. This dependence is explained by the fact that the presence of a layer on the membrane surface creates a barrier for the transfer of both counterions, but this barrier is higher for polyvalent ions due to their greater charge (the Donnan exclusion effect).

Trade-Off Curve between Membrane Permselectivity and Permeability. The Effect of the Concentration of Fixed Groups in the Modification Layer and Its Thickness
Similarly, the value of P1/2 is affected by an increase in the concentration of fixed groups in the modification layer, QML. Figure 9 allows comparing the P1/2 − j/j 0 lim dependencies at different values of QML. An increase in QML (and, as a consequence, the charge density of the fixed ions) leads to an increase in the forces of electrostatic expulsion of cations from the modification layer: the higher the charge number of these ions is, the more significant is the force of expulsion. That is why, with an increase in QML, the permselectivity of the modified membrane with respect to singly charged ions increases. In the previous sections, it was shown that the maximum value of P1/2 ( max 1/ 2 P ) is reached when the sodium ion flux reaches its near-limiting value: when the current density/potential drop increases further, the flux of Na + ions increases slightly compared to In the previous sections, it was shown that the maximum value of P 1/2 (P max 1/2 ) is reached when the sodium ion flux reaches its near-limiting value: when the current density/potential drop increases further, the flux of Na + ions increases slightly compared to the flux of Ca 2+ . The dependency of the sodium ion partial current density along with the total current density as functions of the applied voltage are presented in Figure 10. The values of the partial current densities of Na + and Ca 2+ ions (j' Na and j' Ca respectively) corresponding to the maximum value of the specific permselectivity coefficient, P max 1/2 , at different values of d ML and Q ML are presented in Table S1 (see Supplementary Materials). These results are shown also in Figure 11. Table S1 and Figure 11 testify that an increase in the thickness and concentration of the fixed groups in the modifying layer leads to an increase in the membrane permselectivity but also leads to a decrease in the flux of the preferably transferred ion. The dependence of P max 1/2 on the j' Na , calculated at different values of d ML and Q ML , is the so-called trade-off curve between membrane permselectivity and permeability for the case of monovalent-ion selective bilayer membranes. This kind of curve was first proposed by Robeson [58] as an empirical "upper bound" for gas-separation membranes. Later on this "upper bound" curve was repeatedly used to analyze the transport characteristics of many other types of membranes [13,18,19,22]. A short formulation of this trade-off relationship was formulated by Park et al. [19]: "highly permeable membranes lack selectivity and vice versa". transport characteristics of many other types of membranes [13,18,19,22]. A short form lation of this trade-off relationship was formulated by Park et al. [19]: "highly permea membranes lack selectivity and vice versa".  (a) (b) Figure 10. Total CVC of the system and partial CVCs of Na + at different parameters of the modification layer: at different dML (QML = 0.5 М) (a) and at different QML (dML = 20 nm) (b). Figure 11. Dependence of the maximum value of the specific permselectivity coefficient on the partial current density of Na + . Red dots are the results of simulation at different values of dML and QML (Table S1); the black line is the fitting curve for these results.

Discussion
The application of the quasi-electro-neutral version of the model allows obtaining an analytical assessment of P1/2 max , as mentioned above. However, at high values of the modification layer thickness and the concentration of fixed charges in this layer, the analytical assessment of P1/2 max deviates from the values found in DNS. One of the reasons for the deviation is the use of the local electroneutral assumption. In particular, this  (Table S1); the black line is the fitting curve for these results.

Discussion
The application of the quasi-electro-neutral version of the model allows obtaining an analytical assessment of P 1/2 max , as mentioned above. However, at high values of the modification layer thickness and the concentration of fixed charges in this layer, the analytical assessment of P 1/2 max deviates from the values found in DNS. One of the reasons for the deviation is the use of the local electroneutral assumption. In particular, this assumption stipulates that the ion concentrations at the outer edges of the charged interfaces between the modification layer and both neighboring layers are in local equilibrium. However, this equilibrium is not held when the thickness of these interfaces (which are interfacial electrical double layers) is relatively high. When deriving the equilibrium condition, it is usually assumed that the diffusion and electromigration components of an ion flux through the interface are both large in absolute value and opposite in sign. Then the difference between these components in the Nernst-Planck equation can be neglected, which leads to the Boltzmann equation and constancy of electrochemical potential for a given ion. The latter leads to the Donnan relations [59]. However, in the considered system the situation is rather special: the cations are coions for the modification layer, hence their fluxes usually should be low. However, it is not the case, and the cation fluxes are high; Figure 4a,b shows that the diffusion and electromigration components of the Ca 2+ flux have the same direction in the left-hand SCR of the modification layer. Note also that a linear concentration distribution was assumed for the competing ions in the depleted DL. This assumption is not valid for the Ca 2+ ion as well: due to the high resistance of the modification layer regarding this ion, its concentration increases when approaching the surface of the modification layer. Evidently, the latter effect of the concentration profile non-linearity increases with increasing the diffusion resistance of the modification layer towards cations, that is, with decreasing the value of the a = D ML 2 (1+|z 2 /z 3 |) Q ML d ML parameter entering Equation (14). Therefore, we can see that the application of the local electroneutrality condition allows only an approximate description of the behavior of a bilayer membrane with respect to its permselectivity for monovalent ions. The adequacy of this version deteriorates with increasing modification layer resistance and the P 1/2 max parameter values become significantly less than those given by DNS when using the Poisson equation (Figure 6). The Poisson equation makes it possible to take into account specific behavior of ions within the interfacial space charge regions at high current densities. However, deviations from adequacy are also possible when describing ion transport on a scale of less than 1 nm due to the fact that a model based on the NPP equations rests on the hypothesis of a dilute solution of point-like ions in a continuous medium. However, on such a small scale, ionic size becomes important (e.g., ionic crowding against a blocking surface), which must be taken into account when describing some phenomena such as the differential capacitance of the electrical double layer [60] and electroosmotic flow in (bio)microfluidic systems [61]. For this, modified electrokinetic equations for finite-sized ions can be used. Bazant et al. [60] gave a general framework of approaches involving such equations.
There is a certain similarity in the structure and behavior of artificial ion-exchange membranes and biological cell membranes. Both have pores/channels which are selectively permeable for some types of ions. Charge selectivity of cell membranes can be mediated by the electrostatic interaction of partially dehydrated permeating cations with negatively charged sites within a pore that is formed by side-chain carboxyl groups [62]. We study here an artificial composite bilayer membrane, one (thicker) layer of which bears fixed anions and another (thinner) layer has cations fixed to its matrix. One can imagine this membrane as a system of parallel long pores/channels. A short segment of such a pore would bear positive fixed charges on its walls, and its longer part would bear negative fixed charges on the wall. When two competing cations, such as Na + and Ca 2+ , enter the short segment of the pore, the Ca 2+ ions are repulsed by the positive charge of the walls, while Na + can overcome this segment under the action of the gradient of concentration or electric potential. Furthering its way through the longer part of the pore with negatively charged walls would be rather easy for Na + . This long segment serves to stop the transfer of negatively charged anions, such as Cl − . The model presented in this paper describes monovalent-selective ion transport through such a composite pore/channel. We believe that this model can be useful not only for specialists in artificial membranes but also for those who study selective transport of ions through cell membranes. Mono-and divalent cation voltage-dependent permeation occurs in cyclic nucleotide-gated channels [63][64][65][66]; permselective mass transfer takes place in the plasma membranes of thermophilic bacteria [67] or in the aquaporin channels of cell membranes [68]. In particular, a simple expression for approximative evaluation of mono-divalent permselectivity can be useful for better understanding the causes and membrane parameters which control this selectivity.

Parameters of the Substrate Membrane and Modification Layer
The thickness, d, and the concentration of fixed groups, Q m , in the substrate membrane are taken the same as those for the commercial Neosepta CMX homogeneous cationexchange membrane (produced by Astom, Tokuyama Soda, Japan) [50]. This membrane is preferably permeable to multivalent cations [49].
As mentioned above, the diffusion coefficients of ions in the substrate membrane, D m i , were calculated from the electrical conductivity and diffusion permeability of the Neosepta CMX membrane [50]. In the modification layer, the diffusion coefficients are taken 3 times fewer than in solution. This value was chosen from the consideration that usually the modification layer is not as dense as the substrate membrane, so the diffusion coefficients in it are only slightly less than those in solution [22,36].

Determination of Activity Coefficients
The relationship between the activity coefficients of counterions in solution and membrane follows from the condition of continuity of activities and electric potential at the solution/membrane interface, which leads to the ion-exchange equilibrium equation (also called the Nikolsky's equation [69]): where c i and c i are the molar concentrations of ions i in the membrane and solution, respectively. The thermodynamic equilibrium constant K 21 can be expressed in terms of the ion activity coefficients [45]: However, when considering the laws of ion exchange for the membrane and solution phases, it is advisable to choose equivalent fractions of ions as concentration units, since only in this case are the activity coefficients of the components in the solution always equal to unity [45]. Then Equation (17) takes the form: is ion exchange equilibrium coefficient. For ion exchangers with sulfonate fixed groups K 21 > 1, where subscript 1 refers to monovalent and subscript 2, to divalent cation [45]. A detailed derivation of Equations (17) and (19) is presented in the Supplementary Materials. As was written earlier, the activity coefficients of the ions in the solution and in the modification layer are equal to unity; the activity coefficients of the counterions in the substrate membrane, γ m i , were selected so that the equivalent fraction of Ca 2+ in the substrate membrane is approximately 20 times higher than Na + (K 21 ≈ 20), when their equivalent concentrations in the equilibrium solution are the same and equal to 0.02 eq/L.

Conclusions
Within the framework of the developed 1D model based on the Nernst-Planck-Poisson equations, it was shown that the dependence of the specific permselectivity of the modified membrane on the electric current density passes through a maximum. The P 1/2 − j dependence in this model is caused only by the transitions of the kinetic control of the transfer of competing ions from one layer of the multilayer membrane system to another. No changes in the concentration of fixed groups in the membrane-substrate and in the modification layer occur (these parameters are considered constant); accordingly, structural changes in the modification layer do not take place.
At low current densities, the transfer of cations through a CEM modified with a thin anion-exchange layer is controlled by the substrate membrane due to its relatively large thickness and resistance. Since this layer selectively absorbs doubly charged cations, these cations are selectively transported through the multilayer system. The increase in the specific permselectivity of the modified membrane with increasing current density is due to the transition of the kinetic control of mass transfer from the substrate membrane to the anion-exchange modification layer, which is a very significant barrier to calcium ion transfer and a less significant barrier to sodium ion transfer. This increase is observed until the system reaches the first limiting state due to a drop in the concentrations of Na + and Ca 2+ near the modifying layer/membrane-substrate interface (near the bipolar region) to critically low values by analogy with bipolar membranes. After reaching the first limiting state, the transfer of sodium ions almost does not increase, while a noticeable increase in the flow of calcium ions is observed. In this case, the concentrations of all ions (especially calcium ions) in the solution at the boundary with the modification layer rapidly decrease, and the kinetic control of mass transfer passes from the modification layer to the depleted diffusion layer. The second limiting state occurs when the concentrations of all ions at the depleted DL/modification layer interface reach critically low values. In this case, the value of P 1/2 reaches the same value that it reaches for the unmodified substrate membrane: the value of P 1/2 is determined by the parameters of the diffusion layer and depends very weakly on the reverse transfer of chloride anions from an enriched solution to a depleted solution.
An increase in the thickness and concentration of the fixed groups of the modification layer, on the one hand, allows a significant reduction in the flux of multivalent ions and an increase in membrane monovalent permselectivity. However, this also leads to a decrease in the flux of monovalent ions, and thus to a reduction in membrane permeability. This relationship is expressed by the simulated trade-off curve between membrane permselectivity and permeability for the case of monovalent-ion selective bilayer membranes.