Accurate Determination of Electrical Potential on Ion Exchange Membranes in Reverse Electrodialysis

Reverse electrodialysis is a promising membrane technology to generate energy from controlled mixing of water streams of different salinities. Electrical potentials generate on the ion exchange membranes (IEMs) when selective transport of cations and anions across the membranes driven by concentration difference. The accurate determination of the potentials developed on the IEMs is critical to fairly assess the feasibility of the technology. The Nernst–Planck–Poisson (NPP) equations for IEMs (the membranes with fixed charge) were solved numerically with the boundary updating scheme. The validity of this numerical method was verified by the identical values of Donnan potential obtained with the well-established analytical methods. The suitability and applicability of the classic Teorell–Meyer–Siever (TMS) model were assessed by comparison to the simulation results from the numerical method.


Introduction
Salinity gradient energy (SGE) is a type of clean and renewable energy embodied in water streams of different salinities [1,2]. The amount of SGE is estimated at about 2.6 terawatts globally or approximately 20% of the current worldwide energy consumption rate [3]. Pressure retarded osmosis (PRO) and reverse electrodialysis (RED) are currently under investigation as two promising technologies to harvest SGE [3][4][5][6]. In PRO process, a semipermeable membrane is employed to separate two solutions of different salt concentration. Driven by osmotic pressure, water transports from low concentration side to high concentration side and builds up a hydraulic pressure there. The controlled release of the pressurized water from the high concentration side can generate energy by turning a hydro turbine. The working principle, process features, and technical and economic feasibilities of PRO have been extensively studied.
RED, which is the other promising technology that can be used to harvest SGE, is also a membrane-based process. Instead of semipermeable membranes, ion exchange membranes (IEMs) are used in RED processes [7][8][9][10][11]. IEMs are permselective membranes due to the fixed charges on membrane materials. There are two types of IEMs according to the fixed charges carried by the membranes. Cation exchange membranes (CEMs) carry negative fixed charges and only allow cations to pass while anion exchange membranes (AEMs) carry positive fixed charges and only allow anion to pass. When a CEM and an AEM are arranged on two sides of a concentrated flow, sandwiched with two diluted flows, cations and anions in the concentrated flow would move through membranes to opposite directions. As a result, electrical potentials develop on the membranes. Electrical energy can be produced when the potentials are applied to an external load.
Development of potential on IEMs is a very complex phenomenon that is affected by many factors, such as properties of IEMs, diffusivities of various ions, and configuration of the RED cells. Although it is well established that membrane potential is governed by the Poisson equation coupled with the Nernst-Planck equation for concentration distributions of ions along the membrane thickness, seeking a rigorous solution of the Nernst-Planck-Poisson (NPP) equations remains a serious challenge in the field. The commonly used Nernst equation [5,6] and Teorell-Meyer-Siever (TMS) model [12][13][14] for the potential on IEMs are at best approximations or surrogates for the true membrane potential because both models are obtained with suspicious simplifications and assumptions [15,16].
An optional method for the membrane potential is to directly seek numerical solutions of the NPP equations without the commonly used simplifications and assumptions [17][18][19][20]. However, new challenge emerges that the potentials and ion concentrations at the membrane surfaces cannot be assigned a priori in the rigorous formulation of the problem [21,22]. Inappropriately imposed boundary conditions were often used in the numerical solutions in previous studies until we recently developed a boundary updating scheme to handle this challenge successfully [23]. With this scheme, the NPP equations for a membrane of no-fixed charge was rigorously solved for the first time under the true boundary conditions.
The main objective of this article is to apply the boundary updating scheme to determine the potential developed on the membranes with fixed charges. Numerical solutions of the NPP equations for IEMs with various charge density were obtained and the suitability and accuracy of the numerical method were discussed. The impact of charge density of IEMs on the total potential and potential components were studied. Finally, the classic TMS model was compared and assessed for its applicability and limitations with the numerical method.

