Trade-Off between Operating Time and Energy Consumption in Pulsed Electric Field Electrodialysis: A Comprehensive Simulation Study

Electrodialysis (ED) has been recently introduced in a variety of processes where the recovery of valuable resources is needed; thus, enabling sustainable production routes for a circular economy. However, new applications of ED require optimized operating modes ensuring low energy consumptions. The application of pulsed electric field (PEF) electrodialysis has been demonstrated to be an effective option to modulate concentration polarization and reduce energy consumption in ED systems, but the savings in energy are usually attained by extending the operating time. In the present work, we conduct a comprehensive simulation study about the effects of PEF signal parameters on the time and energy consumption associated with ED processes. Ion transport of NaCl solutions through homogeneous cation-exchange membranes is simulated using a 1-D model solved by a finite-difference method. Increasing the pulse frequency up to a threshold value is effective in reducing the specific energy consumption, with threshold frequencies increasing with the applied current density. Varying the duty cycle causes opposed effects in the time and energy usage needed for a given ED operation. More interestingly, a new mode of PEF functions with the application of low values of current during the relaxation phases has been investigated. This novel PEF strategy has been demonstrated to simultaneously improve the time and the specific energy consumption of ED processes.


Introduction
Electrodialysis (ED) is a mature separation process used to obtain concentrated and diluted streams from saline solutions [1]. In ED, an electric field is applied to set ions in motion within an electrochemical cell, where the use of ion-exchange membranes as separators between different cell compartments creates the concentrate and dilute streams. In addition to the most traditional use of ED as a desalination technology, it has recently been used for new applications in processes where valuable ionic species need to be separated or concentrated from a feed solution [2][3][4][5][6]. ED is a modular technology, where temperature changes such as those needed in evaporative processes are avoided. Also, the addition of chemicals is minimal as compared to competing technologies, such as in those involving the use of ion-exchange resins or chemical precipitation. Thus, the successful production of highly selective and conductive ion-exchange membranes has brought the expansion of electrodialysis into food and biotechnological sectors [7][8][9].
Yet, energy consumption continues to be one of the main challenges of the process. A considerable part of the energy usage is associated with the application of an electric field and the evolving cell electrical resistance. Within an electrochemical cell, the solution layers next to the membranes are the parts where the formation of concentration gradients is more acute; hence, where the major changes in electrical resistance are located. The development of concentration gradients is termed concentration polarization in the field of membrane processes [10]. In ED systems, concentration polarization is an intrinsic consequence of membrane selectivity that depends on the level of applied current density. Apart from the increased cell resistance, in some specific membrane-electrolyte systems, this phenomenon can be accompanied by additional detrimental effects, such as the occurrence of undesired pH changes or the formation of precipitates on the membrane surface [11,12]. Consequently, ED has been commonly limited to working at moderate current densities.
A major research line aiming at improving ED processes is the application of new modes of operation. Enhancement of electroconvection by operating at overlimiting currents is widely studied. When high current densities are applied, electrical instabilities emerge in the depleted solution layers and promote an additional supply of ions from the bulk solution toward the membrane surface [13,14]. Another choice is the operation of ED systems with pulsed electric fields (PEF), which alternate periods of current application with pauses, where the intensity of concentration polarization decreases. Several research works have demonstrated that applying PEF modes can lead to decreased energy consumption and mitigation of fouling phenomena [15][16][17][18]. It has also been found that a combination of PEFs and overlimiting currents can induce increased desalination rates, if electroosmotic mixing lasts during the pauses [19,20].
The choice of the parameters of the pulse functions has a paramount influence on the efficiency of this operation mode because formation and relaxation of concentration profiles near ion-exchange membranes is a dynamic process. Previous experimental works have shown that the duration of the pulse and pause periods under different current density regimes can lead to significant differences in terms of the number of charges transported, energy consumption or fouling mitigation [21][22][23][24]. However, the empirical investigation and selection of PEF parameters can result in a time-and resource-consuming practice. The use of flexible but rather simple models to simulate and predict the outcome of applying different pulse shapes can be extremely advantageous in this optimization process. In the present work, we conduct a comprehensive investigation of the effects of different pulse parameters on the energy and operating time needed to obtain a given charge transport in ED systems. A finite difference method applied to a simplified ion transport model is used to evaluate the interplay between the duration and level of applied current during the pulse and pause periods, when PEFs are implemented in a system formed by a cation-exchange membrane and NaCl solutions. In our opinion, the power of this strategy resides not only on a more efficient optimization procedure, but also on the access to information otherwise difficult to obtain experimentally.

