Mathematical Modeling of the Effect of Pulsed Electric Field on the Specific Permselectivity of Ion-Exchange Membranes

The application of pulsed electric field (PEF) in electrodialysis has been proven to be efficient for a number of effects: increasing mass transfer rate, mitigation of scaling and fouling, reducing water splitting. Recently, the improvement of the membrane permselectivity for specific counterions was discovered experimentally by the group of Laurent Bazinet (N. Lemay et al. J. Memb. Sci. 604, 117878 (2020)). To better understanding the effect of PEF in electrodialysis, simulations were performed using a non-stationary mathematical model based on the Nernst–Planck and Poisson equations. For the first time, it was not only the condition used when the current density is specified but also the condition when the voltage is set. A membrane and two adjacent diffusion layers are considered. It is shown that when applying the regime used by Lemay et al. (the same current density in conventional continuous current (CC) mode and during the pulses in PEF mode), there is a significant gain in specific permselectivity. It is explained by a reduction in the membrane concentration polarization in PEF mode. In the CC mode of electrodialysis, increasing current density leads to a loss in specific permselectivity: concentration profiles in the diffusion layers and membrane are formed in such a way that ion diffusion reduces the migration flux of the preferentially transferred ion and increases that of the poorly transferred ion. In PEF mode, the concentration profiles are partially restored during the pauses when the current is zero. However, if a different condition is used than the condition applied by Lemay et al., that is, when the same average current density is applied in both the PEF and CC modes, there is no gain in specific permeability. It is shown that within the framework of the applied mathematical model, the specific selectivity depends only on the average current density and does not depend on the mode of its application (CC or PEF mode).