Governing Equations
The problem under consideration is to determine the membrane potential at steady state on an IEM as schematically shown in Figure 1. The IEM of fixed charge density X separates two solutions of a monovalent salt with concentrations C b0 and C bL , respectively. There is a transition layer between the membrane surface and the bulk solution on each side of the membrane. While electroneutrality is always maintained in the bulk solutions, charge is usually unbalanced in the transition layers because the attraction of counterions and repulsion of co-ions by the charge on the membrane. As a result, both ion concentrations and potential on the membrane surfaces differ from those in the bulk solutions, as schematically described by the lines for concentrations and potential in Figure 1. There is only a concentration profile for one ion in the figure for the sake of clarity. The concentration profile for the other ion would be different in the transition layers and on the membrane. of the RED cells. Although it is well established that membrane potential is governed by the Poisson equation coupled with the Nernst-Planck equation for concentration distributions of ions along the membrane thickness, seeking a rigorous solution of the Nernst-Planck-Poisson (NPP) equations remains a serious challenge in the field. The commonly used Nernst equation [5,6] and Teorell-Meyer-Siever (TMS) model [12][13][14] for the potential on IEMs are at best approximations or surrogates for the true membrane potential because both models are obtained with suspicious simplifications and assumptions [15,16].
An optional method for the membrane potential is to directly seek numerical solutions of the NPP equations without the commonly used simplifications and assumptions [17][18][19][20]. However, new challenge emerges that the potentials and ion concentrations at the membrane surfaces cannot be assigned a priori in the rigorous formulation of the problem [21,22]. Inappropriately imposed boundary conditions were often used in the numerical solutions in previous studies until we recently developed a boundary updating scheme to handle this challenge successfully [23]. With this scheme, the NPP equations for a membrane of no-fixed charge was rigorously solved for the first time under the true boundary conditions.
The main objective of this article is to apply the boundary updating scheme to determine the potential developed on the membranes with fixed charges. Numerical solutions of the NPP equations for IEMs with various charge density were obtained and the suitability and accuracy of the numerical method were discussed. The impact of charge density of IEMs on the total potential and potential components were studied. Finally, the classic TMS model was compared and assessed for its applicability and limitations with the numerical method.

Governing Equations
The problem under consideration is to determine the membrane potential at steady state on an IEM as schematically shown in Figure 1. The IEM of fixed charge density X separates two solutions of a monovalent salt with concentrations Cb0 and CbL, respectively. There is a transition layer between the membrane surface and the bulk solution on each side of the membrane. While electroneutrality is always maintained in the bulk solutions, charge is usually unbalanced in the transition layers because the attraction of counter-ions and repulsion of co-ions by the charge on the membrane. As a result, both ion concentrations and potential on the membrane surfaces differ from those in the bulk solutions, as schematically described by the lines for concentrations and potential in Figure 1. There is only a concentration profile for one ion in the figure for the sake of clarity. The concentration profile for the other ion would be different in the transition layers and on the membrane. The ion concentrations on a membrane are governed by the Nernst-Planck Equation (1): where C i is the concentration of ion i, D i is the diffusion coefficient of ion i, x is the spatial coordinate perpendicular to the membrane, z i is the valance of ion i, F is the Faraday constant, R is the gas constant, and ψ is the potential. Convection term was not considered in Equation (1) because no water flow across the membrane was assumed. The potential on the membrane is then governed by Poisson equation: where is the permittivity of the membrane material, and X is the fixed charge density. Please be reminded that the variable X can take a positive value for AEMs or a negative value for CEMs. The boundary conditions for the NPP equations in the problem depicted in Figure 1 cannot be specified a priori. There is a dilemma between the boundary conditions and the final solutions of the equations. The ion concentrations and potential on the membrane surfaces are strongly dependent on the free charge developed on the membrane. However, the free charge on the membrane can only be known after the equations have been solved with the specified boundary conditions. We will handle this dilemma with the boundary updating scheme described below.