System Geometry and Characteristics
A membrane-solution system consisting of a homogeneous cation-exchange membrane surrounded by a NaCl solution was simulated with a one-dimensional model using a finite difference method. The geometry of the system is shown schematically in Figure 1, where five differentiated regions can be distinguished:
The total thickness of the simulated system corresponds to 1 mm of the ED cell. The thickness (∆) of each diffusion boundary layer (DBLs) in the simulation was 250 µm, which corresponds to a typical value for ED systems [25][26][27]. A membrane thickness (d) of 50 µm was considered [28], while the remaining 450 µm correspond to the portion of well-stirred bulk solution taken into account in the computation of the system voltage drop. The electrolyte considered in the present work is 0.05 M NaCl solutions. During the pulse periods, the system operates under galvanostatic mode. By working with current density values (A/m 2 ), the results obtained are specific for a membrane area of 1 m 2 . Operation at room temperature of 25 • C was assumed. It must be noted that the results obtained for Membranes 2021, 11, 43 3 of 15 a membrane system would be analogous for any membrane within the central repeating units of an ED stack. Consequently, this study will be focused on the ion transport taking place through a single ion-exchange membrane. drop. The electrolyte considered in the present work is 0.05 M NaCl solutions. During the pulse periods, the system operates under galvanostatic mode. By working with current density values (A/m 2 ), the results obtained are specific for a membrane area of 1 m 2 . Operation at room temperature of 25 °C was assumed. It must be noted that the results obtained for a membrane system would be analogous for any membrane within the central repeating units of an ED stack. Consequently, this study will be focused on the ion transport taking place through a single ion-exchange membrane. Figure 1. Scheme of the system simulated formed by a cation-exchange membrane surrounded by two diffusion boundary layers and portions of bulk solution. Concentration profiles are represented in the y-axis versus the spatial coordinate x for the initial and steady-state conditions.

Model
The flux of ionic species in the system was modeled using an approximation to the Nernst-Planck equation for transport of charged particles: where Ji(x) represents the total flux of the species i at a given position x in the spatial coordinate. The first term on the right-hand side of the equation represents the flux of the species i caused by diffusion, the middle term represents the flux by migration and the third term corresponds to the transport in the x direction caused by convection. In Equation (1), Di, Ci(x), and zi represent the diffusion coefficient, the concentration, and the valence of ion i; while F, R and T represent the Faraday constant, the gas constant and the temperature. φ(x) and v(x) refer to the electric potential and the velocity at x.
In the bulk solution, far enough from the membrane surface, fluid flow implies that the convective term provides sufficient mixing, so that the concentration can be considered homogeneous in the x direction, as shown schematically in Figure 1. On the contrary, in the DBLs (i.e., very thin solution layers at both sides of the membrane surface), the contribution of convection to ion transport can be neglected. Under the application of an electric field, the selectivity of the membrane implies that transport by migration of counterions in the membrane is higher than in the solution, so that concentration gradients are established at both sides of the membrane. Considering a symmetric 1:1 salt (NaCl) and a galvanostatic mode of operation, Equation (1) that models the flux of ions i at any position x within the DBLs (Ji,DBL(x)) can be converted into: Figure 1. Scheme of the system simulated formed by a cation-exchange membrane surrounded by two diffusion boundary layers and portions of bulk solution. Concentration profiles are represented in the y-axis versus the spatial coordinate x for the initial and steady-state conditions.