Introduction
Electrodialysis (ED) is one of the rapidly progressing membrane methods for desalination, concentration and separation of aqueous solutions today. Many years of application of this method in various industries has proven its efficiency [1][2][3][4]. Although ED is a mature method, it has some hindrances. Most of them are associated with the phenomenon of concentration polarization, which is caused by the difference in the transport numbers of ions in the solution and the ion-exchange membrane. The effect leads to a decrease in the electrolyte concentration in a thin layer near the membrane surface. The growth of the system resistance results in the current density decrease with time in the potentiostatic mode or an increase in the potential drop with time in the galvanostatic mode [5,6]. Depending on the mode, the system tends to a certain stationary value of the current density or potential drop. This value can be considered as a maximum of ED productivity under given conditions (constant current or potential drop). The concentration polarization provides such undesirable effects as water splitting [7][8][9][10], membrane scaling [7,11,12] and fouling [13]. 2 of 16 The use of current or voltage pulses alternating with pauses in ED, the so-called PEF mode, leads to the mitigation of concentration polarization and, as a result, improves the performance of ED compared with conventional CC mode [5,7,8,[13][14][15]. In particular, it is found that the use of PEF mode allows increasing mass transfer rate [16,17], can essentially reduce water splitting [8,17], membrane scaling [7,18] and fouling [13,19,20], as well as membrane stack resistance [21].
Mishchuk et al. [5] showed theoretically that the characteristic pulse time should be less than the transition time required to build up the concentration polarization layer near the membrane. The pause duration should be comparable to the pulse duration to enable the concentration profiles that can be restored to their initial state. However, the pause should not be too long; otherwise, diffusion and osmotic fluxes in the membrane counteract the desalination process. The gain in mass transfer in PEF mode mainly occurs at the beginning of the pulse lapse due to the appearance of a current spike (or voltage fall) for a short period of time. Based on a mathematical description using the Nernst-Planck equations, Sistat et al. [16] explain that this current spike (or voltage fall) is due to the partial concentration restoration at the vicinity of the membrane during the pause lapse. Moya and Moleón [22] applied the network simulation method for examining the time variation of the counterion flux leaving the membrane in response to a pulsed electric potential as well as some other parameters in order to establish the conditions under which the application of PEF can improve the performance of electrodialysis.
Uzdenova et al. [23] theoretically studied the effect of a PEF on mass transfer at overlimiting current regimes by applying the Nernst-Planck-Poisson-Navier-Stokes equations. They showed that the intensive electroconvective vortices decay almost immediately after switching off the current. However, the residual weak vortices continue to exist near the membrane surface for a few tenths of a second, fed by a non-uniform concentration field. After the pause, when the voltage pulse resumes, the first weak vortices appear already in 0.01 s. This rapid recovery of electroconvection is explained by the inhomogeneity of the residual concentration field. The inhomogeneity stimulates the development of electroconvection in the same way as it occurs in the case of an electrically heterogeneous surface [24]. Later, similar effects were obtained experimentally and described by Butylskii et al. [25]. Zyryanova et al. [26] showed that in overlimiting PEF modes, the average current density increases up to 33% over the period when applying the same average voltage in PEF and CC modes. Lemay et al. [27] showed that the use of PEF mode allows reducing energy consumption for the demineralization process. Thus, these results showed that the application of the PEF modes in ED is a promising way for the process intensification in the overlimiting current range. It is hypothesized [26,27] that the main cause of improvement of ED performance when applying PEF in intensive current regimes is electroconvection.
Currently, the main effect of PEF on the process of ED treatment of solutions is the mitigation or suppression of the membrane scaling and fouling. This effect was well established experimentally by a number of authors [7,11,13,[18][19][20][28][29][30]. The explanation of the impact of PEF is as follows: -Reduction of concentration polarization. Relaxation of the concentration profile occurs at the membrane surface during the pause lapse. The concentrations of ion species resume partially or completely to the initial values; -Intensification of electroconvection at overlimiting current regimes, which in addition to the increase in mass transfer helps to wash out the scale and foulant components from the membrane surface; -Reduction of water splitting (pH values at which some components can precipitate are not reached).
Recently, Lemay et al. [17] experimentally established another important positive effect of PEF mode on ED characteristics. It was found that the use of PEF mode allows increasing the specific permselectivity of ion-exchange membranes (preferential transport of counterions of one kind over counterions of another kind). In CC mode, a high specific permselectivity (sometimes allowing the flux of one kind of ions to be >100 times higher Membranes 2021, 11, 115 3 of 16 than the flux of another one at the equal bulk concentrations [31,32]) can be obtained only at low current densities, i. When i approaches its limiting value, i lim , the specific selectivity is lost [33,34] due to increasing concentration polarization and the transition of the ion transfer control from the membrane to the depleted diffusion layer [35]. The possibility to have a high permselectivity at elevated current densities is of great practical interest. The authors associated the enhancement of permselectivity by PEF with the partial restoration of ion concentration profiles in the depleted boundary solution during the pause. As a consequence, the membrane partially restores control over the kinetics of ion separation. Earlier [17], was studied only the mode where a constant electric current was applied during the pulses. The concentration profiles were calculated, and it was found that during one pulse, there is a gain in the membrane specific permselectivity compared to the conventional continuous current ED mode. However, the averaged in time-specific permselectivity was not quantified. In this paper, a simulation was carried out using a non-stationary mathematical model based on the Nernst-Planck and Poisson equations. For the first time, two PEF modes are compared when current or voltage pulses are applied.

Experimental Results
Here are some of the experimental results of Lemay et al. [17] showing enhancement of the specific permselectivity by applying PEF in ED of sweet whey. The feed solution was obtained by dissolution of sweet whey powder in distilled water; the total solids content of the obtained solution was 6.5%. The kinetics of desalination of this solution was studied [17] in a batch ED process using a constant current density of 8.0 mA/cm 2 . According to the estimation by the method of Cowan and Brown, the limiting current density was 13.5 mA/cm 2 ; hence the ratio of the current density to its limiting value at the beginning of the desalination process was i/i lim ≈ 0.6. However, as far as the solution was desalted, the value of i lim decreased. When the total degree of desalination was 60%, i approached i lim , and i/i lim became approximately equal to 1.5. The same current density of 8.0 mA/cm 2 was applied during the pulses in PEF mode, while the current was zero during the pauses. Table 1 compares the demineralization rates found for different ions constituting the sweet whey at two different total degrees of desalination (DD). The ratios of DD by individual ions to the total DD of sweet whey, θ k , are given in CC mode and in PEF mode for the pulse/pause combination, 1 s-1 s (0.5 Hz). The greater the θ k value, the greater the transport number of ion k (the fraction of the electric charge transported by this ion through the membrane).
The K + ion is the dominant cation in the sweet whey. Since its mobility in the membrane is relatively high, the demineralization rate for this ion is the highest. The DD by this ion is higher than the total DD, θ K > 1, and this parameter changes a little with decreasing the total concentration and increasing the i/i lim ratio. In CC mode, the values of θ Na and θ Ca are close at i/i lim ≈ 0.6, but θ Na is significantly higher than θ Ca at i/i lim ≈ 1.5. That is, with an increase in the i/i lim ratio, the selectivity of Ca 2+ transfer through the membrane decreases in favor of Na + . A similar situation occurs when comparing the competitive transport of Na + and Mg 2+ . In the study by Lemay et al., a cation-exchange membrane was used whose properties were close to that of the Neosepta CMX membrane [17]. It is a homogenous membrane (produced by Astom, Tokuyama Soda, Japan) made by the "paste method" [36,37]. Initially, the paste contains styrene monomer with functional groups (which are subsequently grafted with ion-exchange groups), divinylbenzene (45-65%) as a crosslinking agent, a radical polymerization initiator and powdered polyvinyl chloride (45-55%). The paste is deposited on the reinforcing polyvinyl chloride fabric. The copolymerization is carried out before the sulfonation. The membrane consists of two interpenetrating phases: ion exchange material and polyvinyl chloride, PVC, having a particle diameter of 100 nm or less [38]. The structural inhomogeneities in the membrane volume do not exceed 1 micron [39]. An exception is a reinforcing fabric having fibers of 25-30 µm in diameter [40]. The results presented in Table 1 show that in CC mode, this permselectivity decreases both in the cases of Na + and Ca 2+ and Na + and Mg 2+ . However, in PEF mode, the membrane permselectivity for the Ca 2+ and Mg 2+ cations is significantly greater than that in CC mode if the same current density is applied in CC mode and in PEF mode during the pulses. Note that the average current density applied in PEF mode was two times lower since, during the pauses, the current was zero. It is important that the total duration of the ED process needed to obtain total DD ≈ 70% is nearly the same in both modes [17], whereas in PEF mode, the current flows across the membranes only for a certain fraction of the total processing time.

Theoretical Part
The system under study has a three-layer geometry (an ion-exchange membrane and two adjacent diffusion boundary layers, DBLs) ( Figure 1). The membrane is 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. The impact of convection transport in solution is taken into account implicitly through the diffusion layer thickness, which is considered independent of the current/voltage applied. Water splitting and electroconvection are not taken into account. θNa and θCa are close at i/ilim ≈ 0.6, but θNa is significantly higher than θCa at i/ilim ≈ 1.5. That is, with an increase in the i/ilim ratio, the selectivity of Ca 2+ transfer through the membrane decreases in favor of Na + . A similar situation occurs when comparing the competitive transport of Na + and Mg 2+ . In the study by Lemay et al., a cation-exchange membrane was used whose properties were close to that of the Neosepta CMX membrane [17]. It is a homogenous membrane (produced by Astom, Tokuyama Soda, Japan) made by the "paste method" [36,37]. Initially, the paste contains styrene monomer with functional groups (which are subsequently grafted with ion-exchange groups), divinylbenzene (45-65%) as a crosslinking agent, a radical polymerization initiator and powdered polyvinyl chloride (45-55%). The paste is deposited on the reinforcing polyvinyl chloride fabric. The copolymerization is carried out before the sulfonation. The membrane consists of two interpenetrating phases: ion exchange material and polyvinyl chloride, PVC, having a particle diameter of 100 nm or less [38]. The structural inhomogeneities in the membrane volume do not exceed 1 micron [39]. An exception is a reinforcing fabric having fibers of 25-30 µm in diameter [40]. The results presented in Table 1 show that in CC mode, this permselectivity decreases both in the cases of Na + and Ca 2+ and Na + and Mg 2+ . However, in PEF mode, the membrane permselectivity for the Ca 2+ and Mg 2+ cations is significantly greater than that in CC mode if the same current density is applied in CC mode and in PEF mode during the pulses. Note that the average current density applied in PEF mode was two times lower since, during the pauses, the current was zero. It is important that the total duration of the ED process needed to obtain total DD ≈ 70% is nearly the same in both modes [17], whereas in PEF mode, the current flows across the membranes only for a certain fraction of the total processing time.

Theoretical Part
The system under study has a three-layer geometry (an ion-exchange membrane and two adjacent diffusion boundary layers, DBLs) ( Figure 1). The membrane is 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. The impact of convection transport in solution is taken into account implicitly through the diffusion layer thickness, which is considered independent of the current/voltage applied. Water splitting and electroconvection are not taken into account.
where J k , c k , D k , z k , and γ k are the flux density, concentration, diffusion coefficient, charge number and activity coefficient of ion k (k = Na + , Ca 2+ , Cl − ), respectively; R is the gas constant; T is the temperature; F is the Faraday constant; ϕ is the electric potential; ε 0 is the vacuum permittivity; ε is the solution relative permittivity; E = − ∂ϕ ∂x is the electric field strength; z m is the charge number of the membrane; Q is the exchange capacity of the membrane; F is the Faraday constant; t is the time.
The electric current density, i, in the system is equal to the sum of ionic fluxes (expressed in appropriate units): Equation (4) does not take into account the displacement electric current. In the present work, the processes with a characteristic time of more than 10 −2 s are considered, while the displacement current occurs in a very short time interval, less than 10 −4 s. This process should be accounted for when considering, for example, a high-frequency impedance of membrane systems [41].
Let δ I and δ II be the thicknesses of the depleted and enriched diffusion layers, respectively, and the origin of coordinates is at the solution/membrane boundary ( Figure 1). Then, the boundary conditions for the system under study will be: where c k0 is the concentration of ion k in the bulk solution. Equation (5) sets the concentration of all ions in the electroneutral bulk solutions on both sides of the membrane. At the membrane/solution interfaces, the continuity of concentrations and electric potential are assumed, i.e., the functions c k (x, t) and ϕ (x, t) are continuous on the interval (x = −δ I , x = d + δ II ). The potential is assumed zero at the left edge of the three-layer system (Equation (6)).
Earlier, we pointed out that our study uses two PEF modes when current or voltage pulses are applied. Equations (1)-(6) are common for both cases.
However, when current pulses are applied, the following boundary condition is used at the right-hand edge of the three-layer system: where the electric current density, i, is set as a known function of time; it takes a given value during a pulse lapse and is zero during a pause lapse. Equation (7) follows from the continuity of current density and Equations (1) and (4), which were proposed by Uzdenova et al. [42].
In the PEF mode, when voltage pulses are applied, the following boundary condition is used at the right-hand edge of the three-layer system during the pulse: where U pulse and U av are the voltage applied during the pulse, and the average voltage, respectively; α is the duty cycle. During the pause, the current density is zero and condition (7) is used, in which i is set to zero, i = 0.
Thus, the mathematical formulation of the problem for PEF mode, when current pulses are applied, is described by Equations (1)-(7); for the case where voltage pulses are applied, it is described by Equations (1)-(6), (8) during the pulse lapse and Equations (1)- (7) during the pause lapse. The numerical solution was obtained using the Comsol Multiphysics software.
For the calculation in the PEF mode, the following parameters were taken: the duty cycle, α, is equal to 0.5, and the PEF frequency, f, is equal to 0.5 Hz and 10 Hz. The authors [17] noted that under these conditions, an increase of the specific permselectivity of the membrane is observed. The frequency of 10 Hz was taken since, as follows from the work of Sistat et al. [16], the mass transfer rate in PEF mode increases with increasing frequency, and at f ≥ 10 Hz, this rate becomes nearly constant.
The model allows calculation of the ion concentration profiles and the electric current density in the membrane system. It is also possible to calculate the partial current densities of ions and their effective transport numbers, T k , in the membrane. T k is defined as a fraction of current transported by ion k: where z k J k F = i k is the partial current density of ion k.

Results and Discussions
All the calculations are carried out using the input parameters presented in Table 2. The output parameters are time-dependent concentrations profiles, total current density (when the voltage is the input parameter) or the voltage (if the current density is given) and the partial current densities, i k (or effective transport numbers, T k ), computed in the membrane at a distance of 1 µm from the depleted solution/membrane interface (left-hand interface). The latter shows the rate of ion transport through the left-hand interface. The point located at 1 µm is chosen since, at a shorter distance, the errors in the calculated values of i k (T k ) are relatively high due to very large concentration and potential gradients. The diffusion coefficients of Na + and Ca 2+ in the membrane are determined from the electrical conductivity of a Neosepta CMX membrane [17]. The activity coefficients of the cations in the membrane are chosen from the following consideration. The activity coefficients in the membrane and solution are linked with the ion-exchange constant, K N (which is also called the Nikolsky's constant [43,44]) as follows [45,46]: K N enters the relation of local thermodynamic equilibrium, which holds at the solution/membrane interfaces: The activity coefficients in the solution (γ 1 and γ 2 ) were taken equal to 1; the activity coefficients in the membrane (γ 1 and γ 2 ) were found in a way that the equivalent fraction of Ca 2+ in the membrane is about 10 times higher than that of Na + (as experiment shows in the case of ion-exchangers with sulfonate fixed groups [45]), when their equivalent fractions in the equilibrium bathing solution are the same at I = 0 and when the concentrations of both ions in the solution were 0.02 eq/L, which approximately corresponds to the experimental conditions [17].
The limiting current density, i lim , of the system was calculated using Equation (12). This equation is applied for systems with ternary electrolyte (e.g., CaCl 2 and NaCl), under the assumption that the membrane is impermeable to coion [47]: here index k is related to counterions (Na + or Ca 2+ ) and A to coion (Cl − ); δ is the thickness of the depleted diffusion layer (δ = δ I ), which was estimated from the hydrodynamic parameters of the experimental system [17]. Since the model does not take into account all the components of the experimental solution, in particular K + , which is the dominant cation [17], the calculated limiting current density i lim = 6.52 mA/cm 2 is lower than the experimental one (13.5 mA/cm 2 ). However, in this paper, we are not aiming for quantitative agreement with the experiment. The goal is to understand whether taking into account electromigration and diffusion is sufficient to obtain a qualitatively correct picture of the effect of PEF on membrane permselectivity for a specific ion.

CC Mode
The case where constant voltage pulses alternate with pauses, during which the current density is zero (i = 0), is considered since it is often used in practice [14,16]. In this case, the results are compared at the same value of the time-average voltage, U av , which is linked with the voltage during the pulse, U pulse , as U av = αU pulse .
A relatively high concentration of Ca 2+ in the membrane (due to great value of K N , Table 2) leads to a higher flux of this ion compared to the Na + flux in conditions where the equivalent concentrations of both ions in the solution bulk are the same; thus, the membrane shows permselectivity for specific ions [48], the divalent Ca 2+ cation in the considered case. However, with increasing current density, the permselectivity decreases, and it is nearly completely lost at I = i lim [33,34,47]. This loss is due to a particular developing concentration polarization of the membrane. Since Ca 2+ is the preferentially transferred ion over Na + , the applied electric current causes a faster decrease in the concentration of Ca 2+ at the depleted solution/membrane interface than the decrease in the concentration of Na + . Moreover, the computation shows (and experimental data [49,50] confirm that) that at low voltages (e.g., at U av = 111 mV in Figure 2a), the concentration of Na + at the depleted interface can be even higher than its bulk value. Thus, under an applied current, the equilibrium at this interface shifts in favor of Na + , whose relative concentration in the depleted solution and in the membrane increases. As a result, the ratio between the fluxes of Ca 2+ and Na + (and between the transport numbers) changes in favor of Na + : preferential transfer of Ca 2+ at low current densities (low voltages) is lost, the transfer rate of Na + at i = i lim becomes even greater than that of Ca 2+ (Figure 3). At the limiting current density, the fluxes of the competing ions are controlled by the depleted diffusion layer and do not depend on the membrane properties Equation (9); the role of the membrane in terms of selective transfer is reduced to being a barrier for co-ions.
of Na + . Moreover, the computation shows (and experimental data [49,50] confirm that) that at low voltages (e.g., at Uav = 111 mV in Figure 2a), the concentration of Na + at the depleted interface can be even higher than its bulk value. Thus, under an applied current, the equilibrium at this interface shifts in favor of Na + , whose relative concentration in the depleted solution and in the membrane increases. As a result, the ratio between the fluxes of Ca 2+ and Na + (and between the transport numbers) changes in favor of Na + : preferential transfer of Ca 2+ at low current densities (low voltages) is lost, the transfer rate of Na + at i = ilim becomes even greater than that of Ca 2+ (Figure 3). At the limiting current density, the fluxes of the competing ions are controlled by the depleted diffusion layer and do not depend on the membrane properties Equation (9); the role of the membrane in terms of selective transfer is reduced to being a barrier for co-ions.  The shape of the Tk vs. Uav dependence is in good agreement with the literature data [33,34]. of Na + . Moreover, the computation shows (and experimental data [49,50] confirm that) that at low voltages (e.g., at Uav = 111 mV in Figure 2a), the concentration of Na + at the depleted interface can be even higher than its bulk value. Thus, under an applied current, the equilibrium at this interface shifts in favor of Na + , whose relative concentration in the depleted solution and in the membrane increases. As a result, the ratio between the fluxes of Ca 2+ and Na + (and between the transport numbers) changes in favor of Na + : preferential transfer of Ca 2+ at low current densities (low voltages) is lost, the transfer rate of Na + at i = ilim becomes even greater than that of Ca 2+ (Figure 3). At the limiting current density, the fluxes of the competing ions are controlled by the depleted diffusion layer and do not depend on the membrane properties Equation (9); the role of the membrane in terms of selective transfer is reduced to being a barrier for co-ions.  The shape of the Tk vs. Uav dependence is in good agreement with the literature data [33,34]. The shape of the T k vs. U av dependence is in good agreement with the literature data [33,34].

PEF Mode (Current Pulses)
Typical time dependences of the current density and voltage for the mode under consideration are shown in Figure 4

PEF Mode (Current Pulses)
Typical time dependences of the current density and voltage for the mode under consideration are shown in Figure 4.  The calculation shows that in the conditions of the experiment [17], when the same current density in CC mode and during the pulse is applied, iCC = ipulse = 0.6ilim, in CC mode, the TCa:TNa ratio is 0.52:0.48, while in PEF mode this ratio is essentially higher: 0.65:0.34; the parameters of the PEF mode are as follows: the duty cycle α = 0.5 and frequency f = 0.5 Hz. This theoretical result is in qualitative agreement with the experiment [17] described above (Section 2). The PEF effect is due to partial restoration of concentrations during a pause lapse. The restoration occurs near the surface; a little increase in concentration produces an important decrease in resistance. As a consequence, after a pause, the profiles are closer to those, which occur at a lower CC current density, which explains a higher permselectivity. An important feature is that the Na + partial current density is negative during the pause lapse, while the Ca 2+ partial current density remains positive (Figure 5a). This explains the increase in permselectivity attained in the PEF mode. However, when calculating the values of the transport numbers in the CC and PEF modes under the condition that the average over the period current density iav in PEF mode is equal to the constant current density in CC mode, no difference was found between the transport numbers in CC and PEF modes. When we have taken iav = iCC = 0.3ilim (in this case the current density during the pulse is ipulse = 0.6ilim), α = 0.5 and f = 0.5 Hz, we find that the TCa:TNa ratio is the same in both modes and equal to 0.65:0.34. The same TCa:TNa ratio (0.65:0.34) is found at a higher frequency f = 10 Hz (α = 0.5), when ipulse = 0.6ilim. The calculations show that the TCa:TNa ratio depends only on the iav value and does not depend on the frequency and duty cycle; it does not matter if PEF or CC mode is used.

PEF Mode (Voltage Pulses)
Typical time dependences of the voltage and current density for the PEF mode, where constant voltage pulses alternate with the pauses during which the current is zero, are shown in Figure 6.
The time dependences of the partial current densities of Na + and Ca 2+ ions at a frequency of 0.5 Hz are shown in Figure 7. The calculation shows that in the conditions of the experiment [17], when the same current density in CC mode and during the pulse is applied, i CC = i pulse = 0.6i lim , in CC mode, the T Ca :T Na ratio is 0.52:0.48, while in PEF mode this ratio is essentially higher: 0.65:0.34; the parameters of the PEF mode are as follows: the duty cycle α = 0.5 and frequency f = 0.5 Hz. This theoretical result is in qualitative agreement with the experiment [17] described above (Section 2). The PEF effect is due to partial restoration of concentrations during a pause lapse. The restoration occurs near the surface; a little increase in concentration produces an important decrease in resistance. As a consequence, after a pause, the profiles are closer to those, which occur at a lower CC current density, which explains a higher permselectivity. An important feature is that the Na + partial current density is negative during the pause lapse, while the Ca 2+ partial current density remains positive (Figure 5a). This explains the increase in permselectivity attained in the PEF mode. However, when calculating the values of the transport numbers in the CC and PEF modes under the condition that the average over the period current density i av in PEF mode is equal to the constant current density in CC mode, no difference was found between the transport numbers in CC and PEF modes. When we have taken i av = i CC = 0.3i lim (in this case the current density during the pulse is i pulse = 0.6i lim ), α = 0.5 and f = 0.5 Hz, we find that the T Ca :T Na ratio is the same in both modes and equal to 0.65:0.34. The same T Ca :T Na ratio (0.65:0.34) is found at a higher frequency f = 10 Hz (α = 0.5), when i pulse = 0.6i lim . The calculations show that the T Ca :T Na ratio depends only on the i av value and does not depend on the frequency and duty cycle; it does not matter if PEF or CC mode is used.

PEF Mode (Voltage Pulses)
Typical time dependences of the voltage and current density for the PEF mode, where constant voltage pulses alternate with the pauses during which the current is zero, are shown in Figure 6.
The time dependences of the partial current densities of Na + and Ca 2+ ions at a frequency of 0.5 Hz are shown in Figure 7. During a pulse lapse, the partial current density of Ca 2+ gradually decreases, while that of Na + increases. Qualitatively, this behavior is similar to that occurring in PEF mode, when constant current pulses are applied. Similar to above, the reason for the Ca 2+ current diminution is a more rapid decrease in its concentration compared to that of Na + in the depleted diffusion layer and membrane. During a pause lapse, the partial current density of Ca 2+ increases, while that of Na + decreases so that iNa becomes negative. It should be noted that the total current density during pulse lapse is significantly higher than that at CC mode at the same average value of the voltage, Uav.
Similar results were obtained at the frequency of PEF of 10 Hz (Figure 8). However, due to a shorter duration of the pulse, the Na + current density during a pulse lapse becomes greater than the current density of Ca 2+ at a higher value of the average voltage ( Figure 8c) than in the case of 0.5 Hz (Figure 7b). During a pulse lapse, the partial current density of Ca 2+ gradually decreases, while that of Na + increases. Qualitatively, this behavior is similar to that occurring in PEF mode, when constant current pulses are applied. Similar to above, the reason for the Ca 2+ current diminution is a more rapid decrease in its concentration compared to that of Na + in the depleted diffusion layer and membrane. During a pause lapse, the partial current density of Ca 2+ increases, while that of Na + decreases so that i Na becomes negative. It should be noted that the total current density during pulse lapse is significantly higher than that at CC mode at the same average value of the voltage, U av .
Similar results were obtained at the frequency of PEF of 10 Hz (Figure 8). However, due to a shorter duration of the pulse, the Na + current density during a pulse lapse becomes greater than the current density of Ca 2+ at a higher value of the average voltage (Figure 8c) than in the case of 0.5 Hz (Figure 7b). During a pulse lapse, the partial current density of Ca 2+ gradually decreases, while that of Na + increases. Qualitatively, this behavior is similar to that occurring in PEF mode, when constant current pulses are applied. Similar to above, the reason for the Ca 2+ current diminution is a more rapid decrease in its concentration compared to that of Na + in the depleted diffusion layer and membrane. During a pause lapse, the partial current density of Ca 2+ increases, while that of Na + decreases so that iNa becomes negative. It should be noted that the total current density during pulse lapse is significantly higher than that at CC mode at the same average value of the voltage, Uav.
Similar results were obtained at the frequency of PEF of 10 Hz (Figure 8). However, due to a shorter duration of the pulse, the Na + current density during a pulse lapse becomes greater than the current density of Ca 2+ at a higher value of the average voltage ( Figure 8c) than in the case of 0.5 Hz (Figure 7b). The calculations show ( Table 3) that the average current density, iav, in PEF mode is higher than that in CC mode at a given average voltage, Uav; moreover, iav increases with increasing frequency that agrees with the results of Sistat et al. [16]. This increase is explained by the fact that only the first moments (a few hundredths of a second) are profitable for obtaining a high current density (Figure 6b). With increasing duration of the pulse, the current density decreases rapidly due to increasing concentration polarization. However, when comparing the TCa:TNa ratios obtained in PEF mode and CC mode, it is found that this ratio is independent of how the PEF mode is applied, whether current or voltage pulses are applied, and what are the frequency and duty cycle. This ratio is a function of only one quantity, the average current density. To summarize, note that our simulation results are consistent with the experiment of Lemay et al. [17] regarding the comparison of CC and PEF modes, provided that the current density during the pulse is the same as in the CC mode. In this case, experiment and calculation show that the membrane permselectivity for Ca 2+ is improved in PEF mode. However, there is controversy about the effect of the frequency. The simulation shows that the membrane specific permselectivity in PEF mode does not depend on the frequency; only the average current density is important. However, the experiment indicates that under conditions that iav is the same, the permselectivity for Ca 2+ in the case of f = 0.5 Hz is greater than in the case of f = 5 Hz. Evidently, the mechanism of the effects observed in [17] is more complicated than the developed model suggests. It is possible that an important role in these phenomena is played by electroconvection. As mentioned in the introduction, at overlimiting average currents, electroconvection can increase the mass transfer rate by about 33% [25,26]. Electroconvective micro-vortices could essentially improve mixing the depleted solution near the membrane surface. The enhancement of convective mass transfer from the bulk to the membrane surface could shift the Ca 2+ :Na + concentration ratio in favor of Ca 2+ . However, modeling competitive ion transport with taking The calculations show ( Table 3) that the average current density, i av , in PEF mode is higher than that in CC mode at a given average voltage, U av ; moreover, i av increases with increasing frequency that agrees with the results of Sistat et al. [16]. This increase is explained by the fact that only the first moments (a few hundredths of a second) are profitable for obtaining a high current density (Figure 6b). With increasing duration of the pulse, the current density decreases rapidly due to increasing concentration polarization. However, when comparing the T Ca :T Na ratios obtained in PEF mode and CC mode, it is found that this ratio is independent of how the PEF mode is applied, whether current or voltage pulses are applied, and what are the frequency and duty cycle. This ratio is a function of only one quantity, the average current density. Table 3. Results of calculation of the average current density and transport numbers in the membrane in CC and PEF modes at different average voltages, U av , and different frequencies in Hz: f = 0 (CC mode), 0.5 and 10; α = 0.5.  To summarize, note that our simulation results are consistent with the experiment of Lemay et al. [17] regarding the comparison of CC and PEF modes, provided that the current density during the pulse is the same as in the CC mode. In this case, experiment and calculation show that the membrane permselectivity for Ca 2+ is improved in PEF mode. However, there is controversy about the effect of the frequency. The simulation shows that the membrane specific permselectivity in PEF mode does not depend on the frequency; only the average current density is important. However, the experiment indicates that under conditions that i av is the same, the permselectivity for Ca 2+ in the case of f = 0.5 Hz is greater than in the case of f = 5 Hz. Evidently, the mechanism of the effects observed in [17] is more complicated than the developed model suggests. It is possible that an important role in these phenomena is played by electroconvection. As mentioned in the introduction, at overlimiting average currents, electroconvection can increase the mass transfer rate by about 33% [25,26]. Electroconvective micro-vortices could essentially improve mixing the depleted solution near the membrane surface. The enhancement of convective mass transfer from the bulk to the membrane surface could shift the Ca 2+ :Na + concentration ratio in favor of Ca 2+ . However, modeling competitive ion transport with taking into account electroconvection is not an easy task; it can be the subject of another publication.

Conclusions
A mathematical model based on the Nernst-Planck and Poisson equations is reported to describe competitive counterion transport through an ion-exchange membrane. CC mode and PEF mode are compared. Two regimes of PEF mode are studied: when constant current pulses or constant voltage pulses alternate with pauses during which the current is zero. It is shown that when the same current density is applied in CC mode and during the pulses in PEF mode, the permselectivity of the membrane for a specific counterion is higher in the case of PEF mode. This result is in agreement with the experimental data obtained by Lemay et al. [17] in the ED demineralization of sweet whey. However, when comparing the results of the simulation, provided that the current density in CC mode is the same as the average current density in PEF mode, it is found that effective transport numbers of two competing ions (T 1 and T 2 ) in PEF mode are equal to their values in CC mode. Moreover, the values of T 1 and T 2 do not depend on whether current or voltage is applied during the pulses; they are also independent of the frequency and duty cycle used in PEF mode. The last results do not agree with the experiment of Lemay et al. [17]. It was found in [17], that the membrane permselectivity for Ca 2+ in the case of f = 0.5 Hz is greater than in the case of f = 5 Hz. Therefore, taking into account only electromigration and diffusion as the mechanisms of ion transport seems insufficient. It is possible that this discrepancy is due to electroconvection. Electroconvective micro-vortices could deliver a fresh solution from the solution bulk to the depleted membrane interface, which can maintain relatively high Ca 2+ :Na + concentration and flux ratios at this interface.
As for the practical application of PEF mode in ED, this seems very promising. The main effect of PEF, proved by numerous authors, is the mitigation of the membrane scaling and fouling. In addition, the use of PEF mode leads to an increase in current efficiency and mass transfer and to a decrease in the membrane stack resistance. However, generally, since the current is zero during the relaxation pauses, the ED process in this mode takes more time than in the classical CC mode. Therefore, the optimization of PEF parameters (frequency, duty cycle) is needed. Some of these problems are discussed in the paper by Martí-Calatayud et al. [51].
The effect of PEF on the membrane permselectivity for specific ions is still not clear enough. The model presented in this paper provides an insight into the understanding of this effect. However, more efforts should be applied to explore the role of electroconvection and perhaps some other effects.

Conflicts of Interest:
The authors declare no conflict of interest. Charge of fixed ions in the membrane Greek symbols α Duty cycle γ k Activity coefficient of ion k in the solution γ k Activity coefficient of ion k in the membrane δ I ,δ I I Thickness of the depleted and enriches diffusion boundary layers, respectively ε Solution relative permittivity ε 0 Vacuum permittivity θ k Ratio of the degree of desalination by ion k to the total degree of desalination ϕ Electric potential