Boundary Updating Scheme
The boundary updating scheme was initially developed in one of our previous papers [23]. For convenience of readers, a more concise description of the method is provided here. The scheme starts with a state of the membrane for which the boundary conditions can be easily specified. For example, there is no potential on a neutral membrane (referring to the bulk solutions) and the ion concentrations at the membrane surfaces are equal to the bulk solutions, i.e.: where C i0 and C iL are the concentrations of ion i on the left surface and right surface of the membrane, respectively, C ib0 and C ibL are the concentrations of ion i in the bulk solutions on the left and right sides of the membrane, respectively, and ψ 0 is the potential on the left surface of the membrane, and  where J i is the flux of ion i. The derivatives of concentration and potential at two membrane surfaces are calculated from the numerical solutions. The cumulative net charges in the two transition layers can be calculated with: and where Q 0 and Q L are the cumulative charges in the transient layers on left and right sides of the membrane, respectively, Q 0 and Q L are the cumulative charges at the previous time step in the conresping transient layers, respectively, and ∆t is timestep size. At the first timestep, Q 0 = 0 and Q L = 0 for the initially neutral membrane.
The boundary conditions at any timestep are updated with the cumulative charges in the transition layers by: where λ 0 and λ L are the Debye lengths in the solutions on the left and right sides of the membrane, which are calculated by: The numerical solution of the NPP equations can be obtained for the next timestep with the updated boundary conditions. The above procedure is repeated with the newly obtained numerical solution to update the boundary conditions for the next timestep until the steady state is reached, which is indicated by the null current condition: With boundary updating scheme, the concentration boundary conditions are always consistent with potential boundary conditions at any time step. Therefore, the true solution of NPP equations at the steady state is guaranteed.

Calculation of Membrane Potentials
With the numerical solution at the steady state, the ion fluxes can be calculated by: The total membrane potential can be calculated with: where ψ T is the total potential of the membrane, ∆ψ m is the potential difference across the membrane thickness, and ∆ψ 0 and ∆ψ L are the potential differences across the transition layers on the two sides of the membrane, respectively. The potential difference across the membrane thickness is available from the numerical solution of the Poisson equation. The two potential differences across the transition layers can be calculated from the cumulated charges in the layers with:

Simulations and Discussion
A C++ program for the numerical procedure described above was developed on Visual Studio 2019 for calculation of the potential on IEMs under various conditions. All simulations reported below were done at a PC with CPU of Intel i7-9700 at 3.00 Ghz, on which a numerical solution can be obtained in a few seconds. Unless other specified, the parameters listed in Table 1 were used in the simulations.

The Effectiveness of Boundary Updating Schemne
The numerical solution of NPP equations for parameters given in Table 1 is presented in Figure 2. The membrane in this case is an AEM because of the positive fixed charge on it. The numerical solution for a membrane of no-fixed-charge for the same parameters is also shown in the figure. Figure 2a shows that concentrations of both cations and anions at the membrane surfaces (boundaries) differ significantly from the bulk concentrations, which are 10 mol/m 3 and 5 mol/m 3 , respectively, in the solutions on the left and right sides of the membrane. The different cation and anion concentrations at the membrane surfaces are an essential feature of ion transport across the membranes, which is especially true for IEMs. The boundary updating scheme is particularly effective for the IEMs because of the bigger differences between cation and anion concentrations at the surfaces than those for the membranes of no-fixed-charge.
The anion concentration is higher than the cation concentration throughout the entire membrane thickness because of the positive fixed charge on the AEM. The electroneutrality can only be roughly maintained in a middle section of the membrane thickness, i.e., the difference between anion concentration and cation concentration is nearly equal to the fixed positive charge of the AEM. The lower anion concentrations and higher cation concentrations (unbalanced charges) in the regions near the membrane surfaces are the results of ion transfer between the membrane and solutions. The distributions of ions for the membrane of no-fixed-charge are relatively simple. Concentrations of anions and cations are equal through almost the entire membrane thickness with exceptions in two narrow regions near the membrane surfaces. The concentration differences cannot be clearly seen for the membrane of no-fixed-charge in the figure. Interested readers can find more details about the distributions of ion concentrations on the membrane of no-fixed-charge in one of our previous papers [23].  The corresponding potential distributions are presented in Figure 2b. The potential on the membrane of no-fixed-charge decreases monotonically from the left surface to the right surface. The potential on the AEM behaves more complexly than that on the membrane of no-fixed-charge. The middle section of potential on the AEM declines with a rate like that on the membrane of no-fixed-charge. This behavior of potential is obviously due to the differential diffusion of ions with different mobilities. A jump in an initial region and a drop in a final region of the potential on the AEM can be attributed to the repulsion of cations and the attraction of anions by the fixed positive charge on the AEM, i.e., the Donnan effect of the fixed charge on the membrane. The potentials at the membrane surfaces were again determined by the boundary updating scheme, which are consistent with the ion concentrations at the membrane surfaces. Such rigorous numerical solutions of NPP equations with the well-defined consistent boundary conditions have not been obtained before with any methods other than the boundary updating scheme.