Model
The flux of ionic species in the system was modeled using an approximation to the Nernst-Planck equation for transport of charged particles: where J i (x) represents the total flux of the species i at a given position x in the spatial coordinate. The first term on the right-hand side of the equation represents the flux of the species i caused by diffusion, the middle term represents the flux by migration and the third term corresponds to the transport in the x direction caused by convection. In Equation (1), D i , C i (x), and z i represent the diffusion coefficient, the concentration, and the valence of ion i; while F, R and T represent the Faraday constant, the gas constant and the temperature. ϕ(x) and v(x) refer to the electric potential and the velocity at x.
In the bulk solution, far enough from the membrane surface, fluid flow implies that the convective term provides sufficient mixing, so that the concentration can be considered homogeneous in the x direction, as shown schematically in Figure 1. On the contrary, in the DBLs (i.e., very thin solution layers at both sides of the membrane surface), the contribution of convection to ion transport can be neglected. Under the application of an electric field, the selectivity of the membrane implies that transport by migration of counterions in the membrane is higher than in the solution, so that concentration gradients are established at both sides of the membrane. Considering a symmetric 1:1 salt (NaCl) and a galvanostatic mode of operation, Equation (1) that models the flux of ions i at any position x within the DBLs (J i,DBL (x)) can be converted into: In Equation (2), the convective term is neglected, the coupled transport of anions and cations is considered by using the salt diffusion coefficient (D salt ), i op refers to the operating The model used in the present study considers the ion concentration and ion flux as the variables of the system, while the rest of the terms are assumed constant.
The membrane is regarded as a dense phase and has a high concentration of fixed charges, so that migration of ions is the only mechanism of ion transport considered in this part of the system. Accordingly, the flux of ions through the membranes (J i,m (x)) can be expressed as a function of the operating current density: In Equation (3), T i refers to the transport number of ions i in the membrane phase. To obtain the time evolution of the concentration in the two DBLs, the continuity equation was also considered: The system of equations is solved for a system formed by a 0.05 NaCl solution and a cation-exchange membrane, with transport number values of Na + and Cl − ions in the membrane phase of T Na + = 0.95 and T Cl − = 0.05, respectively. The transport number of Na + and Cl − ions in the solution were obtained from the equivalent conductivity data at infinite dilution, resulting in t Na + = 0.396 and t Cl − = 0.604 [29]. The salt diffusion coefficient (D salt = 1.611·10 −9 m 2 /s) was also calculated from the equivalent ionic conductivity data at infinite dilution using the Nernst-Einstein equation [30]. The simulation problem is solved for different functions of operating current considering partitions in time (t) and space (x) coordinates. Initial and boundary conditions are needed to solve the simulation problem. The initial conditions (t = 0) assumed in our work involve the constant concentration of NaCl in both DBLs equal to C bulk = 0.05 M and an initial flux of ions equal to 0: The system of equations is solved in each DBL starting from the space node located at the membrane solution interface until reaching the boundary with the bulk solution at a distance ∆ from the membrane. The boundary conditions at the edges between the DBL and the bulk solution imply that the concentration is equal to C bulk , while the boundary conditions at both membrane-DBL interfaces are represented by the equality of fluxes at both sides of the interface: Equations (2) and (4) solved in this model are sensitive in terms of stability, so they were converted into dimensionless forms by dividing the variables c i , J i , t i and x i by reference values cº i , Jº i , tº i and xº i , according to the rules of Langtangen and Pedersen, where the subindex i refers to either Na + or Cl − ions [31]. The system is solved for one spatial dimension, x, using the finite difference method, which allows obtaining the profiles of the dependent variables, J i and C i , after their conversion in the dimensional forms. The system of equations was applied to each DBL subsystem, while the subsystem membrane and bulk solution can be assumed as constant ohmic resistances. The space and time variables are discretized establishing grids over both domains at space and time steps, ∆x and ∆t, respectively.
Once the concentration profiles across the system at different times have been obtained, the potential difference between two consecutive spatial nodes (i−1 and i) can be calculated by using the following equation [27]: which in its integral form for the investigated system can be expressed as: The function k(C NaCl ), which represents the conductivity of NaCl solutions as a function of the salt concentration, was obtained from the interpolation of experimental data reported in the literature [29]. ∆U(i − 1 : i, t) values are computed between all adjacent nodes of the DBLs and are summed to obtain the voltage drop for the diluted and concentrated DBL subsystems, ∆U dil (t) and ∆U conc (t), respectively. The voltage drop in the portions where a constant resistance is assumed, i.e., the solution bulk and the membrane phase, can be calculated according to Equation (12): where L represents the thickness of the subsystem and k subsystem was calculated as k(0.05 M) with the interpolated function k(C NaCl ) for the bulk solution, and a characteristic value was taken for the membrane conductivity. The total voltage drop of the membrane system ∆U m was calculated as the sum of the contributions of all subsystems: Finally, the specific energy consumption (kJ/m 2 ) for a given degree of desalination or number of charges transported can be calculated by integrating the product between the voltage drop through the membrane system and the input function of operating current density over time: Here, we want to emphasize that this model entails several assumptions and simplifications. For instance, it is assumed that the local electroneutrality condition holds in all parts of the system, and the dissociation of water has not been taken into account. These phenomena may become relevant, especially at overlimiting currents, but are very specific of the membrane and solution properties. Also, the use of 1D models implies neglecting the effects of membrane heterogeneity and electroconvection. However, the agreement between the ohmic regions of current-voltage curves obtained using 1D and 2D models was confirmed by M.K. Urtenov et al. in a previous work [32].

Implementation of PEF Currents
The calculation of the dynamic system response with the model described above was accomplished for varying PEF input functions of current density, i op (t). As can be seen in Figure 2, the input functions alternate periods where a constant current density is applied with periods of pauses or relaxation of the membrane system, where a zero or low level of current density is applied. Thus, four primary parameters of the input functions were varied in this work: the level of current density applied during the pulses (i high ), the level of current density applied during the pauses or relaxation periods (i low ), the time of pulse duration (t on ) and the time of pause or relaxation duration (t off ). The application of current pulses was also compared with the behavior of the system under galvanostatic mode without pauses.
Membranes 2021, 11, x FOR PEER REVIEW 6 of 15 duration (ton) and the time of pause or relaxation duration (toff). The application of current pulses was also compared with the behavior of the system under galvanostatic mode without pauses. The variation of the primary pulse function parameters (ihigh, ilow, ton, toff) was implemented to simulate the system under varying values of frequency and duty cycle, which are parameters more commonly used in the literature to define the characteristics of PEF modes of operation. The frequency of a cycle (f) is defined as: The duty cycle (α) correlates the time during which the pulses are applied with the total duration of a cycle: In the following section, a comprehensive investigation of the effect of the frequency and duty cycle of the pulses on the membrane voltage drop at different levels of pulse current density (ihigh) is presented. The values of current during the relaxation periods (ilow) were also investigated. It must be noted that the simulation time was adapted in order to compare the results for a value of transported charge equal to 16,800 A·s in all cases.

