Analysis of Membrane Transport Equations for Reverse Electrodialysis (RED) Using Irreversible Thermodynamics

Reverse electrodialysis (RED) is an electro-membrane process for the conversion of mixing energy into electricity. One important problem researchers’ face when modeling the RED process is the choice of the proper membrane transport equations. In this study, using experimental data that describe the membrane Nafion 120 in contact with NaCl aqueous solutions, the linear transport equation of irreversible thermodynamics was applied to calculate the power density of the RED system. Various simplifying assumptions about transport equation (i.e., four-, three-, and two-coefficients approaches) are proposed and discussed. We found that the two-coefficients approach, using the membrane conductivity and the apparent transport number of ions, describes the power density with good accuracy. In addition, the influence of the membrane thickness and the concentration polarization on the power density is also demonstrated.


Introduction
Global population growth contributes to the escalation of energy demand for households, transportation, and industry. The use of fossil fuels results in various environmental threats. Therefore, the search is increasingly widespread for alternative methods of energy conversion, e.g., osmotic into electrical [1] and mechanical energy (so called "osmotic pump" [2,3], or the conversion of electrokinetic energy [4]. Reverse electrodialysis (RED) is an electro-membrane process for the conversion of mixing energy into electricity. The concept of RED was invented in 1950 [5,6]; however, it was not heavily investigated until the 21st century, when this technique regained the attention of researchers looking for novel, unconventional, and renewable energy sources [7][8][9].
One of the important RED issues is the cost-effectiveness of the process. Turek and Bandura [10] suggested that energy conversion using RED is not profitable; however, Giacalone et al. [11] presented a much more optimistic perspective. There are a number of parameters that influence the efficiency of RED, like ionic shortcut currents [12], type of electrodes [13], kind of membrane [14], membrane thickness [15], fouling issues [16,17], efficiency of solutions mixing in the stack [18,19], and geometry of spacers [20].
In this section, we present a description of membrane transport in RED. In the literature discussing the modeling of RED, apart from very simplified membrane transport (i.e., considering only the transport numbers of ions [21]), usually the Nernst-Planck equation or its transformations are used [22][23][24][25][26][27][28]. In some research, the volume flow caused by the osmotic and electroosmotic transport of water through the membrane was considered [25,29]. Whereas in other research, only osmotic transport was considered [22,24]. In other studies, no volume transport was anticipated [27]. The problem of membrane transport is probably neglected because of the many other important aspects of RED. A general (i.e., not exactly related to RED) review of transport through ion-exchange membranes can be found elsewhere [30].
The Nernst-Planck equation is a simplified version of the linear transport equation of irreversible thermodynamics. The transport equation is comprised of six transport coefficients for the system membrane in contact with a single electrolyte aqueous solution, under the condition that the Onsager symmetry of the coefficients is fulfilled [31]. The determination of all these coefficients requires extensive work especially because their concentration dependence should be considered [32][33][34][35]. Consequently, the practical use of the full transport equation is limited.
Therefore, the aim of this work was to discuss which of the simplifications of the irreversible thermodynamics transport equation (ITE) is useful for calculating the power density in RED. The following three approaches based on four, three, and two coefficients (to be determined experimentally) were analyzed and discussed in detail: (1) the four-coefficients approach-the pressure difference in the full set of ITE is neglected, (2) the three-coefficients approach-the coupling of ion flows with the gradient of chemical potential of water is neglected, (3) the two-coefficients approach (i.e., Nernst-Planck equation)-the ion flux is a function of the gradient of its electrochemical potential only.
The results of these simplifications (i.e., values of power density vs. electric current density) are compared with the values obtained from the full set of ITE. The transport coefficients used in these latter equations were determined for the cation-exchange membrane Nafion 120 in contact with NaCl solutions [34,35] in a wide range of concentrations (0.05-3.7 M). To consider the concentration dependence, the differential form of the transport equation was used. This method was applied for the first time by McCallum and Meares to describe transport through an ion-exchange membrane [36]. The Nafion 120 membrane is relatively thick (0.3 mm); therefore, the transport coefficients are rescaled to smaller thickness values to demonstrate the influence of that parameter on the power density.

Theory
In the following, one basic part of the RED module is considered, i.e., one ion-exchange membrane with adjacent concentration polarization layers ( Figure 1). The steady-state conditions in the system are assumed, i.e., the fluxes of ions in the membrane and both polarization layers are the same.

Transport Equations of Irreversible Thermodynamics (ITE)
The starting point for the ITE approach to isothermal conditions is the dissipation function, Φ, the product of entropy production, and temperature, which can be expressed as the sum of driving forces, X i , and coupled with fluxes, J i : For discontinuous systems, the forces are represented by the differences in chemical potentials of transported species. In continuous systems, the negative values of gradients of these potentials are used.
According to irreversible thermodynamics, it is assumed that fluxes and forces are linearly related: and the Onsager reciprocity relation L ik = L ki holds if the forces are not too high. One can choose many sets of fluxes and forces; all of them have to fulfill the condition of the entropy production invariance: From this condition, we obtain various forms of the fluxes and forces, among which is the set of practical equations useful for experimental determination of transport coefficients [37]: where J 1 is the flux of ions 1 (cations), ν 1 is the number of ions 1 in the electrolyte molecule, s is volume flux, I is electric current density, ∆p is the pressure difference, ∆π is the osmotic pressure difference, c s is the mean concentration of electrolyte solutions contacting the membrane (c s , c s "), and L αβ is the mean transport coefficient for that concentration range. ∆E is the electric potential difference measured with electrodes reversible to anions (2): The power density is calculated from the electric potential difference ∆ϕ, determined from Equation (4) assuming that ∆ ln a 2 ≈ ∆ ln a ± .
The determination of all L αβ coefficients is a laborious task. In Table 1, we list the experiments we performed to determine L αβ in Equation (3).
Having these coefficients determined for the vanishingly small concentration differences (∆π, ∆µ s → 0) allowed us to assume the Onsager reciprocity relation (L αβ = L βα ) and determine the concentration dependence of six L αβ (α, β = π, p, E) [34,35]. The method of such determination was described in detail elsewhere [32]. The finite concentration differences (L αβ ) in Equation (3) are not symmetric because of their strong concentration dependence [32,38]. the real transport number of ion 1 the apparent transport number of ion 1 (electromotive force) Calculating the ion flux, volume flux, and electric potential difference, considering the concentration dependence of L αβ is the most exact approach. To calculate these quantities for the assumed conditions (concentrations on both sides of membrane, electric current density, membrane thickness), we applied the method proposed by McCallum and Meares [36]. According to this method, the membrane is divided into n m slices and for each slice the transport is given by the set of Equation (3) with local variables (electrolyte 1:1): where l m is the thickness of the real membrane for which the coefficients L αβ = f (c s ) were determined and w sl is the width of slice equal to the assumed membrane thickness (l m,assumed ) divided by the number of slices (w sl = l m,assumed /n m ). The concentration difference is expressed as ∆π and the formula ∆π = ∂π/∂c s ( c) c s,k−1 − c s,k . The fluxes across each slice are the same, assuming the stationary state in the considered system. Therefore, for n m slices, we have 3n equations and 3n m unknowns: c s,k , p k , E k , k = 1, . . . , n − 1, J 1 , and J v , which we solve for given values of parameters (electric current, concentration, pressure on both sides of membrane, and the assumed membrane thickness). All of these refer to the hypothetical (virtual) solution being in equilibrium with the membrane at coordinate x = kw sl , k = 1, . . . , n − 1; such a solution is defined by equality of chemical potentials of the transported species in that solution and in the membrane at a given x.
The concentration polarization layers adjacent to membrane surfaces were treated in a similar way as the membrane. The experimental data for NaCl solutions were taken from [39]. The used equations are provided in Appendix A.
A comment should be made on the extended Nernst-Planck equation, frequently used in the literature for describing the membrane transport of ions. Written in terms of hypothetical solution parameters, it takes the following form: where K p,i is the partition coefficient of ion I and α i is the convective coupling coefficient [40,41]. There are 2 coefficients per ion (K p,i α i is treated as one coefficient); thus, for two ions and J v described by 3 coefficients (L pπ , L pp , L pE ), there are 7 coefficients-more than in the set of Equation (3) assuming the symmetry of coefficients. For this reason, the extended Nernst-Planck equation will not be discussed here.

The Four-Coefficients Approach
In this case, the pressure difference across each slice is assumed to be zero. Although on both sides of membrane the pressure is the same, it may change across the membrane to maintain the volume flux across each slice the same. As the number of unknowns is reduced, the equation for J v in Equation (5) can be neglected (J 1 is needed in the concentration polarization layers). Thus Equation (5) is reduced to: The four coefficients (L ππ , L πE , L Eπ , L EE ) can be determined from P s,dif , t 1 , k m , and t 1,app -the Electromotive Force (EMF) measurements. However, practically, it is better to determine the electroosmotic coefficient W (Table 1) and to calculate t 1 from W and t 1,app .

The Three-Coefficients Approach
In the equation of ion transport, the coupling with water, i.e., the term L i0 X 0 is neglected. Thus, Equation (1) with L ii replaced by P i c i /RT reduces to: where X i is: We assume that L 21 = L 12 ; therefore, there are 3 parameters that can be obtained from (1) specific conductivity, k m , (2) electrolyte permeability P s , and (3) the apparent transport number, t 1,app . Instead of t 1,app , we could choose the real transport number of cation, t 1 ; however, because t 1,app includes the electroosmotic number of water and is much easier to determine than t 1 , we choose t 1,app . Notably, t 1,app and t 1 obtained from Equation (7) are the same because the flux of water is neglected here (see formulae for t 1 and t 1,app in Table 1). The EMF measurements allow for the determination of the apparent transport number of ion 1. Using Equation (7), we obtain the following set of practical equations (for the 1:1 electrolyte): where P 1 is given by: and P s by:

The Two-Coefficients Approach
In this approach, we assume that the flux of ion I depends only on the gradient of electrochemical potential of that ion: The following practical equations are derived from Equation (10): There are 2 transport coefficients only, P 1 and P 2 , with four experimentally available parameters: (1) specific conductivity k m , (2) transport number of cation t 1 , (3) electrolyte permeability P s , and (4) t 1,app . The relation between these quantities and P 1 , P 2 is as follows (assuming 1:1 electrolyte, and denoting ion 1 as a cation): and Similar to the 3-coefficients approach, the expressions for t 1 and t 1,app are the same because the water transport is also neglected.
The determination of t 1,app by the EMF method is much easier and therefore much more popular than the determination of the real transport number, t 1 . We chose the following pairs of experiments to determine P 1 and P 2 : (1) k m , t 1,app , (2) k m , P s , (3) t 1,app , P s . Additionally, contrary to t 1 , the value of t 1,app includes the electroosmotic transport of water; thus, to some extent, the transport of water molecules is considered.

Results and Discussion
As mentioned above, in this work, only one membrane (cation-exchange membrane) with adjacent concentration polarization layers is discussed.
In Figure 2, a typical dependence of the power density P on the current density I, calculated as P = −I∆ϕ, based on the full transport equation (Equation (5)), is shown for two concentration differences (0.1-1.0 M and 0.05-3.7 M) and two values of membrane thickness (0.1 and 0.3 mm). Two common trends were observed: the higher the concentrations difference, the higher the value of P; the lower the membrane thickness, the higher the value of P.
The four-coefficients approach, neglecting the pressure changes (Equation (6)), yields similar results. For 0.1-1.0 M NaCl, the difference in P max is below 1% irrespective of the membrane thickness. For higher concentration differences (0.05-3.7 M), this difference increases to ca. 5-6%.  [36]. They suggested that the membrane thickness could decrease in that case. At higher I values, the pressure becomes positive. This can be explained by different changes in the profiles of electroosmotic and osmotic parts of J v across the membrane caused by the change of concentration profiles (Figure 4). To keep J v constant, the pressure part of J v should change accordingly, which demands various profiles of pressure. Assuming the constant pressure value (the four-coefficients approach), the volume flux across the membrane is not constant-it can deviate from the actual steady state value even several times, which is unrealistic.  The influence of the membrane thickness on the maximum of P and the current density at P max is shown in Figure 5. On both curves (P max = f (l m ) and −I(P max ) = f (l m )), maxima are observed at l m < 0.1 mm, and the maximum on P max is at a slightly higher l m value than that of −I(P max ). The higher concentration difference demands a thicker membrane (maxima are shifted toward thicker membranes). The thinner the membrane, the higher I, but also the higher the permeability of electrolytes; the impact of diffusion layers (of the assumed constant thickness) also becomes relatively stronger. The dependence of the maximum of power density P max and the current density at P max on the diffusion boundary layer thickness l p is a monotonically decreasing function ( Figure 6). Due to a substantial decrease in P max with l p , the reduction of the thickness of that layer by ensuring adequate stirring of the solutions at the membrane surface is an important practical issue. Even more important is the decrease in l p on the lower concentration side of membrane. To reduce the concentration polarization effect, an extra energy for pumping is needed. Figure 6. Dependence of the maximum of power density, P max , and the current density, −I, at P max on the diffusion boundary layer thickness, l p , for 0.1-1.0 M and 0.05-3.7 M NaCl, T = 298 K, Equation (5), full ITE.

The Two-Coefficients Approach
The results for the two-coefficients approach, for which P 1 and P 2 are calculated either from t 1,app and k m (TC), or from t 1,app and P s (TP), or from k m and P s (CP), are shown in Figure 7. Among these three combinations of coefficients, the smallest deviation from the power density calculated from the full ITE was observed for the TC combination (t 1,app , k m ). It yielded even better results than the three-coefficients approach based on all three parameters (TCP: t 1,app , k m and P s ). Power density vs. current density for full model and approximate approaches: TCP-3-coefficients approach (P 1 , P 2 , and L12 calculated from t 1,app , k m , and P s ), (t 1app , P s ); TP, CP, TC-2-coefficients approach (P 1 and P 2 calculated from t 1,app and P s , k m and P s , t 1,app and k m , respectively); the Nafion 120 membrane separates (a) 0.1 and 1 M NaCl, (b) 0.05 and 3.7 M NaCl; l m = 0.1 mm, l p = 0.05 mm, T = 298 K.
Because of the quite good fit of the TC approach, another question arose: is it necessary to determine the concentration dependence of t 1,app and k m in the whole range c -c", or is it enough to determine t 1,app in that range and k m for the mean concentration of c and c"?
In that case (P 1 , P 2 = const), we can determine the integrated form for J 1 from Equation (10): ∆E can be expressed as a function of a certain mean concentration of c and c": P 1 and P 2 are determined from k m (Equation (12)), with c s replaced by c s , and t 1,app calculated from (15) with dE/dµ s replaced by (E(c") − E(c ))/(µ s (c") − µ s (c )). To check this simplification, the calculations were performed for c s being the arithmetic and logarithmic means of c , c". The boundary diffusion layers were treated as in the previous cases to keep the same conditions for all discussed approaches. To simplify calculations, in the calculations of P 1 , P 2 , and c s , the external bulk concentrations were used, not those at the membrane surfaces. A comparison of this approach (TC const.) with the full ITE equation and with the TC approach is shown in Figure 8. This simplification of the logarithmic mean of concentration (marked in Figure 8 as "TC const ln") yields comparable results to the other approaches for the lower concentration difference (0.1-1.0 M). For 0.05-3.7 M the predictions of "TC const ln" are too low at the current density range 0-j(P max ).
In this case, the concentration dependence of t 1,app and k m is needed to improve the predictions of power density. The arithmetic mean of concentration yields much worse results especially for high concentration difference and is not shown in Figure 8.

Conclusions
Among the tested simplifications of the full set of transport equations of irreversible thermodynamics (six independent coefficients), the four-coefficients approach (k m , t 1 , t 1,app , P s,dif ) with only the pressure effect neglected yields almost the same results for 0.1-1 M as the full ITE approach.
The three-coefficients approach (k m , t 1,app , P s,dif ), where only ion fluxes are considered with their couplings, slightly overestimates the power density at lower concentration difference. However, at higher concentration differences, it strongly overestimates and becomes much less accurate than the two-coefficients approach (equation of the Nernst-Planck type, assuming no coupling between ion fluxes) based on k m and t 1,app . Other combinations of the two-coefficients approach are substantially worse. The model based on k m and t 1,app can be further simplified by using constants P 1 and P 2 calculated from the mean values of k m and t 1,app in a given concentration range c -c" and using the logarithmic mean of these concentrations. These coefficients are easily accessible experimentally.
None of the discussed approximations considered the volume flow effect on the concentrations of the solutions contacting the membrane. To do so, we need to estimate the osmotic volume permeability and the electroosmotic coefficient W.
To obtain the optimal power density for given concentrations of feed solutions, the appropriate membrane thickness should be chosen. The higher the concentration difference, the thicker the membrane should be. As the position of P max depends on the membrane properties, the optimal thickness should be independently determined for cation-and anion-exchange membranes.
The thickness of the diffusion boundary layer should be as low as possible especially on the low concentration side of the membrane-it strongly decreases the power density.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Transport Equation in the Concentration Polarization Layers
In [39], the transport coefficients, l ik , are determined for the set of equations: in the solvent-fixed reference frame. Equation (A1) can be transformed with forces similar to those in Equation (3) (electrolyte 1:1): where l 1E = F(l 11 − l 12 ), κ e = F 2 (l 11 − 2l 12 + l 22 ), and a = (1 − c s v s )/c s . Due to volume flow through the membrane, we assume that the flux of ion 1 in the boundary layer is given by: In the practical equations (Equations (3) and (5)), J v is defined as J v = v 0 J 0 + v s J 1 /ν 1 ; thus, for our case (electrolyte 1:1), Equation (A3) can be rewritten as: To simultaneously calculate the c and E profiles in the boundary layers with the membrane, we divide both layers into n e slices for which the transport equations are as follows: J 1 =(l 11 ( c s )a( c s )∆π + L 1E ( c s )∆E)(1 − c s v s )/w sl,e + c s J v I =(L 1E ( c s )a( c s )∆π + κ e ( c s )∆E)/w sl,e (A4) where w sl,e is the width of the slices in those layers. Dividing each layer into n p slices we have (4n p + 3n m ) equations (including equations for the membrane) and the same number of unknowns: (2n p + n m − 1) of c s,k , (2n p + n m ) of E k , (n m − 1) of p k , and J 1 , J v .