Impact of the Charge Density of IEMs
Potential distributions along the membrane thickness for IEMs with different charge densities are presented in Figure 3. The thick line in the middle presents the potential profile on the membrane of no-fixed-charge. The potentials for AEMs are all above the line for the uncharged membrane while those for CEMs are below. This feature can be attributed to the Donnan effect of the fixed charge on the IEMs. The parallel middle sections of all the potential lines decline from the left to the right reflect impact of the differential diffusion of cations and anions on the potential profiles.
Another useful feature of the boundary updating scheme is that the potential differences across the transition layers are determined simultaneously with the numerical solution of the NPP equations. The potential differences across the two transition layers and across the membrane calculated with parameters in Table 1 are presented in Figure 4, together with the total membrane potential. The potential difference on the membrane Δϕm varies in a range of ~5mV under the simulation conditions. The potential difference across the left transition layer Δϕ0 varies in a range that is more than doubled the range for the The corresponding potential distributions are presented in Figure 2b. The potential on the membrane of no-fixed-charge decreases monotonically from the left surface to the right surface. The potential on the AEM behaves more complexly than that on the membrane of no-fixed-charge. The middle section of potential on the AEM declines with a rate like that on the membrane of no-fixed-charge. This behavior of potential is obviously due to the differential diffusion of ions with different mobilities. A jump in an initial region and a drop in a final region of the potential on the AEM can be attributed to the repulsion of cations and the attraction of anions by the fixed positive charge on the AEM, i.e., the Donnan effect of the fixed charge on the membrane. The potentials at the membrane surfaces were again determined by the boundary updating scheme, which are consistent with the ion concentrations at the membrane surfaces. Such rigorous numerical solutions of NPP equations with the well-defined consistent boundary conditions have not been obtained before with any methods other than the boundary updating scheme.

Impact of the Charge Density of IEMs
Potential distributions along the membrane thickness for IEMs with different charge densities are presented in Figure 3. The thick line in the middle presents the potential profile on the membrane of no-fixed-charge. The potentials for AEMs are all above the line for the uncharged membrane while those for CEMs are below. This feature can be attributed to the Donnan effect of the fixed charge on the IEMs. The parallel middle sections of all the potential lines decline from the left to the right reflect impact of the differential diffusion of cations and anions on the potential profiles.
Another useful feature of the boundary updating scheme is that the potential differences across the transition layers are determined simultaneously with the numerical solution of the NPP equations. The potential differences across the two transition layers and across the membrane calculated with parameters in Table 1 are presented in Figure 4, together with the total membrane potential. The potential difference on the membrane ∆φ m varies in a range of~5 mV under the simulation conditions. The potential difference across the left transition layer ∆φ 0 varies in a range that is more than doubled the range for the membrane potential (~10 mV). The potential difference across the right transition layer ∆φ L varies in a range of about 25 mV. Under the simulation conditions, the potential differences in the transition layers are primarily controlled by Donnan potential. The larger absolute value of the potential difference in the right transition layer is reasonable because Donnan potential of the same fixed charge density is greater in the solution of lower concentration. The total potential difference of the membrane ∆φ T decreases monotonically throughout the entire range of charge density. throughout the entire range of charge density.
Ion fluxes at various charge densities are listed in Table 2 with the corresponding membrane potentials and total potentials. The impact of the charge density on flux is much more moderate than on the potentials. Ion flux changes less than 50% in magnitude for the range of charge density used in the simulations. It clearly demonstrates that the ion flux through the membrane is largely determined by the concentration difference across the membrane with relatively minor impact of charge density on the membrane. This observation is reasonable because concentration difference across the membrane is the original driving force for ion transport. The fixed charge density on the membrane can only modify the ratio of cation and anion concentrations inside the membrane.    Ion flux increases as the charge density decreases in CEMs and maximizes at charge density of −5 mol/m 3 . Ion flux then declines as charge density further decrease to zero in CEMs and as charge density increase in AEMs. In this example, because cations have Ion fluxes at various charge densities are listed in Table 2 with the corresponding membrane potentials and total potentials. The impact of the charge density on flux is much more moderate than on the potentials. Ion flux changes less than 50% in magnitude for the range of charge density used in the simulations. It clearly demonstrates that the ion flux through the membrane is largely determined by the concentration difference across the membrane with relatively minor impact of charge density on the membrane. This observation is reasonable because concentration difference across the membrane is the original driving force for ion transport. The fixed charge density on the membrane can only modify the ratio of cation and anion concentrations inside the membrane. Ion flux increases as the charge density decreases in CEMs and maximizes at charge density of −5 mol/m 3 . Ion flux then declines as charge density further decrease to zero in CEMs and as charge density increase in AEMs. In this example, because cations have smaller mobility than the anions, the appropriate increase in cation concentration would improve ion flux due to the zero current condition at steady state. The positive fixed charge in the AEMs will reduce cation concentration and therefore further reduce the capacity of ion transport. On the other hand, the negative fixed charge of the membrane can increase cation concentration in the CEMs so that the capacity of ion transport can be increased. Of course, the anion concentration would be reduced in the CEMs at the same time. Because of the larger mobility of anions, the impact of the reduced anion concentration in a certain range on ion flux would be over-compensated by the benefit of the increased cation concentration. Beyond this range, anion would become the limiting factor for ion transport. Further decrease in anion concentration would reduce the ion flux.