Results Obtained Applying a Constant Current in Continuous Mode
First, the transport of ions through the membrane system was simulated under the application of a continuous current at different levels, so that these reference conditions could be used for comparison purposes with the PEF modes. Considering the characteristics of the simulated system described in the previous section, the limiting current density can be calculated theoretically by means of the Peers' equation [33]: The current density for the membrane-electrolyte system selected for this study corresponds to a value of 56.12 A/m 2 . Table 1 summarizes the parameters considered for the simulations in continuous galvanostatic mode. Apart from the operation at ilim, two levels of underlimiting current density were also simulated to evaluate the system response at varying degrees of concentration polarization. For evaluation purposes, the application of overlimiting currents was also investigated. Here, it must be noted that at overlimiting The variation of the primary pulse function parameters (i high , i low , t on , t off ) was implemented to simulate the system under varying values of frequency and duty cycle, which are parameters more commonly used in the literature to define the characteristics of PEF modes of operation. The frequency of a cycle (f ) is defined as: The duty cycle (α) correlates the time during which the pulses are applied with the total duration of a cycle: In the following section, a comprehensive investigation of the effect of the frequency and duty cycle of the pulses on the membrane voltage drop at different levels of pulse current density (i high ) is presented. The values of current during the relaxation periods (i low ) were also investigated. It must be noted that the simulation time was adapted in order to compare the results for a value of transported charge equal to 16,800 A·s in all cases.

Results Obtained Applying a Constant Current in Continuous Mode
First, the transport of ions through the membrane system was simulated under the application of a continuous current at different levels, so that these reference conditions could be used for comparison purposes with the PEF modes. Considering the characteristics of the simulated system described in the previous section, the limiting current density can be calculated theoretically by means of the Peers' equation [33]: The current density for the membrane-electrolyte system selected for this study corresponds to a value of 56.12 A/m 2 . Table 1 summarizes the parameters considered for the simulations in continuous galvanostatic mode. Apart from the operation at i lim , two levels of underlimiting current density were also simulated to evaluate the system response at varying degrees of concentration polarization. For evaluation purposes, the application of overlimiting currents was also investigated. Here, it must be noted that at overlimiting currents, mechanisms of mass transfer different from those contemplated in the model come into play. More specifically, the dissociation of water and the onset of electroconvection become relevant when the electrolyte concentration at the membrane-DBL interface approaches zero [34][35][36]. Thus, the results obtained at i op > i lim will be considered only for short time spans during which the concentration of species at the depleting membrane-DBL interface is still far from zero. The simulation time was chosen to be 10 min for the 0.5·i lim case, while it was adapted for the rest of simulations so as to simulate the same amount of charge transported through the membrane.  Figure 3 shows the concentration profiles resulting from the simulation results under the continuous application of current. Note that electroneutrality conditions across the DBLs involve that the concentration profiles for Cl − ions have the same shape as those shown in the figures for Na + . Figure 3a shows the evolution with time of the concentration in the two DBLs when a current density equal to i lim is applied. It can be seen how the concentration of Na + ions decreases gradually with time in the depleting DBL until reaching a value close to zero at the membrane interface after 60 s of current application. On the contrary, in the concentrate compartment, the concentration rises as a consequence of the faster transport of Na + through the membrane as compared with the transport through the DBL towards the cathode. Figure 3b shows the steady-state concentration profile obtained under the application of different levels of current density. As expected, when the steady state is achieved, the fully developed concentration profiles across the Nernst diffusion layer adopt the shape of straight lines. While the Na + concentration at the depleting membrane-DBL interface takes the value of zero for i lim , this value is significantly higher in the case of underlimiting currents.   The evolution with time of the Na + concentration at the depleting membrane-DBL interface for the system under the application of different values of current density is represented in Figure 4a. Here, it can also be seen that the concentration at the membrane surface reaches a value of zero when i op = i lim and the steady state is achieved. This result confirms that the simulation is consistent with the Peers' equation. Figure 4b also shows the evolution of the voltage drop across the membrane system for the three values of applied current. Here, the initial jump in U m registered at the beginning of the simulation is approximately proportional to the applied current, as it is expected from a quasi-ohmic relationship between current and voltage when the influence of concentration polarization is still small. However, as the time increases, a jump in the U m value is registered before reaching a stationary value, especially for the i lim case. The time needed to reach the steady state also increases with the applied current: the value of U m reaches a constant value at around 50 s when 0.5·i lim is applied, while it needs more than 100 s in the case of i lim . From the results obtained with the continuous application of current, we can already expect that the implementation of current pulses will be more beneficial at high current densities, when the concentration profiles in the DBLs are more abrupt and the associated increase in U m is bigger.

Effect of Frequency
The effect of the frequency of PEF input functions was analyzed at three different levels of current density: 0.7·ilim, ilim and 1.2·ilim. As deduced from the previous subsection, the almost flat evolution of Um with time for 0.5·ilim indicates that the benefits of applying pulse current functions may be minor at very low values of current density. In contrast, the results obtained at 1.2·ilim are expected to be more beneficial in terms of energy savings; although the model can only be used in the unsteady condition for this current, i.e., under the condition that the electrolyte concentration at the depleting membrane-DBL is higher than zero. The duration of the pulses was also selected from the evolution of Um versus

Effect of Frequency
The effect of the frequency of PEF input functions was analyzed at three different levels of current density: 0.7·i lim , i lim and 1.2·i lim . As deduced from the previous subsection, the almost flat evolution of U m with time for 0.5·i lim indicates that the benefits of applying pulse current functions may be minor at very low values of current density. In contrast, the results obtained at 1.2·i lim are expected to be more beneficial in terms of energy savings; although the model can only be used in the unsteady condition for this current, i.e., under the condition that the electrolyte concentration at the depleting membrane-DBL is higher than zero. The duration of the pulses was also selected from the evolution of U m versus time shown in Figure 4b, considering times before the attainment of the steady state. In this section, the pause duration (t off ) was set to the same value as the pulse duration (t on ), so the values of α are 0.5 for all simulations. The pulse function parameters are summarized in Table 2.  Figure 5a shows some characteristic concentration profiles at times corresponding to the pulse and pause periods. When the current is on, the concentration drops through the diluting DBL and increases through the concentrating DBL. In both cases the concentration decreases from the anodic side of the electrochemical cell (left part of the graphs) towards the cathodic side (right part of the graphs). The gradient is also sharper near the membrane-DBL interfaces. However, when the current ceases in a pause phase, the relaxation of the concentration gradients is also more notorious near the membrane interfaces. This effect can be related to the driving force for diffusive mass transfer existing in the system right after the application of a current pulse: the differences of concentration are bigger next to the membrane surface, and the diffusive mass transfer between two locations where the differences in concentration are bigger takes place faster. The concentration of Na + at the diluting membrane-DBL interface is also plotted at different frequencies for i op = i lim in Figure 5b. In contrast with the results achieved when the current is continuously applied, the application of pulses promotes an oscillatory evolution of the concentration with time. The most remarkable outcome from the application of pulses is that the concentration does not reach the value of 0 at the membrane interface. Thus, concentration polarization can be effectively attenuated. For a duty cycle of α = 0.5, the mean value of concentration at the membrane interface is approximately invariant with the frequency; however, the amplitude of the oscillations becomes smaller as the frequency increases. time shown in Figure 4b, considering times before the attainment of the steady state. In this section, the pause duration (toff) was set to the same value as the pulse duration (ton), so the values of α are 0.5 for all simulations. The pulse function parameters are summarized in Table 2.  Figure 5a shows some characteristic concentration profiles at times corresponding to the pulse and pause periods. When the current is on, the concentration drops through the diluting DBL and increases through the concentrating DBL. In both cases the concentration decreases from the anodic side of the electrochemical cell (left part of the graphs) towards the cathodic side (right part of the graphs). The gradient is also sharper near the membrane-DBL interfaces. However, when the current ceases in a pause phase, the relaxation of the concentration gradients is also more notorious near the membrane interfaces. This effect can be related to the driving force for diffusive mass transfer existing in the system right after the application of a current pulse: the differences of concentration are bigger next to the membrane surface, and the diffusive mass transfer between two locations where the differences in concentration are bigger takes place faster. The concentration of Na + at the diluting membrane-DBL interface is also plotted at different frequencies for iop = ilim in Figure 5b. In contrast with the results achieved when the current is continuously applied, the application of pulses promotes an oscillatory evolution of the concentration with time. The most remarkable outcome from the application of pulses is that the concentration does not reach the value of 0 at the membrane interface. Thus, concentration polarization can be effectively attenuated. For a duty cycle of α = 0.5, the mean value of concentration at the membrane interface is approximately invariant with the frequency; however, the amplitude of the oscillations becomes smaller as the frequency increases. A more detailed examination of the evolution of U m for different frequencies is shown in Figure 6. When the pulse frequency is low (Figure 6a), the duration of the pulse is longer, and U m reaches higher values during the pulse phase. However, the longer duration of the relaxation phases at lower frequencies implies the restoration of U m to lower values before the start of a new pulse. Conversely, at high frequencies (Figure 6b), concentration profiles undergo a relaxation process before reaching critical U m values, so that concentration polarization can be controlled more efficiently. relaxation phases; (b) evolution with time of the NaCl concentration at the interface between the membrane surface and the diluting diffusion boundary layer (DBL) (position x = 0) under the application of ilim in continuous mode and in PEF mode at different frequencies. The arrows in the first plot indicate the evolution of concentration profiles during the pulses and relaxation phases.
A more detailed examination of the evolution of Um for different frequencies is shown in Figure 6. When the pulse frequency is low (Figure 6a), the duration of the pulse is longer, and Um reaches higher values during the pulse phase. However, the longer duration of the relaxation phases at lower frequencies implies the restoration of Um to lower values before the start of a new pulse. Conversely, at high frequencies (Figure 6b), concentration profiles undergo a relaxation process before reaching critical Um values, so that concentration polarization can be controlled more efficiently. To provide a global picture of the effect of pulse frequency on the specific energy consumption, the different simulated results are compared in Figure 7. The effect of the application of PEFs is very moderate at 0.7·ilim, where concentration polarization is not very intense. However, the reduction in energy consumption becomes noticeable when increasing the frequency from 0.01 Hz to a value of 0.02 Hz, and then reaches a plateau at higher frequencies. At current densities equal to ilim, when the operation is changed to PEF modes, the reduction in specific energy consumption represents more than a third with respect to the value obtained in continuous mode. The biggest change is observed when frequency increases from 0 to 0.01 Hz, and then energy consumption decreases more gradually with frequency until reaching a plateau around 0.05 Hz. Thus, for the investigated system, the application of frequencies higher than 0.05 Hz does not provide a further improvement of the process. For the operation at 1.2·ilim, similar trends confirm the effect of frequency on the specific energy consumption. However, the minimum value is shifted towards higher frequencies. To provide a global picture of the effect of pulse frequency on the specific energy consumption, the different simulated results are compared in Figure 7. The effect of the application of PEFs is very moderate at 0.7·i lim , where concentration polarization is not very intense. However, the reduction in energy consumption becomes noticeable when increasing the frequency from 0.01 Hz to a value of 0.02 Hz, and then reaches a plateau at higher frequencies. At current densities equal to i lim , when the operation is changed to PEF modes, the reduction in specific energy consumption represents more than a third with respect to the value obtained in continuous mode. The biggest change is observed when frequency increases from 0 to 0.01 Hz, and then energy consumption decreases more gradually with frequency until reaching a plateau around 0.05 Hz. Thus, for the investigated system, the application of frequencies higher than 0.05 Hz does not provide a further improvement of the process. For the operation at 1.2·i lim , similar trends confirm the effect of frequency on the specific energy consumption. However, the minimum value is shifted towards higher frequencies. Figure 7. Effect of frequency of PEF implementation on the specific energy consumption at different levels of current density and duty cycle α = 0.5. Empty symbols represent the operation in continuous mode.