Comparison of Donnan Potentials Determinded Analytically and Numerically
As a demonstration of the boundary updating scheme, Donnan potential was calculated by solving the NPP equations with solutions of equal concentrations on both sides of an AEM. Concentrations C b0 = C bL = 10 mol/m 3 and X = 10 mol/m 3 were used in the simulation. The distributions of ion concentrations and potential at the steady state from a numerical solution is presented in Figure 5. Anion concentration is higher than cation concentration throughout the membrane thickness. In a middle section of membrane thickness (~10-40 nm), electroneutrality is roughly satisfied because difference between anions and cations is about equal to the fixed charge density on the membrane. The results also show that there are some ion exchanges between the initially neutral membrane and solutions. As a result, lower anion concentration and higher cation concentration compared to those in the middle section were observed in regions near the membrane surfaces. The unbalanced charges in these regions are the cause for the Donnan potential. The boundary updating scheme is capable to catch this subtle feature of the numerical solution for the NPP equations. Potential distribution is also presented in Figure 5 for the AEM. The potential increases inward from the membrane surfaces on the membrane. The value of the plateau on the potential curve is the Donnan potential, which is 12.36 mV in this case.
Donnan potential for a membrane of a fixed change density X in contact with solution of a monovalent salt can be calculated analytically by [14,24]: where ∆ψ D is the Donnan potential and C is the salt concentration in the solution. Equation (21) was derived with the electroneutrality assumption of the membrane. Donnan potentials determined both numerically with our method and analytically with Equation (21) are listed in Table 3. It is a good surprise to find that the two methods yield Donnan potentials that are identical up to 3 effective digits! It is convincing evidence that the numerical method for membrane potential is very accurate. At the same time, it also shows that the use of electroneutrality assumption to reach the analytical expression for Donnan potential is acceptable though it is unphysical fundamentally.
thickness (~10-40nm), electroneutrality is roughly satisfied because difference between anions and cations is about equal to the fixed charge density on the membrane. The results also show that there are some ion exchanges between the initially neutral membrane and solutions. As a result, lower anion concentration and higher cation concentration compared to those in the middle section were observed in regions near the membrane surfaces. The unbalanced charges in these regions are the cause for the Donnan potential. The boundary updating scheme is capable to catch this subtle feature of the numerical solution for the NPP equations. Potential distribution is also presented in Figure 5 for the AEM. The potential increases inward from the membrane surfaces on the membrane. The value of the plateau on the potential curve is the Donnan potential, which is 12.36 mV in this case. Donnan potential for a membrane of a fixed change density X in contact with solution of a monovalent salt can be calculated analytically by [14,24]: where ∆ is the Donnan potential and C is the salt concentration in the solution. Equation (21) was derived with the electroneutrality assumption of the membrane. Donnan potentials determined both numerically with our method and analytically with Equation (21) are listed in Table 3. It is a good surprise to find that the two methods yield Donnan potentials that are identical up to 3 effective digits! It is convincing evidence that the numerical method for membrane potential is very accurate. At the same time, it also shows that the use of electroneutrality assumption to reach the analytical expression for Donnan potential is acceptable though it is unphysical fundamentally.