Effect of Duty Cycles
The significant reduction in energy consumption observed with the application of PEFs at different frequencies could be further optimized by varying the duty cycle, that is, by implementing asymmetric pulse functions with differing ton and toff values. Moreover, the evaluation of the effect of duty cycle on the performance of ED systems may also show improvements, not only related to the energy consumption, but also to the time needed to transport a given number of ions through the membranes. Since in the previous section, the specific energy consumption reached an almost minimum plateau value at frequencies close to 0.05 Hz, this value has been taken as the central one around which the effect of duty cycle is investigated. It must be mentioned that, by varying the duty cycle, the values of frequency also change slightly around the central value of 0.05 Hz. Table 3 summarizes the conditions of current density pulse functions. In essence, in this section, the study is equivalent to varying the pause duration for a fixed pulse duration of 10 s. Some of the most representative results obtained by varying α are presented in Figure 8. For pulse periods with ton = 10 s, values of α higher than 0.5 imply pauses shorter than 10 s, while lower values of α imply extending the pause duration. Regarding the evolution of Na + concentration at the diluting membrane-DBL interface (Figure 8a), higher duty cycles lead to a more intense depletion of the DBL, with concentrations at x = 0 oscillating around 0.01 mol/L when ihigh = ilim. On the contrary, a lower α value of 0.44 leads to Na + concentrations at the membrane-solution interface oscillating around 0.025 mol/L. It can also be seen that the oscillation amplitudes become broader when α decreases, which is caused by the longer relaxation times leading to an enhanced restoration of the concentration profiles within the DBLs. Sosa-Fernandez et al. already demonstrated in a previous study, that extending the pause periods (toff values) during the desalination by PEF ED of

Effect of Duty Cycles
The significant reduction in energy consumption observed with the application of PEFs at different frequencies could be further optimized by varying the duty cycle, that is, by implementing asymmetric pulse functions with differing t on and t off values. Moreover, the evaluation of the effect of duty cycle on the performance of ED systems may also show improvements, not only related to the energy consumption, but also to the time needed to transport a given number of ions through the membranes. Since in the previous section, the specific energy consumption reached an almost minimum plateau value at frequencies close to 0.05 Hz, this value has been taken as the central one around which the effect of duty cycle is investigated. It must be mentioned that, by varying the duty cycle, the values of frequency also change slightly around the central value of 0.05 Hz. Table 3 summarizes the conditions of current density pulse functions. In essence, in this section, the study is equivalent to varying the pause duration for a fixed pulse duration of 10 s. Table 3. Pulse parameters considered for the investigation of the effect of duty cycle on the ED system. Some of the most representative results obtained by varying α are presented in Figure 8. For pulse periods with t on = 10 s, values of α higher than 0.5 imply pauses shorter than 10 s, while lower values of α imply extending the pause duration. Regarding the evolution of Na + concentration at the diluting membrane-DBL interface (Figure 8a), higher duty cycles lead to a more intense depletion of the DBL, with concentrations at x = 0 oscillating around 0.01 mol/L when i high = i lim . On the contrary, a lower α value of 0.44 leads to Na + concentrations at the membrane-solution interface oscillating around 0.025 mol/L. It can also be seen that the oscillation amplitudes become broader when α decreases, which is caused by the longer relaxation times leading to an enhanced restoration of the concentration profiles within the DBLs. Sosa-Fernandez et al. already demonstrated in a previous study, that extending the pause periods (t off values) during the desalination by PEF ED of polymer-flooding produced water yielded lower energy consumptions [16].
However, this reduction in energy consumption occurs at the cost of extending the time of operation significantly. Figure 8b shows the correlation between specific energy consumption and time needed to achieve a charge transfer of 16,800 A·s for different duty cycles. As found in previous works, decreasing α implies a decrease in the energy consumption, being this reduction more notorious at high current densities. On the contrary, for α values higher than 0.5, the times saved by decreasing the pause periods implied an increase in energy consumption. Thus, the application of asymmetric pulses can be used to obtain a trade-off between energy and time needed to desalinate a feed stream, being a useful strategy depending on the specific priorities of the given application of ED. polymer-flooding produced water yielded lower energy consumptions [16]. However, this reduction in energy consumption occurs at the cost of extending the time of operation significantly. Figure 8b shows the correlation between specific energy consumption and time needed to achieve a charge transfer of 16,800 A·s for different duty cycles. As found in previous works, decreasing α implies a decrease in the energy consumption, being this reduction more notorious at high current densities. On the contrary, for α values higher than 0.5, the times saved by decreasing the pause periods implied an increase in energy consumption. Thus, the application of asymmetric pulses can be used to obtain a tradeoff between energy and time needed to desalinate a feed stream, being a useful strategy depending on the specific priorities of the given application of ED.

Effect of the Current Density Applied during the Relaxation Phases
In addition to changing the frequency and duty cycle, which depend on the values of ton and toff, and evaluating their effect on the system's performance at different current density regimes (ihigh), an alternative means to reduce the time needed to achieve the transport of a given number of charges through ED membranes is to implement periods of relaxation at currents different from zero. To our knowledge, this type of PEF implementation has not been evaluated before, and it could simultaneously entail improvements in both the time and the energy needed to reach a specific desalination degree. To explore this new type of pulse functions, the simulation with symmetric pulses at frequencies equal to 0.05 Hz (ton = toff = 10 s) was performed with values of ilow of 0, 0.25·ilim and 0.5·ilim. Figure 9 shows an overview of the simulation results obtained varying ilow. In Figure 9a, the evolution with time of the concentration of Na + ions at the diluting membrane-DBL interface are shown for two different values of ilow. The evolution of the concentration with time clearly shows that applying a value of current density different from zero during the relaxation phases of the PEF functions is effective in attenuating concentration polarization, because the Na + concentration increases during this phase of the input signals. In fact, the shape of the curves is similar to the case where the current is turned off during