The Applicability and Accuracy of TMS Model
Membrane potential on an IEM with monovalent salt solutions of different concentration on two sides of the membrane can be calculated analytically with the classic TMS model [14]: where U is the mobility coefficient that is defined as: where D + and D − are the diffusion coefficients of cations and anions, respectively. The potentials on an AEM calculated analytically with Equation (22) for varying charge density are presented (as dots) in Figure 6. In the calculation, concentrations C b0 = 10 mol/m 3 and C bL = 1 mol/m 3 were employed on two sides of the membrane. The numerically calculated membrane potentials under the same conditions are also presented (as lines) in Figure 6.
The mobility coefficients are indicated on the graphs. The results for a CEM would produce similar graphs as Figure 6 but with the mobility coefficients flipped on the figure.
where and are the diffusion coefficients of cations and anions, respectively. The potentials on an AEM calculated analytically with Equation (22) for varying charge density are presented (as dots) in Figure 6. In the calculation, concentrations Cb0 = 10 mol/m 3 and CbL = 1 mol/m 3 were employed on two sides of the membrane. The numerically calculated membrane potentials under the same conditions are also presented (as lines) in Figure 6. The mobility coefficients are indicated on the graphs. The results for a CEM would produce similar graphs as Figure 6 but with the mobility coefficients flipped on the figure. Figure 6. Membrane potentials on an AEM calculated analytically (symbols) and numerically (lines) as a function of ratio of the fixed charge density to salt concentration on the right side of the membrane. The mobility coefficients are indicated on the graphs. Figure 6 shows that the membrane potentials from TMS model agree well with the numerically calculated values when the ratio of X/Cb0 is smaller than 0.03 or greater than 30. However, there are significant discrepancies between the two values in the range of 0.03 < X/Cb0 < 30. The TMS potentials are systematically lower than the true membrane potentials, with the biggest discrepancy at X/Cb0 = 1. For example, the TMS potential Figure 6. Membrane potentials on an AEM calculated analytically (symbols) and numerically (lines) as a function of ratio of the fixed charge density to salt concentration on the right side of the membrane. The mobility coefficients are indicated on the graphs. Figure 6 shows that the membrane potentials from TMS model agree well with the numerically calculated values when the ratio of X/C b0 is smaller than 0.03 or greater than 30. However, there are significant discrepancies between the two values in the range of 0.03 < X/C b0 < 30. The TMS potentials are systematically lower than the true membrane potentials, with the biggest discrepancy at X/C b0 = 1. For example, the TMS potential (8.71 mV) is only about 40% of the numerical potential (22.0 mV) for mobility coefficient of 0.9 at X/C b0 = 1.
TMS model is not a rigorous solution of the NPP equations because it is simply impossible. The assumptions and simplifications [12,13] used to reach the analytical solution would certainly cause some deviations from the true solution of the equations. On the other hand, the numerical potentials can be taken as the true values because no assumptions and simplifications are made in the solution seeking process. Furthermore, the TMS model can be only used for limited cases where only two types of ions of equal charges exist in solutions while the numerical method can be applied to the general cases with multiple ions of various charges.

Conclusions
The potential of an IEM in RED can be rigorously calculated with the boundary updating scheme. The validity and accuracy of the numerical method is strongly supported by the identical Donnan potentials calculated with the well-established analytical method and the numerical method presented in this study. It is demonstrated that the classic TMS model tends to underestimate the membrane potentials systematically for IEMs with fixed charge density in the intermediate range.