Effect of the Current Density Applied during the Relaxation Phases
In addition to changing the frequency and duty cycle, which depend on the values of t on and t off , and evaluating their effect on the system's performance at different current density regimes (i high ), an alternative means to reduce the time needed to achieve the transport of a given number of charges through ED membranes is to implement periods of relaxation at currents different from zero. To our knowledge, this type of PEF implementation has not been evaluated before, and it could simultaneously entail improvements in both the time and the energy needed to reach a specific desalination degree. To explore this new type of pulse functions, the simulation with symmetric pulses at frequencies equal to 0.05 Hz (t on = t off = 10 s) was performed with values of i low of 0, 0.25·i lim and 0.5·i lim . Figure 9 shows an overview of the simulation results obtained varying i low . In Figure 9a, the evolution with time of the concentration of Na + ions at the diluting membrane-DBL interface are shown for two different values of i low . The evolution of the concentration with time clearly shows that applying a value of current density different from zero during the relaxation phases of the PEF functions is effective in attenuating concentration polarization, because the Na + concentration increases during this phase of the input signals. In fact, the shape of the curves is similar to the case where the current is turned off during the pauses. While the oscillations take place at a slightly lower value of concentration, their amplitude is somewhat smaller for values of i low = 0.25·i lim . Moreover, during the relaxation phases, current is still being effectively transported through the membranes, which implies a reduction in the operating time. Figure 9b summarizes the results obtained in this section. The application of 0.25·i lim during the relaxation phases of the pulse functions induce a simultaneous reduction of the time and energy needed to achieve the same number of charges transported through the membranes, thus proving the convenience of adopting this strategy in ED systems. This improvement is general for all values of i high , but more significant for i high values close or higher than i lim . A further increase in the i low values up to 0.5·i lim shows to be also beneficial for the time needed for a given desalination objective, but it represents a further increase in the specific energy consumption as compared with the cases where i low = 0.25·i lim . These improvements observed in the ED performance by applying low values of current during the relaxation phase of PEF modes are even more notorious when compared with the results obtained operating in continuous mode. For example, for the case when i high = i lim , a substantial decrease in specific energy is achieved, changing from 4.53 kJ/m 2 in the continuous mode to 2.42 kJ/m 2 in the PEF mode with i low = 0.25·i lim . However, the operating time only increases from five to eight minutes, thus representing a significant improvement as compared with conventional PEF modes where the current ceases completely during the relaxation phases.
Membranes 2021, 11, x FOR PEER REVIEW 13 of 15 the pauses. While the oscillations take place at a slightly lower value of concentration, their amplitude is somewhat smaller for values of ilow = 0.25·ilim. Moreover, during the relaxation phases, current is still being effectively transported through the membranes, which implies a reduction in the operating time. Figure 9b summarizes the results obtained in this section. The application of 0.25·ilim during the relaxation phases of the pulse functions induce a simultaneous reduction of the time and energy needed to achieve the same number of charges transported through the membranes, thus proving the convenience of adopting this strategy in ED systems. This improvement is general for all values of ihigh, but more significant for ihigh values close or higher than ilim. A further increase in the ilow values up to 0.5·ilim shows to be also beneficial for the time needed for a given desalination objective, but it represents a further increase in the specific energy consumption as compared with the cases where ilow = 0.25·ilim. These improvements observed in the ED performance by applying low values of current during the relaxation phase of PEF modes are even more notorious when compared with the results obtained operating in continuous mode. For example, for the case when ihigh = ilim, a substantial decrease in specific energy is achieved, changing from 4.53 kJ/m 2 in the continuous mode to 2.42 kJ/m 2 in the PEF mode with ilow = 0.25·ilim. However, the operating time only increases from five to eight minutes, thus representing a significant improvement as compared with conventional PEF modes where the current ceases completely during the relaxation phases.

Conclusions
In this work, the effect of the input signal parameters on the performance of pulsed electric field electrodialysis has been systematically evaluated using a facile simulation model solved applying finite element methods. Different parameters, such as pulse frequency, duty cycle or applied current density in both the pulse and relaxation phases have been evaluated in terms of the specific energy consumption and operating time needed to achieve a constant transport of charges through the membranes. In general, it can be concluded that changing from continuous to PEF modes implies a general decrease in specific energy consumption; being this improvement more significant at high current densities, where high degrees of polarization are reached. There is a threshold frequency value from

Conclusions
In this work, the effect of the input signal parameters on the performance of pulsed electric field electrodialysis has been systematically evaluated using a facile simulation model solved applying finite element methods. Different parameters, such as pulse frequency, duty cycle or applied current density in both the pulse and relaxation phases have been evaluated in terms of the specific energy consumption and operating time needed to achieve a constant transport of charges through the membranes. In general, it can be concluded that changing from continuous to PEF modes implies a general decrease in specific energy consumption; being this improvement more significant at high current densities, where high degrees of polarization are reached. There is a threshold frequency value from which further improvements in energy consumption become almost negligible. This frequency depends on the level of applied current i high , taking a value of 0.02 Hz for 0.7·i lim and increasing up to 0.05 Hz for 1.2·i lim . Decreasing the duty cycle provides moderate savings in energy consumption, while increasing it implies significant reductions in the operating time needed for a given desalination process. For the first time, a novel type of PEF function with application of current during the relaxation periods has been evaluated. These new PEF modes are found to be a promising strategy to effectively reduce the specific energy consumption, without compromising the operating time of ED processes. Application of low levels of current density during the relaxation periods is useful for attenuating concentration polarization, while maintaining ion transport during the entire ED operating time.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.