Modelling of Ion Transport in Electromembrane Systems: Impacts of Membrane Bulk and Surface Heterogeneity

Artificial charged membranes, similar to the biological membranes, are self-assembled nanostructured materials constructed from macromolecules. The mutual interactions of parts of macromolecules leads to phase separation and appearance of microheterogeneities within the membrane bulk. On the other hand, these interactions also cause spontaneous microheterogeneity on the membrane surface, to which macroheterogeneous structures can be added at the stage of membrane fabrication. Membrane bulk and surface heterogeneity affect essentially the properties and membrane performance in the applications in the field of separation (water desalination, salt concentration, food processing and other), energy production (fuel cells, reverse electrodialysis), chlorine-alkaline electrolysis, medicine and other. We review the models describing ion transport in ion-exchange membranes and electromembrane systems with an emphasis on the role of microand macroheterogeneities in and on the membranes. Irreversible thermodynamics approach, “solution-diffusion” and “pore-flow” models, the multiphase models built within the effective-medium approach are examined as the tools for describing ion transport in the membranes. 2D and 3D models involving or not convective transport in electrodialysis cells are presented and analysed. Some examples are given when specially designed surface heterogeneity on the membrane surface results in enhancement of ion transport in intensive current electrodialysis.


Introduction: Heterogeneity and Multiscale Nature of Ion Transport in Membrane Systems
The main aim of this paper is to give a review of the well-established models as well as new approaches for describing the transport phenomena in charged membranes. These materials include ion-exchange (IEM), nanofiltration (NF) and reverse osmosis (RO) membranes, which are widely used in separation and energy production processes.
The peculiarity of mass transfer phenomenon is membrane systems is the multiscale character of ion and water transport, which is due to the heterogeneity at multiple scales [1]. The charged membranes are self-assembled nanostructured materials constructed from macromolecules [2]. It is remarkable that the structure and, hence, the properties of both biological and artificial membranes show considerable similarities [3,4]. The main components of the structure are amphiphilic compounds, which have both hydrophilic and hydrophobic constituents. The mutual interaction of the hydrophobic segments and the interaction of polar hydrophilic groups with aqueous medium lead to a self-organized stable structure. In particular, nanoscale hydrophilic channels which are selectively permeable to water and ions are formed within a hydrophobic matrix [5]. The size of transport channels is of a few nanometres, so that the electrical double layers (EDL) formed at the opposing pore walls overlap. Overlapping EDLs leads to an essentially higher concentration of the ions whose sign of the charge is opposite to the that of the fixed ions (the counterions) in comparison with the concentration of the ions having the same sign of the charge as the fixed ions (the co-ions).
The second scale is due to the diffusion boundary layers (DBLs) at the membrane surface formed under the action of the electric current flowing across the membrane. This scale is of the order of 100 µm. When an overlimiting current density flows, an intermediate scale determined by the thickness of the extended space charge (of the order of 1 µm) appears.
Finally, there is a scale of the order of 10 cm determined by the length of a compartment in an electrodialysis (ED) stack used in industrial process.
The processes occurring at one scale are coupled with the processes occurring at other scales. For example, selective ion transport in membrane pores causes concentration polarization, that is, formation of concentration gradients in the DBLs of bathing solution.
In this paper, we consider the transport at all the mentioned above scales and examine the possibilities of multiscale modelling. The latter is rather difficult, as such modelling needs large computer resources. We give the mathematical description of ion transport in the membranes taking into account the structure-properties relationships. Concentration polarization occurring at the scale of DBL and related surface phenomena such as water splitting along with current induced convection are described as well. Their mechanisms are analysed in close relation with the membrane surface properties. Here we focus on the role of the membrane surface heterogeneity, which has been subject to intensive attention in the literature of recent years. Effects associated with desalination of the solution in the longitudinal direction along the length of an ED compartment are also studied.

Structure of Ion-Exchange Membranes
Synthetic IEMs are built using flexible polymer chains with aliphatic, aromatic or perfluorinated residues, which contain functional groups such as -SO 3 H, -COOH, -NH 3 OH and others. Protons or OH − -ions of these groups can be replaced by cations or anions, respectively, from bathing solutions.
The diversity of practical applications of membranes conditions a variety of specific requirements expected from them and, eventually, promotes the development of a wide range of membrane materials. A considerable number of membranes and membrane materials are now available [6][7][8].
Depending on the material structure and the method of preparation, IEM are usually divided into homogeneous and heterogeneous membranes. Homogeneous membranes (Figure 1a) are generally obtained by copolymerization of monomers [9] that produces homogeneous material at least at the scale of 100 nm or more, which is determined by the size of polyvinyl chloride (PVC) particles used as a non-conductive binder in the Neosepta ® membranes [10]. Heterogeneous membranes include macroparticles (with the size from 1 to 50 µm) of different polymer materials. For example, the MA-40 anion-exchange membranes (AEMs) (Shchekinoazot, Tula Region, Russia) contain particles of polystyrene crosslinked with divinylbenzene (DVB) anion-exchange resin and polyethylene used as the binder (Figure 1b). Both homogeneous and heterogeneous membranes contain reinforcing cloth (Figure 1a,b). This cloth not only plays an important role as a means to improve the mechanical properties but also causes a certain undulation of the membrane surface, which seems important for initiating electroconvection in overlimiting current regime [11]. Recently, Svoboda et al. [13] developed a specific scanning cell for the analysis of heterogeneous IEM structure by micro-computed tomography, which provides volume reconstruction of the studied material. The unique working principle of the cell allows one to scan membranes not only in their dry state but also in fully swollen state. The micro-computed tomography showed that there are profound changes in the structure of the IEM associated with their transition from the dry state into the swollen one. This study applied to the heterogeneous cation-exchange membrane (CEM) coming from Mega a.s., Czech Republic showed that swelling of the membrane is caused by the swelling of the ion-exchange resin. Its volume increased nearly three times after exposure to deionized water, while the overall membrane volume increased by 50%. If the membrane contained approximately 28%vol. of the resin in the dry state, this content increased to 52%vol. after swelling.
Micro-computed tomography not only allows one to analyse volumetric composition of the membrane bulk but also to quantify the composition of the surface. By using two different methods, Vobecka et al. [14] evaluated the fraction of the conductive area (e.g., the ion-exchange resin) on the surface of an AEM and CEM in the swollen state. The obtained results showed that the volume fraction of the conductive domains does not correspond to the surface fraction of these domains and thus the volume composition is not reflected on their surface. For example, the volume fraction of the conductive domains for the studied AEM was 0.45 and its surface fraction only 0.26. Similar results were obtained earlier by Volodina et al. [15] and other authors [16]. The reason is that during the hot pressing of heterogeneous membranes, polyethylene, which is more fluid than polystyrene, moves from the membrane bulk on the surface.
As it is mentioned above, the terms "homogeneous" and "heterogeneous" IEM are commonly used in the field of IEMs to distinguish between the membranes produced using the paste method and the membranes made via the hot pressing of powdered resins and inert binder. However, "homogeneous" membranes are not strictly speaking homogeneous. For instance, Nafion® is a perfluorinated sulfonated ion-exchange dense polymer generally considered as homogeneous [28]. Nevertheless, when swelled in water, it contains hydrophilic pores/channels confined within a Recently, Svoboda et al. [13] developed a specific scanning cell for the analysis of heterogeneous IEM structure by micro-computed tomography, which provides volume reconstruction of the studied material. The unique working principle of the cell allows one to scan membranes not only in their dry state but also in fully swollen state. The micro-computed tomography showed that there are profound changes in the structure of the IEM associated with their transition from the dry state into the swollen one. This study applied to the heterogeneous cation-exchange membrane (CEM) coming from Mega a.s.,Česká Lípa District, Liberec Region, Czech Republic showed that swelling of the membrane is caused by the swelling of the ion-exchange resin. Its volume increased nearly three times after exposure to deionized water, while the overall membrane volume increased by 50%. If the membrane contained approximately 28%vol. of the resin in the dry state, this content increased to 52%vol. after swelling.
Micro-computed tomography not only allows one to analyse volumetric composition of the membrane bulk but also to quantify the composition of the surface. By using two different methods, Vobecka et al. [14] evaluated the fraction of the conductive area (e.g., the ion-exchange resin) on the surface of an AEM and CEM in the swollen state. The obtained results showed that the volume fraction of the conductive domains does not correspond to the surface fraction of these domains and thus the volume composition is not reflected on their surface. For example, the volume fraction of the conductive domains for the studied AEM was 0.45 and its surface fraction only 0.26. Similar results were obtained earlier by Volodina et al. [15] and other authors [16]. The reason is that during the hot pressing of heterogeneous membranes, polyethylene, which is more fluid than polystyrene, moves from the membrane bulk on the surface.
As it is mentioned above, the terms "homogeneous" and "heterogeneous" IEM are commonly used in the field of IEMs to distinguish between the membranes produced using the paste method and the membranes made via the hot pressing of powdered resins and inert binder. However, "homogeneous" membranes are not strictly speaking homogeneous. For instance, Nafion ® is a perfluorinated sulfonated ion-exchange dense polymer generally considered as homogeneous [28].
Nevertheless, when swelled in water, it contains hydrophilic pores/channels confined within a hydrophobic matrix [5]. These pores/channels provide for the transport of ion and water through a Nafion membrane. The presence of internal solution within a hydrophobic matrix allows one to speak about phase separation [5]. From this point of view, this membrane is multiphase and heterogeneous, however, this heterogeneity appears on a scale of the order of 10 nm.
The main structural feature of the functional part of IEMs (which is an ion-exchange resin in the case of heterogeneous membranes) is the hydrophobic backbone of hydrocarbon or perfluorinated chains, which antagonizes the hydrophilic character of the charged functional groups. Due to the flexibility of the polymer chains, this combination provides self-organization during the membrane formation [29]. Polymer chains form the matrix of the membrane, whereas the functional groups are aggregated in clusters. The size of the clusters is of the order of a few nanometres. Consequently, the hydrophobic components contribute to the morphological stability of the membrane. The hydrophilic domains create a system of conducting channels. When placed in an aqueous solution, hydrophilic functional groups are hydrated and the membrane swells. The size of clusters increases, conducting channels appear between them and they form a percolation system when a certain water content is reached [30]. The water content and degree of swelling depend on the external solution concentration: with increasing solution concentration water content of an IEM decreases [30,31]. When a membrane contains weakly basic functional groups, such as secondary and tertiary amino groups, there is a dependence of water content on the external solution pH [31][32][33]. The model of clusters and channels first proposed by Gierke [34] describes these main features and explains correctly the transport properties of IEMs [5,30,35]. A rather realistic presentation of two-dimensional structure of Nafion ® , which generally complies with the Gierke model, is made by Kreuer [36]. The clusters, channels, some defects of the structure, spaces between ion-exchange resin particles, the binder and the fabrics of the cloth form a system of pores in an IEMs, whose size varies from a few nm to 1-2 µm [13,[37][38][39]. The investigations of the pore size distribution made using the standard contact porosimetry (Divisek et al. [40]; Berezina et al. [37] and Kononenko et al. [38,39]) and differential scanning calorimetry (DSC)-based thermoporosimetry (Bryk et al. [41]; Berezina et al. [42]; Kononenko et al. [39]) methods have shown that the homogeneous IEMs (namely, perfluorinated MF-4SK membranes (Plastpolymer, Russia) [39,41,42]), as well as Nafion-112, Nafion-115, Nafion-117 membranes [40]) do not contain macropores (the pores, whose effective radius size is larger than 100 nm), while heterogeneous membranes have such pores. It was found [37,39,42] that the pore size distribution of Russian heterogeneous MK-40 and MA-40 membranes has two maxima, the first one is at about 10 nm (which is the range of micro-and mesopores) and the second one relates to macropores with size of about 1000 nm. According to [39], the pores of the first type are situated within the ion-exchange resin particles and the pores of the second type are formed by the spaces between different membrane constituents (resin, binder, cloth).
The solution filling the clusters where the fixed functional groups are concentrated is charged. It contains mainly counterions coming from the dissociation of functional groups. The counterions being in thermal motion are attracted by the fixed charged groups, forming in this way an EDL. A pore also contains a small amount of co-ions carrying a charge of the same sign as that of the fixed ions. The co-ions are repelled from the pore by electrostatic forces; this effect is called the Donnan (co-ion) exclusion [31]. The distribution of ions within the EDL formed at inner or outer interfaces is described by the Gouy-Chapman theory, which takes into consideration both electrostatic and thermal interaction. Figure 2 schematically shows the ion distribution in the proximity of a pore wall represented as a charged plane. The concentration of counterions in the pore solution increases and that of co-ions decreases with decreasing the ratio of the pore radius to the Debye length. The latter characterizes the EDL thickness. Therefore, the less the pore radius, the higher the membrane selectivity (and the higher the membrane resistance). When the pore radius is higher than the Debye length, there is an electrically neutral solution in the pore centre. When the water molecules are close to the functional groups, they are structured and lose to a great extent their mobility [3,5]. As a result, the relative dielectric permittivity, ε, decreases when approaching a functional group. The mobility of counterions concentrated near the functional groups is very low mainly due to strong electrostatic interaction [3,5]. The side-chains also contribute to their low mobility by encumbering the space [43]. Due to the finite size of hydrated counterions, their concentration decreases when the distance from a fixed group becomes lower than about 0.5 nm. When getting close to the pore centre, the counterion concentration decreases as well, according to the Gouy-Chapman law of ion distribution in EDL ( Figure 2).
The transport properties of the central zone of a hydrophilic domain are different from those of a zone, which is close to the internal interface. The water molecules located in the centre of sufficiently large pores behave as in free solution. In particular, ε approaches the value of 81 there, as in free solution [44]; the proton transport mechanism in the pore centre is identical to that in free solution. The Grotthuss shuttling dominates over the vehicular transport (solvation diffusion): proton jumps from one H 3 O+ (or H 5 O 2 + or even greater) complex to a neighbouring one rather than moves within its aqueous environment-hydration shell [44]. The contribution of the vehicular transport increases with increasing the degree of water structuring, which takes place in the vicinity of the fixed charged groups. The volume fraction of bulk-like water depends on the membrane water content. Electronic structure calculations [45] show that 2-3 water molecules (n) per sulfonic acid group are needed in perfluorinated membranes (such as Nafion ® ) for proton dissociation. The dissociated proton is separated from the sulfonate anion, when 6 water molecules are added in a membrane. According to Kreuer et al. [5], only when n > 14, one can speak of a two-phase system where the bulk-like water is clearly distinguishable in the pore.
Schematically, the structure of a functional part of an IEM is shown in Figure 3. It is a fragment including hydrophobic domains of the matrix, micropores (where the EDLs are overlapped) and mesopores (where the EDLs are not overlapped but occupy an essential portion of a pore). This fragment is typical for a homogeneous IEM. Heterogeneous IEMs contain also macropores (spaces between ion-exchange resin particles and reinforcing binder and tissue). To pass from one fixed ion to a neighbouring fixed ion, a counterion has to overcome a potential barrier. The height of this barrier depends on the energy of electrostatic interactions between the hydrated counterion and the fixed ion, on the chemical interaction between this pair of ions and on the energy required to move the polymer chains to form a sufficiently large ion transport channel in the matrix. The average distance between two neighbouring fixed ions, b, depends on the exchange capacity and membrane morphology. For common membrane materials, b is in the range 0.5-1.0 nm [3]. This distance for the fixed ions in a hydrated cluster of perfluorosulfonic membranes is about 0.8 nm [3,5].
Kamcev et al. [47][48][49] propose to use the Bjerrum length as the distance from the fixed ions, which determines the boundary dividing the electrolyte solution in membrane pores into two parts whose properties are different. The Bjerrum length, λB, is the distance at which the electrostatic potential energy between two charges equals the thermal energy scale [50]. According to Bjerrum, if a counter-ion is closer to the fixed charge than λB, it is regarded as "associated": its electrostatic interaction with the fixed ion is so great that it cannot diffuse away from the fixed charge [48]. If the distance between two ions is less than λB, the energy of electrostatic interaction exceeds the kinetic energy of thermal motion, which reduces the mobility of ions. Manning [51,52] proposed a model for polyelectrolyte known as the Manning condensation theory [53]. Kamcev et al. [47][48][49] have applied recently this theory to ion-exchangers. According to [49], if the distance between two fixed charges (b) is less than λB, counterions condense (or localize) on the polymer chain, which is a result of strong electrostatic interaction between the fixed ions and mobile counterions. When b > λB, counterion condensation does not occur. The value of λB in water is 0.71 nm, in IEMs with the IEC close to 2 meq. (g dry membrane) −1 ] its value is 1.2-1.3 nm [49] (due to lower dielectric permittivity in IEMs). The value of b in the IEMs studied by Kamcev et al., was evaluated between 0.43 and 0.57 nm [49].

Irreversible Thermodynamics Approach
There are several approaches for describing ion and water transport in membranes. Some books [54][55][56][57] and reviews [58][59][60] deal with this topic. Phenomenological approach based on the irreversible thermodynamics provides a rather simple mathematical description, which, however, is not directly related to the membrane structure. This approach will be considered in this section. The . Schematic representation of the IEM structure with its main elements: fixed ions (shown as circles with "−"), EDL formed at the internal interfaces and electroneutral solution in the centre zone of inter-gel spaces. The gel phase includes the polymer matrix bearing the fixed ions and the EDL. Adapted with permission from [46], Elsevier, 2008.
To pass from one fixed ion to a neighbouring fixed ion, a counterion has to overcome a potential barrier. The height of this barrier depends on the energy of electrostatic interactions between the hydrated counterion and the fixed ion, on the chemical interaction between this pair of ions and on the energy required to move the polymer chains to form a sufficiently large ion transport channel in the matrix. The average distance between two neighbouring fixed ions, b, depends on the exchange capacity and membrane morphology. For common membrane materials, b is in the range 0.5-1.0 nm [3]. This distance for the fixed ions in a hydrated cluster of perfluorosulfonic membranes is about 0.8 nm [3,5].
Kamcev et al. [47][48][49] propose to use the Bjerrum length as the distance from the fixed ions, which determines the boundary dividing the electrolyte solution in membrane pores into two parts whose properties are different. The Bjerrum length, λ B , is the distance at which the electrostatic potential energy between two charges equals the thermal energy scale [50]. According to Bjerrum, if a counter-ion is closer to the fixed charge than λ B , it is regarded as "associated": its electrostatic interaction with the fixed ion is so great that it cannot diffuse away from the fixed charge [48]. If the distance between two ions is less than λ B , the energy of electrostatic interaction exceeds the kinetic energy of thermal motion, which reduces the mobility of ions. Manning [51,52] proposed a model for polyelectrolyte known as the Manning condensation theory [53]. Kamcev et al. [47][48][49] have applied recently this theory to ion-exchangers. According to [49], if the distance between two fixed charges (b) is less than λ B , counterions condense (or localize) on the polymer chain, which is a result of strong electrostatic interaction between the fixed ions and mobile counterions. When b > λ B , counterion condensation does not occur. The value of λ B in water is 0.71 nm, in IEMs with the IEC close to 2 meq. (g dry membrane) −1 ] its value is 1.2-1.3 nm [49] (due to lower dielectric permittivity in IEMs). The value of b in the IEMs studied by Kamcev et al., was evaluated between 0.43 and 0.57 nm [49].

Irreversible Thermodynamics Approach
There are several approaches for describing ion and water transport in membranes. Some books [54][55][56][57] and reviews [58][59][60] deal with this topic. Phenomenological approach based on the irreversible thermodynamics provides a rather simple mathematical description, which, however, is not directly related to the membrane structure. This approach will be considered in this section. The structure-properties relationships may be described using other approaches, which take into account parameters of membrane structure. These approaches, considered below, involve the "solution-diffusion" and "pore-flow" models.

The Onsager Phenomenological Equations
The strength but also the weakness of the description of ion and water transport across the membranes in the framework of irreversible thermodynamics is that information about the membrane structure and transport mechanism is not used. This approach gives general relationships between fluxes and driving forces, which are valid for any type of membrane. The specificity of a given membrane is taken into account through a set of phenomenological coefficients, which cannot be found within this approach. Usually, these coefficients are determined in specially designed experiments, which is sufficient for engineering purposes. However, a more theoretical way is possible, when the transport equations of irreversible thermodynamics are completed with other equations/boundary conditions, taking into account the details of membrane structure.
At the equilibrium state, the electrochemical potential, µ i , of any mobile species i (an ion or a molecule) is the same in every point of the membrane material. If nonzero gradients dµ i /dx are applied, fluxes of different species appear. If the system is close to equilibrium, linear relations between the thermodynamic forces, F j = −dµ i /dx and resulting fluxes, J i , can be written [56,61]: where L ij are the phenomenological conductance coefficients. In addition to the material properties, L ij , can depend on the species concentrations, temperature and pressure but not on fluxes and forces. It is assumed that no net transfer occurs in the direction parallel to the membrane surface, and, therefore, the value of µ i is the same at every point of the plane normal to this direction, regardless of the phase through which the plane passes ( Figure 3). dx in Equation (1) can be interpreted as a distance between two planes perpendicular to the transport axis x, the first one corresponding to the value of µ i on the plane at the distance x, µ i (x), the second one to µ i (x + dx) [62,63].
Note, that Equation (1) shows that J i depends not only on the force applied to species i, F i but also on all forces applied to other species.
It is suitable to present the L ij coefficients as functions of the "virtual" electroneutral solution concentration. This solution is assumed to be in local equilibrium with a small volume of the membrane. The notion of virtual solution was first introduced by Kedem and Katchalsky, who named it "corresponding solution" [64]. This solution may be physically present in the central zone of meso-and macropores (where the pore radius is higher than the Debye length) ( Figure 3); if not (the EDLs at the opposite walls of the pore are overlapped), it may be considered as a hypothetical solution. Since usually local equilibrium at external membrane interfaces takes place, the virtual solution is identical to the external bathing solution at each membrane side.
The electrochemical potential µ i in Equation (1) can be present as a function of activity a i of species i, electrical potential φ and pressure p in the virtual solution: where V i is the partial molar volume and z i is the charge of species i; F, R and T are the Faraday constant, the gas constant and the temperature, respectively. The choice of thermodynamic forces (F i ) and conjugated fluxes (J i ) in Equations (1) may be different but not arbitrary: when the choice of J i and F i is correct, the sum of products J i F i gives the dissipation function [56,61]. The transport coefficients' values depend on the choice of forces and fluxes. Narebska et al. [65] carried out delicate and laborious experiments and found coefficients L ij of the Onsager equations, Equation (1), for a Nafion ® 120/NaCl membrane system, as functions of NaCl solution concentration. For the case of use other forms of transport equations, with other set of the forces and fluxes, the phenomenological coefficients were also reported in Refs. [66,67].
Along with the Onsager equations, Equation (1), there are other equation systems (Stefan-Maxwell, Spiegler, Kedem-Katchalsky and others [68][69][70][71]); their review is given in References [54,56]. The difference in these systems is in the choice of fluxes and forces, and, as consequence, in the transport coefficients. It is possible, by a simple transformation of variables, to pass from one system to another [67], hence, these systems are mathematically equivalent.
The famous Onsager reciprocal relations allows reduction of the number of independent phenomenological coefficients. If three different species (a counterion, a co-ion and a solvent) are present in the membrane or solution, the number of independent coefficients is 6.

The Kedem-Katchalsky Equations
The system of Kedem-Katchalsky equations [68] is of great interest for the practical description of ion and water transfer through the membranes. These equations in differential form are written as follows [67]: where three thermodynamic fluxes are presented by J v , J i and i (which are the volume and solute flux densities and the current density, respectively); c s and c i are the molar concentration of salt and ion i, respectively, in the virtual solution of the membrane; p and π are the hydrostatic and osmotic pressures and φ is the electric potential in the virtual solution, respectively; ν = ν + + ν − is the stoichiometric number; d is the membrane thickness; subscripts "+" and "−" relate to cation and anion, respectively. Equations (4)- (6) can be applied to a CEM and an AEM, in which cations and anions act as counterions, respectively. Three thermodynamic forces in Equations (4)-(6) are as follows: mechanical, dp/dx, electric, dφ/dx and chemical. The latter is expressed through the gradient of osmotic pressure or that of concentration, which are linked between them [68]: where a ± is the mean ionic activity of the electrolyte. The major interest of the Kedem-Katchalsky equations is in the fact that their transport coefficients are those usually used in the characterization of membrane transport properties. That is why these six transport coefficients are called practical: hydraulic permeability (L p ); electro-osmotic permeability (β); diffusion permeability (P); transport number of the counter-ion (t + ); electrical conductivity (κ); and Staverman reflection coefficient (σ). The meaning of the latter can be understood from considering two limiting cases: σ = 1, when the membrane completely reflects the solute transported with convective flow across the membrane and σ = 0, when the retention is zero. The physical meaning of coefficients P, κ and L p is the same as that in the linear equations of Fick, Ohm and Darcy, respectively. These linear equations are easily derived from Equations (4)- (6), if only one of three driving force is kept non-zero. t + is the fraction of electric charge carried by cation under the condition that the gradients of concentration and pressure are zero; β, the electroosmotic permeability, should be determined as the proportionality coefficient between the volume flux and the current, when the gradients of concentration and pressure are zero.
A form of Kedem-Katchalsky equations different from Equations (4)- (6) can be also used [68]: the pressure gradient in Equations (5) and (6) is replaced by the volume flux. This slightly changes the meaning of coefficients P, κ and t + [68].
Equations (4)-(6) form the basis of membrane characterization [37,65,[71][72][73][74][75][76], since each of the coefficients refers to one of the membrane transport properties. A transport coefficient should be measured when only one driving force is applied. For example, diffusion permeability P is to be obtained when ∇p = i = 0. More details about measuring protocols can be found in Refs. [37,67,76].
As it was mentioned previously, one can pass from the Onsager to the Kedem-Katchalsky equations by a mathematical transformation, since the coefficients in these two systems of equations are linked [67,77]. The expressions linking κ, t i and P with the Onsager coefficients are as follows [67]: where refer to the individual ions "+" and "−" respectively; g = 1 + d ln y ± /d ln c is the activity factor written for the virtual solution [67,78], y ± is the mean (molar) activity coefficient of electrolyte in the virtual solution; t +app is the apparent transport number of counter-ion (cation in a CEM). According to the Scatchard equation, t +app = t + − mM w t w , where m is the virtual solution molality (in eq/kg H 2 O), M w is the water molar mass and t w is the water transport number in the membrane. The frequently applied approximation (presented in Equation (10) by its right-hand side after "≈") is found when assuming L +− − m s M w L −w = 0, t +app = t + and g = 1. Note that the central part of this equation (after sign "=") is presented for 1:1 electrolyte, whereas the right-hand side is for any single electrolyte; here c = |z i |c i is the electrolyte concentration in the virtual solution expressed in eq L −1 . Similar expressions for L p , β and σ are presented in Ref. [67].
The laborious and delicate experimental determination of Kedem-Katchalsky and Onsager coefficients can be facilitated when using an incomplete set of phenomenological coefficients. For instance, the Staverman coefficient σ, which is very close to 1 in IEMs, is rarely used in characterization. Besides, additional approximate relations between the coefficients can be applied [67,71,77]. Such a relationship, presented by Equation (12), can be obtained if we assume that the difference between the L −+ coefficient and the m s M w L −w term in Equation (10) is small compared to the L −− diagonal coefficient [67,77]: Equation (12) is a generalization of the Nernst-Einstein equation for the case where the osmotic and electro-osmotic water transport are implicitly taken into account. This and similar equations have been applied and compared with experiment in different conditions [67,72]. A slight deviation from Equation (12) was found [67] in the case of perfluorinated Nafion ® 120 and MF-4SK membranes. It was interpreted as a consequence of a relative hydrophobicity of these membranes according to Kedem and Perry [72]. Equation (12) and similar relations are successfully implemented in the characterization of IEMs [67,71,73,75,79].
In practice, Equation (12) is more convenient to be applied in the form: The values of P and κ are easier to measure than the transport numbers. The apparent transport number, t +app , is found from the membrane potential measurements [55,67]: where E is the potential difference (PD), measured between two reversible electrodes situated from both sides of the membrane; a I ± and a I I ± are the electrolyte activities in the left-and right-hand side bathing solutions; E max is the theoretical value of E for an ideal membrane not permeable for co-ions. The activity factor, g, can be found from the literature [75,76,78]. For LiCl, NaCl and NaCl solutions, g equals 1 at zero concentration, it passes through a minimum close to 0.9 (in the range 0.1-0.2 M) and regains the value 1 at 1 M. Thus, if we use the approximation t +app = g = 1, when applying Equation (13) in the case where t +app is not too low (>0.9), we find a slightly underestimated (by 10-20%) value of t − . The true counterion transport number, t + = 1 − t − , will be find with a reasonable accuracy: it would be overestimated by 1-3%. The accuracy can be improved when using approximation t +app ≈ t + = 1 − t − In this case Equation (13) is reduced to: Note that Equations (8)- (15) link local values relating to the (virtual) solution with a concentration c. It is relatively easy to measure the conductivity of a membrane equilibrated with a solution of concentration c. The membrane integral (or global) permeability, P int , is measured in conditions where one of the bathing solutions is of concentration c s and the other one is water. To find the differential (local) permeability, P, the concentration dependence of P int should be measured. Then [37,76]:

The Nernst-Planck Equation
The Nernst-Planck equation is a partial case of the transport equations presented above. It is most simple to derive this equation from the Onsager Equation (1) [80]. If one neglects the cross phenomenological coefficients in Equation (1), sets p = const and assumes that the activity coefficients are constant, one immediately finds [78,80]: where is the diffusion coefficient of ion i. The electrochemical potential gradient can be expressed through the ion activity (a i ) or concentration (as in Equation (17)) and the electrical potential (φ) in any constituent phase of the membrane material. The use of the virtual solution is often convenient, since in this case a i (c i ) and φ have no jumps at the external membrane boundaries (if only the condition of local equilibrium is satisfied at these boundaries, which is verified for the under-limiting current densities [81][82][83]). More often, Equation (17) is applied to a solution and/or a membrane considered as a homogeneous medium [31,57,84]. In this case, the concentration and electric potential refer to this medium (which can be imagine as a swollen sponge with a charged matrix [31]) and not to a particular phase. To indicate the fact that a quantity refers to the membrane material (or its part) considered as a homogeneous medium, we will use an over-bar: c i , ϕ and D i . The models, which use the approach where the permeants are assumed to be dissolved in the membrane material and transported there under the action of concentration gradient and/or external electric field, are called the "solution-diffusion" models [85,86]. This term is used in contrast to the term "pore-flow" models, in which the transport of species is considered as occurring within a single membrane pore. Usually, this terminology is used in the context of the pressure-driven processes [86,87], while it is also applied to the electromembrane processes [88,89], where the transport occurs under the action of the electrochemical potential gradient.
The Nernst-Planck equation is used as the governing equation in both types of the models mentioned above [90]. In the "solution-diffusion" models, it is applied in the membrane material considered as quasi-homogeneous medium (effective medium approach [91]), where the assumption of local electroneutrality (LEN) is commonly used. In the "pore-flow" models, this equation is applied within a single liquid phase, which is the pore solution. Usually, it is coupled there with the Poisson equation, thus taking into account the electric space charge in the EDL at the charged pore wall.
It is also often convenient to use the transport equation in the membrane in the reduced Kedem-Katchalsky form: Equation (18) is obtained from Equation (5) of the Kedem-Katchalsky set of equations when neglecting the term containing dp/dx. Equation (18) is written in two ways. The first one is applied when the membrane is considered as a quasi-homogeneous gel, in this case the quantities are denoted with an over-bar; the concentrations refer to the number of moles per unit volume of the gel. When the concentrations refer to the virtual electroneutral solution (c i ) (the second way of writing), the differential diffusion permeability coefficient, P, is used instead of the diffusion coefficient, D. In both cases, the transport number has the same meaning: it is the fraction of electric charge transported by ion i in the membrane in conditions where the concentration gradient is zero. It is possible also to write this equation for a solution. In this case, the electrolyte diffusion coefficient in solution is used instead of P.
Equation (18) is equivalent to the Nernst-Planck equation, if the LEN assumption is accepted [80]. While Equation (18) does not describe the pressure-driven transport, the impact of current-induced (electroosmosis) or diffusion-induced (osmosis) water transport can be implicitly taken into account [88]. Water transport gives rise to the ion convection transport. Therefore, electromigration and/or diffusion of ions occur in a mobile medium, which affects the effective value of the diffusion coefficient. This effect in the thermodynamics of irreversible processes is called the frame of reference effect [92].
It is easy to see from Equation (18), which is valid in the framework of the "solution-diffusion" approach, that the electrolyte permeability coefficient, P, is related to the electrolyte diffusion coefficient, D and the electrolyte partition coefficient, K s , as follows [85]: The definition of K s (i.e., K s = c − /c − ) is clear from Equation (19). It is taken into account that the electrolyte concentration in the membrane is equal to the co-ion concentration c − co-ions. In strongly charged IEMs, K s is essentially less than 1 due to Donnan exclusion of co-ions. With increasing external concentration, the Donnan exclusion gets less strong and K s increases [46,48,93,94].
The electrolyte diffusion coefficient and the ion transport number in the membrane are linked with the individual ion diffusion coefficients via the following relations: Equations (20) can be simplified, if the LEN assumption is used: where Q is the concentration of functional fixed groups (in meq. cm −3 ), that is, the ion-exchange capacity (IEC).
There are several generalizations of the Nernst-Planck equation including described above the Kedem-Katchalsky and Onsager equations. The following equation takes into account the concentration variation of the activity coefficients (through the introduction of the activity factor, g, [defined just after Equation (11)] and the convective flow (by adding the term with the centre of mass velocity, V) [80]: Equation (22) is called the extended Nernst-Planck equation.
Convective transfer is the main transport mechanism in pressure-driven processes [95,96], in dialysis across biological or biomimetic membranes [97]. In ED, water transport has to be taken into account in the applications aimed at the production of concentrated solutions (brines), since water moving into the concentration compartment dilutes the brine. Water transport across IEMs affects also the characteristics of energy generation from salinity differences via reverse ED [98].
When applying the approach of the virtual solution, g can be found from the literature data [78] on the activity coefficient as a function of concentration [46,75,76]. Another approach is developed by Kamcev et al. [48,88], who calculated g for a membrane considered homogeneous by applying Manning's counter-ion condensation theory [51]. The theory provides expressions for individual ion activity coefficients as functions of ion concentration in the membrane.

Modelling the Structure-Property Relationships
As it was mentioned above, the models linking the structure and transport properties can be divided in two main groups: the "solution-diffusion" and "pore-flow" models. The "solution-diffusion" models involve the models considering the membranes as homogeneous medium and those, which assume that the membrane may be represented as a multiphase system. According to the multiphase models, the membrane properties are functions of the properties of the constituent phases as well as their volume fractions and their relative disposition within the membrane. The "pore-flow" models operate with such parameters as the pore radius and sometimes length. During the last years, in addition to these kinds of models, a significant progress was achieved by applying molecular dynamics (MD) to describe macroscopic behaviour of ions and water in a membrane pore [45,[99][100][101][102].

Teorell-Meyer-Sievers (TMS) Model
Historically, the first model successfully applied to explain the main properties of charged membranes, was the model developed in thirties years of the last century independently by Teorell [103] and Meyer and Sievers [104], called the TMS model. It provides the basis of mathematical description of the ion transport in membranes.
The main component of the TMS model is the Nernst-Planck equation, Equation (17). The membrane is considered as a single homogeneous phase (charged gel), representing an aqueous solution of matrix polymer chains together with mobile and fixed ions. Along with Equation (17), the model involves the LEN assumption in the membrane, Equation (21), the equation expressing the electric current as a function of the individual ionic flux densities: and the Donnan equilibrium at the left-hand (I) and right-hand (II) interfaces used as the boundary conditions: where K D is the Donnan equilibrium coefficient; c k i and c k i are the concentrations of ion i at interface k, from the side of the membrane and solution, respectively [105].
Equations (17), (21)-(24) form a boundary-value problem for the electrodiffusion of ions across an IEM placed between two solutions of known component concentrations.
The Donnan relation, Equation (24), is obtained from the condition of continuity of electrochemical potential written for each ion species i at the membrane/solution interface k: The ion activity coefficients implicitly entering into Equation (25) are generally functions of the pressure in each phase. It is easy to derive the following two relations from Equation (25): one for the activities (a k + ) and the other for the electric PD between two phases: When introducing ionic concentrations (mmol cm −3 ) in Equation (26) instead of activities (a i = c i y i and a i = c i y i ), we arrive at Equation (24) with the Donnan equilibrium coefficient, K D , expressed through the ratio of mean ionic activity coefficients: where y ± = y ν + + y For the mean ionic activity coefficient in the membrane y ± , a similar expression can be written.
The activity coefficients in water solutions can be determined from the data in the literature [78,106]. Recently, Kamcev, Paul and Freeman [47] have calculated the ionic activity coefficients using Manning's counterion condensation theory. The theory was originally developed for polyelectrolyte solutions [51] and contains no adjustable parameters. It is one of the most recognized theories, which describe colligative properties of polyelectrolyte solutions [51]. In this theory, a polyelectrolyte chain is modelled as a linear, infinitesimally thin and infinitely long rod with fixed charges evenly spaced along the backbone.
The mathematical description involves a single parameter, ξ, related to the linear charge density of the polyelectrolyte chain [51]. ξ is a ratio of two length scales, the distance between neighbouring fixed charges and the Bjerrum length, λ B . The counterions located closer than λ B to a functional charged group do not have sufficient thermal energy to freely diffuse away from this group. They are considered as "associated" as opposed to "free" ions, which are separated from this group by a distance greater than λ B [51,107,108].
The critical linear charge density characterizing the onset of counterion condensation in polyelectrolyte solutions was found experimentally by a number of researchers confirming Manning's counterion condensation theory. A comprehensive review of these publications is given by Manning [109].
While the model is rather simple and uses no adjustable parameters, there is a good agreement between theoretical and experimental values of activity coefficients for many polyelectrolyte systems [51,109,110]. The projected length of one monomer unit was evaluated as 0.25 nm [111]; the average distance between fixed charges for a Nepton ® CR61 membrane, as 0.73 nm [47].
The TMS models gives an adequate qualitative description of such fundamental membrane properties as electric conductivity, transport number and membrane potential [112]. When two kinds of counterions, 1 and 2, are present in the membrane, the Nernst-Planck equation together with Equations (24) and (27) allow description of counterion competitive transfer. For more information on the TMS model, see Refs. [31,55,57,93,113].
Similar assumptions are accepted in some recent models of ion sorption and transfer in charged membranes regarded as homogeneous medium [114][115][116][117][118][119]. Filippov et al. [114][115][116][117] use the so-called fine-pore model, which allows description of some special cases, for example, when an organic IEM contains inorganic nanoparticles. Fila and Bouzek [118,119] use the same assumption as the TMS model but apply the extended Nernst-Planck equation, which contains the convective term. That is important when considering the application of IEMs in chlor-alkali electrolysis.
There are models considering the membrane as a single phase (charged gel), in which the concentration of fixed charges and/or diffusion coefficients continuously change along the normal coordinate [120,121]. In this case, the TMS theory can be applied locally. Generally, these models show that a heterogeneity in the fixed charge distribution leads to an increase of permselectivity in comparison with a membrane wearing homogeneously distributed fixed charges of the same average concentration. The explanation is that the permselectivity is controlled by the layer having the highest local concentration of fixed charge, while the layers with low fixed charge concentrations have a minor impact on the global membrane behaviour. For example, the use of conical and double-conical nanopores in NF membranes allows improving the performance of these membranes (the salt rejection) [122]. Two different phenomena are at the basis of the filtration rectification properties of conical nanopores: the co-ion exclusion at the pore mouth and the sign and the magnitude of the (pressure-induced) electric field arising through the nanopores. Hourglass-shape nanopores show improved separation performance compared with cylindrical and conical nanopores having identical average diameter, length and surface charge density. Notably, they exhibit a higher salt rejection with a smaller driving force than the other pore geometries. That makes them attractive candidates for the design of advanced NF membranes.

Multiphase Models
As it was discussed previously (Section 2), charged membranes are porous materials having generally micro-, meso-and macropores. For the definition of these types of pores, several characteristics can be applied (such as the length of action of adsorption forces) [123]. However, here we will use only one parameter to distinguish the pore type: the ratio of the pore radius, r, to the Debye length, λ, characterizing the EDL thickness at the inner pore wall. In the case of r < λ, we have a micropore (r < 1 nm); if r > λ but r and λ being of the same order of magnitude, we deal with a mesopore and if r >> λ (r > 50 nm), it is a macropore. Despite the complicated structure of charged membranes, it is important that mesopores and macropores contain electrically neutral solution in their central part, which can be considered as a separate phase ( Figure 3). The properties of this solution (dielectric permittivity, ion and water diffusion coefficients) are very close to the corresponding properties of the external free solution [5,29,44]. Note that heterogeneous IEMs contain macropores, while homogeneous IEMs do not [37]. Therefore, one can expect that the volume fraction of electroneutral solution in the homogeneous membranes is essentially lower than in the heterogeneous membranes. As soon as the first electroneutral phase is defined, the remaining volume can be attributed to the second phase.
This phase, called "gel phase" [63,124], includes the polymer matrix with fixed charged groups together with mobile ions and water within the EDLs, which compensate the charge of the fixed ions. Since the central parts of the relatively large pores are excluded from the gel phase, it can be considered as a microporous relatively homogeneous medium. Sometimes, within the gel phase, one distinguishes the third phase, being the aggregate of interlaced hydrophobic polymer chains, which is not permeable for ions and solvent [63]. This hydrophobic phase may also involve the (non-conducting) inert reinforcing binder, such as polyethylene in heterogeneous membranes.
There is no distinct boundary between the inter-gel electroneutral solution and the EDL, which belongs to the gel phase. Nevertheless, the separation of the membrane material into different phases is a useful method, which allows one to simplify the mathematical description. Mafé et al. [125] carried out such a separation when applying a continuous space charge model for calculation of the membrane conductivity.
The main idea for modelling the transport within a membrane considered as a multiphase medium is in attributing some physico-chemical properties to each phase and then to describe the properties characterizing the whole membrane as functions of the single phase properties. This idea is the central one in the approach called effective medium theory (or model) [91]. Let us consider a macroscopic volume in the form of a layer of thickness dx. This layer should be sufficiently large to be "representative" and contain all phases of the membrane ( Figure 3). Simultaneously, it should be small enough to allow differentiation. We neglect the detailed variations of concentrations and potential within dx and operate with the averaged quantities in the phase elements; the phases are assumed in equilibrium with one another within the dx layer.
When ion and water transfer is described in accordance with irreversible thermodynamics, by Equations (1) or (17), the problem is to obtain the effective transport coefficient L i , which characterizes a membrane layer of thickness dx as a function of coefficients L j i (characterizing an individual phase j), as well as the structural and geometric parameters, describing the shape and mutual position of the phases. This formulation addresses the effective-medium approach [91,126] intensively developed in relation to a great variety of systems and transfer phenomena: electrical conduction [127], diffusion [128], heat transfer [129], optics and other [130,131]. This approach was firstly applied by Maxwell [132], who, as early as in 1873, has theoretically found the conductivity of a system containing a conductive continuum medium and dispersed spheres, the conductivity of which differed from that of the medium. Later, Rayleigh, Lichtenecker, Bruggeman, Landau and Lifshitz, Mackie and Meares [133] and other [134,135] contributed to the development of effective-medium approach. The hypothesis called "principle of generalized conductivity" [135] is often applied for generalization of experimental results. According to this hypothesis, the function relating L i to L j i does not depend on both the nature of applied force and the nature of the transported species.
Regarding the application of this approach to IEM, several forms of functions were proposed, linking the effective conductivity of the medium with that of individual phases. Some of them relate to the diffusion [133,134], the others, the electrical conductance [124,134,136,137]. Gnusin et al. [63,74,135] developed a comprehensive model called the "microheterogeneous model" which deals with conductance coefficients L i and L j i instead of particular transport coefficients. Along with the membrane electrical conductivity [54,[138][139][140][141][142][143], this approach also allows one to find electrolyte diffusion permeability [54,114], permselectivity (ion transport numbers) [54,144] and some other properties [74,[145][146][147]. In the case of electrolyte sorption, conductivity, diffusion permeability and transport numbers, it is possible to quantitatively describe all these four characteristics as functions of bathing solution concentration starting from a single set of structural and kinetic parameters [147].
In the framework of the microheterogeneous model, the effective membrane conductance coefficient, L i , is expressed through L j i in the following form [63,135]: where L g i refers to the gel phase and L s i to the inter-gel electroneutral solution; f 1 and f 2 are the volume fractions of the corresponding phases: and V m are the volumes of gel, inter-gel solution and the membrane, respectively, see Figure 3); α is the structural parameter depending on the position of the phases with respect to the axis of transport: when the phases are parallel to this axis, α = 1; when they are in serial disposition, α = −1; in other cases −1 < α < 1.
When describing the transport in the gel phase, all the assumptions made in the TMS model are accepted. The properties of the inter-gel solution are assumed the same as those of the bathing equilibrium solution; the Nernst-Planck and the LEN assumption are applied there. in the corresponding phase [54,63]: Concentrations c s i and c g i are linked by the Donnan relations, Equation (24), the local equilibrium being assumed between both phases.
Usually the concentration of fixed ions (IEC, Q) in ion-exchange materials is rather high, close to 1 mol (L wet membrane) −1 or higher. Co-ions are therefore strongly excluded, so that their concentration in the gel phase is small compared to Q, at least in relatively dilute solutions. When assuming c − << Q and c + ≈ Q, Equation (24) can be simplified. In the case of a symmetrical electrolyte (z + = −z − = z) one gets [46,148]: where subscript "1" is used for the counterion; "A," for the co-ion; c A refers to the co-ion concentration in the inter-gel solution (assumed to be the same as the external solution). The ions are present in both phases; the partition coefficient, K s , involving the co-ion concentration in the membrane, c * A , can be found using Equation (31) [63]: in which the first term on the right-hand side stands for the contribution of the gel phase and the second term, for that of the inter-gel solution. Q g and K D are the IEC of the gel phase and the Donnan coefficient, respectively. Despite of small value of f 2 (typically less than 0.1 in homogeneous membranes and close to 0.2 in heterogeneous ones), the sorption of electrolyte by the inter-gel spaces is dominant, especially in diluted solutions, because of the co-ion exclusion from the gel involving EDL in micro-and mesopores [46,148,149] is small. For conventional IEMs, Equation (31) is verified for external concentrations up to 1−2 M. Note that an equation similar to Equation (32) was proposed by Geise et al. [150], where the second term on the right-hand side of this equation was attributed to an "uncharged material," which can sorb electrolyte stronger than the charged gel. Equation (32) predicts a linear increase of K s with increasing the external concentration. At low external electrolyte concentration K s is close to f 2 , that is, it value is in the vicinity of 0.1 (or lower) for the homogeneous membranes and 0.2 for the heterogeneous ones [46,148,149].
To summarize, Equations (29)-(32) form the microheterogeneous model. There are six parameters (in the case of a single electrolyte): two static, K D and Q g (thermodynamics coefficients); two structural, f 1 and α; and two kinetic ones, the diffusion coefficients of counterion, D g 1 and co-ion, D g A , in the gel phase. The diffusion coefficients in the inter-gel solution, D s i , are assumed the same as in free solution. When these parameters are known, L i coefficients can be calculated. Then the transport membrane characteristics κ, t i and P, can be found by using Equations (8)- (10) giving the links between the Onsager and Kedem-Katchalsky conductance coefficients. Equations (8)-(10) can be applied not only in the case, where the membrane is equilibrated with a bathing solution but also when there is a concentration gradient across the membrane. Then Equations (8)- (10) are applied at any coordinate in the membrane for a local concentration c of the inter-gel solution. In this way, it is possible to incorporate the microheterogeneous model into a boundary-value problem for modelling ion and solvent transport in membrane systems [46,151].
To obtain the microheterogeneous model's parameters, some experiments should be carried out and the data obtained processed. To get all the six parameters listed above, the IEC, electrolyte uptake, conductivity and diffusion permeability should be determined as functions of the electrolyte concentration. The data processing algorithm is described in References [54,152]. However, often only some among the six parameters are of interest. For instance, the volume fraction of the gel phase (or inter-gel solution) is a quite important membrane characteristic. It may be relatively easily obtained from the concentration dependence of the membrane conductivity (κ). Equations (29) and (8) allows a simple approximation at α → 0 [63]: Since the conductivity of gel, κ g , only slightly depends on electrolyte concentration (since the co-ion uptake by the gel is quite low-when the external concentration is not high, usually <1 M), the lgκ − lgκ s correlation, according to Equation (33), should be linear with f 2 as the coefficient. Numerical calculations indicate that near the "isoconductance point" (where κ = κ g = κ s ), κ is almost independent on α [63]. In the range 0.1c iso < c < 10c iso , the lgκ − lgκ s dependence can be approximated by a straight line with slope f 2 , if |α| ≤ 0.2. For most IEMs, α lays in the range 0.1−0.3 [37,54], hence, Equation (33) can be confidently applied, when the concentration is sufficiently close to the isoconductance concentration c iso . Figure 4 shows Lgκ − Lgc dependencies for different IEMs taken from different papers. Since κ s is approximately proportional to the concentration of bathing solution, c, the slope of these curves should be close to f 2 . Taking into account that some approximations have been made, we will call the value of f 2 found as the slope of Lgκ − Lgc linear regression the "apparent volume fraction of the inter-gel solution, f 2app ". To obtain the microheterogeneous model's parameters, some experiments should be carried out and the data obtained processed. To get all the six parameters listed above, the IEC, electrolyte uptake, conductivity and diffusion permeability should be determined as functions of the electrolyte concentration. The data processing algorithm is described in References [54,152]. However, often only some among the six parameters are of interest. For instance, the volume fraction of the gel phase (or inter-gel solution) is a quite important membrane characteristic. It may be relatively easily obtained from the concentration dependence of the membrane conductivity (κ). Equations (29) and (8) Since the conductivity of gel, κ g , only slightly depends on electrolyte concentration (since the coion uptake by the gel is quite low-when the external concentration is not high, usually <1 M), the lgκ − lgκ s correlation, according to Equation (33) (33) can be confidently applied, when the concentration is sufficiently close to the isoconductance concentration iso c . Figure 4 shows Lgκ − Lgс dependencies for different IEMs taken from different papers. Since s κ is approximately proportional to the concentration of bathing solution, c, the slope of these curves should be close to f2. Taking into account that some approximations have been made, we will call the value of f2 found as the slope of Lgκ − Lgс linear regression the "apparent volume fraction of the inter-gel solution, f2app".
It can be seen from Figure 4 that the slope of linear Lgκ − Lgс approximations for the homogeneous membranes is significantly lower than this slope for the heterogeneous ones. The membranes are identified in Table 1, which contains also the values of f2app and the conductivity of the gel phase, g κ . To a first approximation (justified when the contribution of co-ions in the conductivity of the gel can be neglected), the value of g κ may be considered as a constant equal to the conductivity of the membrane at the isoconductance point.  It can be seen from Figure 4 that the slope of linear Lgκ − Lgc approximations for the homogeneous membranes is significantly lower than this slope for the heterogeneous ones. The membranes are identified in Table 1, which contains also the values of f 2app and the conductivity of the gel phase, κ g . To a first approximation (justified when the contribution of co-ions in the conductivity of the gel can be neglected), the value of κ g may be considered as a constant equal to the conductivity of the membrane at the isoconductance point. AMV anion-exchange
Zabolotsky et al. [41,148] have compared the values of f 2 found from conductivity measurements when applying Equation (33) and when using other methods, namely electrolyte sorption, differential scanning calorimetry (DSC) and standard contact porosimetry [37,38,156]. Relatively close values of f 2 were reported. For example, for a heterogeneous MK-40 membrane, the authors found f 2 equal to 0.17 ± 0.02 from conductivity measurements [41], 0.10 ± 0.02 [148] and 0.26 ± 0.02 [41] from sorption, 0.21 from DSC [41] and 0.23 ± 0.05 from porosimetry [41]. These results show that the microheterogeneous model adequately describes the membrane properties when operating with reasonable structure and kinetic parameters. The ease of applying Equation (33) for processing experimental conductivity data, the possibility to determine the parameters having a clear sense, the ability to describe a set of transport properties, involving the diffusion permeability and ion transport numbers, makes this model suitable for IEMs characterization [8,37,138,[142][143][144][157][158][159][160][161].

Modelling of Properties of Ion-Exchange (IEMs) Containing Nanoparticles
Recently, IEMs containing immobilized nanoparticles have attracted a keen interest in the literature. This type of membrane modification was found promising for improving IEMs used in fuel cells. It was established that Nafion based composite membranes with incorporated SiO 2 , ZrO 2 , TiO 2 and other nanoparticles showed better water retention and proton conductivity, especially at low water vapor pressure, in a rather large temperature interval up to 140 • C [162][163][164][165][166][167], which sometimes was accompanied with increasing IEC [168].
Yaroslavtsev et al. [169][170][171] reported that along with increasing water content and conductivity, the hybrid nanoparticle-containing membranes showed a decrease in diffusion permeability and an increase in permselectivity.
It has been shown that IEMs modification with different nanoparticles is very important not only for a better performance of fuel cells but also for improving separation processes in different industries [172]. Namely, the modification of RO, NF and ultrafiltration membranes by immobilization of inorganic nanoparticles increases their antifouling resistance [172].
Yaroslavtsev [169] suggested an interpretation of the change of membrane transport characteristics, in particular, the increase in κ by proposing a physicochemical model named the model of semielasticity of membrane pore walls. The main impact of nanoparticles on the membrane properties, according to this model, is due to the expansion of the pores and channels caused by the nanoparticles. The increase in the size of membrane channels leads to growing ion conductivity. Nevertheless, when the volume fraction of nanoparticles is too large, they block the ion conducting channels, which results in a decrease of the conductivity.
Quantitative description of the effect of nanoparticles on the transport properties of IEMs is a rather complicated problem. Filippov et al. [115][116][117] developed a fine-pore model of ion transfer in an IEM with immobilized nanoparticles. The membrane is considered as quasi-homogeneous medium as in the TMS theory and applies the approach known as the "fine porous membrane" model suggested by Starov, Churaev and Martynov [173]. The simulation was compared with experimental data on diffusion permeability of a hybrid MF-4SK membrane containing silica nanoparticles. This comparison allowed to the authors [115][116][117] to find in what extent the ion effective diffusion coefficients change in the membrane due to incorporation of silica. By applying the "fine porous membrane" model, the main membrane physicochemical parameters can be determined. In addition to the diffusion coefficients, it is possible to find the equilibrium distribution coefficients of electrolyte [117].
Porozhnyy et al. [174] proposed a mathematical model describing the effect of nanoparticles on the transport properties of membranes on the basis of the microheterogeneous model [54,63]. In addition to two phases (the gel phase and the electroneutral inter-gel solution) present in commercial membranes, there are nanoparticles. The core of a nanoparticle is non-permeable to ions and water. However, the EDL, which forms around the particle in the case where its surface is charged (Figure 5), contributes to the selective permeability of the membrane.
Porozhnyy et al. [174] proposed a mathematical model describing the effect of nanoparticles on the transport properties of membranes on the basis of the microheterogeneous model [54,63]. In addition to two phases (the gel phase and the electroneutral inter-gel solution) present in commercial membranes, there are nanoparticles. The core of a nanoparticle is non-permeable to ions and water. However, the EDL, which forms around the particle in the case where its surface is charged ( Figure  5), contributes to the selective permeability of the membrane. It is assumed that the nanoparticles are localized in the meso-and macropores; they replace there the bulk-like electroneutral solution. If the surface of nanoparticle is charged (e.g., in the case of SiO2 or ZrO2, where it has a negative charge), the EDL around the particle is considered to be complementary to the EDL occurring at the charged pore walls ( Figure 5).
To take into account the presence of nanoparticles in a pore, the inter-gel domain is considered as a "two-phase" system containing a solution with the volume fraction = and a solid nonconductive inclusion (the core of the nanoparticle) with the volume fraction = , where + = 1 . The effective conductance of the inter-gel domain containing a solution and a nanoparticle, , may be presented, in accordance to the approach expressed by Equation (29), as It is assumed that the nanoparticles are localized in the meso-and macropores; they replace there the bulk-like electroneutral solution. If the surface of nanoparticle is charged (e.g., in the case of SiO 2 or ZrO 2 , where it has a negative charge), the EDL around the particle is considered to be complementary to the EDL occurring at the charged pore walls ( Figure 5).
To take into account the presence of nanoparticles in a pore, the inter-gel domain is considered as a "two-phase" system containing a solution with the volume fraction f sin = V sol V intergel and a solid nonconductive inclusion (the core of the nanoparticle) with the volume fraction f pin = V particle V intergel , where f sin + f pin = 1. The effective conductance of the inter-gel domain containing a solution and a nanoparticle, L s+p i , may be presented, in accordance to the approach expressed by Equation (29), as where the particle core conductance L p i = 0 (the non-conductive medium), β is a structural parameter, whose role is similar to that of α in Equation (29).
L s1 i characterizes the solution within a pore, which can be considered again as a two-phase domain, containing the bulk electroneutral solution and the EDL, formed by the nanoparticle. L s1 i is found by using the same two-phase approach: where f EDLp = V EDLp V sol+EDLp is the volume fraction of the nanoparticle's EDL within this domain, 1 − f EDLp is the volume fraction of the electroneutral solution whose conductivity is, as before, assumed equal to that of the free solution.
To find V EDLp , it is assumed that each nanoparticle is a sphere with diameter d; the size of all particles is the same. The growth of their volume fraction is achieved by increasing the number of pores where they are present.
The thickness of the diffuse part of the EDL is presented as follows: where ε is the relative dielectric permittivity, ε 0 is the vacuum permittivity, z i is the ion charge number, F is the Faraday constant. It is assumed that the diffusion coefficients in the nanoparticle's EDL are the same as in the free solution, while the concentrations of cation and anion are found using the Donnan relation, as in the gel phase.
When applying Equations (34)- (36), it is possible to calculate the conductance coefficients L i and then the membrane conductivity, electrolyte diffusion permeability, transport numbers.
The EDL formed around a nanoparticle plays an important role. It contributes to an increase in the membrane counterion conductivity, since the counterion concentration there is essentially higher than in the bulk pore solution. Increasing counterion conductivity and reduction of the concentration of co-ions (which are excluded from the EDL) results in higher membrane permselectivity. As the EDL thickness increases with decreasing the salt concentration in the external solution, its role increases; hence, the difference in membrane conductivity and permselectivity between commercial membranes and the membranes doped with nanoparticles increases with decreasing the concentration of the external solution.
However, the solid core of the nanoparticle is an insulator and this contributes to loss of the membrane conductivity at relatively high concentrations (>1 M), where the EDL thickness and, consequently, its role decrease. Thus, according to the model [174], the conductivity of nanoparticles containing membranes should be higher than that of the pristine membranes in dilute solutions. On the contrary, their conductivity should be lower than that of the pristine membranes in concentrated solutions. However, the factor of increasing pore size and water content (determined experimentally and not taken into account in the model [174]) could contribute in additional increase in conductivity. Though, as far as we know, there are no publications reporting experimental data suitable for comparison of the membrane conductivity at different concentration, in the cases where the nanoparticles are present or absent. Known publications [166, 170,175] provides the conductivity of the membrane equilibrated in pure water. As for the membrane permselectivity (quantified by the counterion transport number, t 1 ), this parameter should be higher in the membranes doped with nanoparticles, since even when the EDL thickness is low, the increase in t 1 would be conditioned by the fact that a nonconductive nanoparticle replaces the electroneutral solution present in the inter-gel spaces. Therefore, contribution of the non-selective ion transport in this solution will be reduced in favour of highly selective transport in the charged gel phase.
The content (volume fraction) of nanoparticles in the membrane, f p , is an important parameter. With increasing f p , the membrane conductivity initially increases (because the total EDL volume increases), then it passes through a maximum and declines. The latter is due to the saturation of the pore space available for nanoparticles. When f p becomes large enough, the nanoparticles block ion conductivity in the pores, which they occupy. Then ion transport occurs only across the microporous medium (the "gel phase"), which is assumed not available for nanoparticles. An example of such a dependence is given in Figure 6a. It involves experimental data for a Nafion ® membrane containing nanoparticles of caesium acid salt of phosphotungstic heteropolyacid [170] and the results of calculation according to the model described above [174]. Note that a similar dependence was found for a perfluorinated sulfocation-exchange MF-4SK membrane (analogue of Nafion ® ): the proton conductivity curve as a function of additive content has a maximum at 2.6 vol%, when it is doped with SiO 2 and 1.5 vol%, when it is doped with PAni [175]. Parizian et al. [176] studied an experimental IEM modified with ZnO nanoparticles. Similarly as Safronova et al. [170] and Novikova et al. [175], they found that the dependence of water content, membrane potential, transport number and selectivity of this membrane initially increased with increasing the content of added ZnO nanoparticles, then these characteristics reach their maximum and declined when the nanoparticles' content was rather large.
The presence of nanoparticles within the pores reduces their permeability towards the electrolyte diffusion (Figure 6b). Both the particle's core and the EDL around it contributes to a decrease in the diffusion permeability. However, there is another factor, the expansion of pores (not taken into account in the model [174]), which gives rise to increasing membrane diffusion permeability. Which factor is determinant, depends on the membrane material (which can be more or less rigid) and the incorporated nanoparticles. Filippov, Safronova and Yaroslavtsev [115] reports an increase in the diffusion permeability of a perfluorinated sulfocation-exchange MF-4SK membrane after its doping with silica nanoparticles. Similar data are provided in Refs. [168,175].
An increase in conductivity and a decrease in diffusion permeability essentially improve the permselectivity of modified membranes towards counterions [170,174].
content was rather large.
The presence of nanoparticles within the pores reduces their permeability towards the electrolyte diffusion (Figure 6b). Both the particle's core and the EDL around it contributes to a decrease in the diffusion permeability. However, there is another factor, the expansion of pores (not taken into account in the model [174]), which gives rise to increasing membrane diffusion permeability. Which factor is determinant, depends on the membrane material (which can be more or less rigid) and the incorporated nanoparticles. Filippov, Safronova and Yaroslavtsev [115] reports an increase in the diffusion permeability of a perfluorinated sulfocation-exchange MF-4SK membrane after its doping with silica nanoparticles. Similar data are provided in Refs. [168,175].  [170]; the curves are calculated by applying the model [174] for different particle diameters shown in nm near the curves. The best fitting is obtained, when assuming that the nanoparticle size is 16 nm. Reproduced with permission from [174], Elsevier, 2016.
An increase in conductivity and a decrease in diffusion permeability essentially improve the permselectivity of modified membranes towards counterions [170,174].

"Pore-flow" Models
The "pore-flow" models consider the transport of species as occurring within the membrane pores forming ion and water conducting channels [85,86]. This approach originates from the description of ion and water transport inside microcapillaries with charged walls [90,[177][178][179][180]. Namely, the main goal of these models is the description of electrokinetic phenomena, that is, coupled transfer of a liquid relative to a charged solid (or liquid) surface or the transfer а charged solid particles relative to a liquid [181,182]. The main role in these phenomena is played by the EDL on charged surfaces.
The term "space-charge models" [183][184][185][186][187][188][189][190] is often applied in the literature for this kind of models. In these models, the ion and water transport is described on the basis of the Nernst-Planck-  [170]; the curves are calculated by applying the model [174] for different particle diameters shown in nm near the curves. The best fitting is obtained, when assuming that the nanoparticle size is 16 nm. Reproduced with permission from [174], Elsevier, 2016.

"Pore-flow" Models
The "pore-flow" models consider the transport of species as occurring within the membrane pores forming ion and water conducting channels [85,86]. This approach originates from the description of ion and water transport inside microcapillaries with charged walls [90,[177][178][179][180]. Namely, the main goal of these models is the description of electrokinetic phenomena, that is, coupled transfer of a liquid relative to a charged solid (or liquid) surface or the transfer a charged solid particles relative to a liquid [181,182]. The main role in these phenomena is played by the EDL on charged surfaces.
The term "space-charge models" [183][184][185][186][187][188][189][190] is often applied in the literature for this kind of models. In these models, the ion and water transport is described on the basis of the Nernst-Planck-Poisson-Navier-Stokes equations. The main parameters characterizing the membrane are pore size and wall charge density [191].
The basic version of a space-charge model deals with a pore (schematically shown in Figure 7). It involves the extended Nernst-Planck equation containing the convective term, Equation (22) written in 2D geometry. The diffusion coefficients are usually assumed to be the same as in free solution. The cross-sectional distribution of local concentration, c i , is described by the Poisson-Boltzman equation. The fluid flow (assumed to occur only in the axial direction) is described by the Navier-Stokes equation, where the axial pressure gradient and the electrical body force provided by the tangential electric field applied to the space charge in EDL are taking into account [177,183]. The most studies dealt with the pores whose walls are uniformly charged, such as shown in Figure 7, bottom. The impact of a non-uniform distribution of the pore wall charge on NF performance has been investigated by the group of Szymczyk in Refs. [186,192,193].
The research field making accent on the fluid transfer in narrow channels was only recently termed as nanofluidics but it has a quite rich history, in particular, in membrane science [194]. Nanoand microfluidics research as well as the applications in this field have experienced extensive growth in recent years [195][196][197][198][199]. 100 nm enables the occurrence of the phenomena, which are impossible at bigger length scales.
Boltzman equation. The fluid flow (assumed to occur only in the axial direction) is described by the Navier-Stokes equation, where the axial pressure gradient and the electrical body force provided by the tangential electric field applied to the space charge in EDL are taking into account [177,183]. The most studies dealt with the pores whose walls are uniformly charged, such as shown in Figure 7, bottom. The impact of a non-uniform distribution of the pore wall charge on NF performance has been investigated by the group of Szymczyk in Refs. [186,192,193]. The research field making accent on the fluid transfer in narrow channels was only recently termed as nanofluidics but it has a quite rich history, in particular, in membrane science [194]. Nanoand microfluidics research as well as the applications in this field have experienced extensive growth in recent years [195][196][197][198][199]. 100 nm enables the occurrence of the phenomena, which are impossible at bigger length scales.
The "pore-flow" models involve the effects, which are specific for interfaces: change of the dielectric constant as a function of the distance from the pore wall [183], finite ion sizes, ion hydration effects [200], adsorption [185] and other [196]. These effects are taken into account when computing the streaming potential [201], diffusion permeability [186] and permselectivity [202], pore conductivity [180,203]. This type of models has been shown as very useful for describing the ion and fluid transport in and around nanometre-sized objects [204]. Bazant et al. [182,205] paid attention to the steric effects near a non-permeable wall at high voltages taking into account that solvated counterions are crowded there.
Cwirko and Carbonell [206] and later on Szymczyk et al. [207] have calculated using spacecharge models the macroscopic Onsager's Lij coefficients for a Nafion membrane as functions of the membrane nanostructure parameters (which are the pore radius, pore wall charge density and the tortuosity factor). Their comparison with the experimental coefficients determined by Narebska et al. [65] has shown a rather good agreement. Thus, it becomes possible to bridge the gap between two different approaches, that is, the microscopic model description and the irreversible thermodynamics. The "pore-flow" models involve the effects, which are specific for interfaces: change of the dielectric constant as a function of the distance from the pore wall [183], finite ion sizes, ion hydration effects [200], adsorption [185] and other [196]. These effects are taken into account when computing the streaming potential [201], diffusion permeability [186] and permselectivity [202], pore conductivity [180,203]. This type of models has been shown as very useful for describing the ion and fluid transport in and around nanometre-sized objects [204]. Bazant et al. [182,205] paid attention to the steric effects near a non-permeable wall at high voltages taking into account that solvated counter-ions are crowded there.
Cwirko and Carbonell [206] and later on Szymczyk et al. [207] have calculated using space-charge models the macroscopic Onsager's L ij coefficients for a Nafion membrane as functions of the membrane nanostructure parameters (which are the pore radius, pore wall charge density and the tortuosity factor). Their comparison with the experimental coefficients determined by Narebska et al. [65] has shown a rather good agreement. Thus, it becomes possible to bridge the gap between two different approaches, that is, the microscopic model description and the irreversible thermodynamics.

Current-Induced Concentration Gradients
The permselectivity of IEMs towards the ions bearing the charge of a certain sign (counterions) is the main functional property [208]. Namely this property conditions the applications of this kind of membranes in ED, where the feed electrolyte solution is simultaneously demineralized and concentrated in different ED compartments (Figure 8).
When a direct current of density i flows across a membrane, concentration polarization occurs: concentration gradients arise in the boundary solution layers adjacent to the membrane. Immediately after the current is turned on, there are no diffusion fluxes, while the counterion migration flux through the membrane is higher than their flux from the solution bulk to the membrane surface. As a result, the concentration of the counterions (as well as that of co-ions) decreases at one membrane side and increases at the other. The changes of concentration near the membrane continue until the diffusion flux, which grows with time, completely compensates for the difference in migration fluxes in solution and membrane.

Current-Induced Concentration Gradients
The permselectivity of IEMs towards the ions bearing the charge of a certain sign (counterions) is the main functional property [208]. Namely this property conditions the applications of this kind of membranes in ED, where the feed electrolyte solution is simultaneously demineralized and concentrated in different ED compartments (Figure 8). When a direct current of density i flows across a membrane, concentration polarization occurs: concentration gradients arise in the boundary solution layers adjacent to the membrane. Immediately after the current is turned on, there are no diffusion fluxes, while the counterion migration flux through the membrane is higher than their flux from the solution bulk to the membrane surface. As a result, the concentration of the counterions (as well as that of co-ions) decreases at one membrane side and increases at the other. The changes of concentration near the membrane continue until the diffusion flux, which grows with time, completely compensates for the difference in migration fluxes in solution and membrane.
At any time, the flux density of species i coming from the solution and entering the membrane/solution interface is equal to the flux density leaving the interface into the membrane bulk. This condition of flux continuity is written as: where ( ) i s J is the flux density of ion i through the interface; the middle part of the Equation (37) refers to the solution: D is the electrolyte diffusion coefficient in solution, Ci and ti are the concentration and transport number of ion i in solution; the right-hand part refers to the membrane: Ti is called integral [31,209] or effective [54,210] transport number of this ion in the membrane. Ti is defined as the current fraction carried by ion i through the interface or the membrane in steady state where no restrictions on the driving forces are imposed: At any time, the flux density of species i coming from the solution and entering the membrane/solution interface is equal to the flux density leaving the interface into the membrane bulk. This condition of flux continuity is written as: where (J i ) s is the flux density of ion i through the interface; the middle part of the Equation (37) refers to the solution: D is the electrolyte diffusion coefficient in solution, C i and t i are the concentration and transport number of ion i in solution; the right-hand part refers to the membrane: T i is called integral [31,209] or effective [54,210] transport number of this ion in the membrane. T i is defined as the current fraction carried by ion i through the interface or the membrane in steady state where no restrictions on the driving forces are imposed: T i can differ from the (electromigration) transport number in the membrane, t mb i (t i in Equation (5) written for a membrane), since according to Equation (5), the (electromigration) transport number must be determined under the condition that the concentration gradient and the convective flow across the membrane are zero. Nevertheless, if the external solution concentration is not too high, T i is quite close to t mb i . Normally, commercial membranes are highly permselective to counterions [76]: the values of T 1 and t mb i approach 1 when the solution is diluted. According to classical electrochemistry, the occurrence of concentration gradients near the surface of a membrane (or an electrode) results from the limitation of the current density i. With increasing i, the electrolyte concentration at the interface, c s , decreases (Figure 8). When c s approaches zero (becomes much smaller than the bulk concentration, c 0 ), the current density reaches its limiting value (i lim ). The expression for i lim can be obtained from Equations (37), (38), when writing the concentration gradient of the counterion (1) at the membrane surface as: where δ is the thickness of the Nernst DBL in the depleted solution. According to Equation (39), δ is the distance from the membrane to the intersection point of the tangent drawn to the concentration profile at the interface with the line corresponding to the bulk solution concentration. Substituting Equation (39) into Equation (37) and setting c 1s = 0 yields: where c 0 = z + c 0 + = −z − c 0 − is the bulk electrolyte concentration in eq L −1 . Equation (40) was first obtained by Peers in 1956 [211].
According to the classical theory [57] based on the LEN assumption, when i tends to its limiting value (i lim ), the PD over a membrane surrounded by two DBLs tends to infinity. However, in real membrane or electrode systems, the limiting current density can be exceeded in several times [212][213][214][215]. The causes of the overlimiting current are vividly discussed in literature [212,214,[216][217][218][219]. The main effects are current-induced convection (electro-and gravitational convection) and generation of additional charge carriers (e.g., water splitting at the depleted solution/membrane interface, which produces the H + and OH − ions). While current-induced convection enhances the mass transfer across the membrane, water splitting in most cases is undesirable effect causing pH changes and the loss in current efficiency.
Along with these two main effects, additional effects, which may occur in some conditions, are discussed in literature. Andersen et al. [234] described theoretically the effect of current-induced membrane discharge. This effect can occur in a system with an AEM, which has weakly basic functional groups. In this case, OH − ions generated in water splitting at the depleted membrane interface may deprotonate these functional groups, that is, discharge them. This produces: (1) a decrease in the surface charge density, which reduces electroconvection developing as electroosmosis of both the 1st and 2nd kind [235,236] and (2) a decrease in IEC, which results in increasing co-ion back diffusion through the membrane. The latter leads to a decrease in T 1 , hence, in an increase of the limiting current density according to Equation (40). Note that electroosmosis of the 2nd kind is named also induced-charge electroosmosis [237,238], since this effect is due to the extended space charge region [239], which emerges next to the equilibrium EDL at overlimiting currents.
In the early period of ED, the use of underlimiting current was prescribed [210,240]. However, nowadays the overlimiting current regimes are widely used in ED and especially in electro-deionisation of diluted electrolyte solutions [212,241]. The use of such intensive current regimes allows reduction of expensive membranes' area, which can cause the reduction of overall costs of ED desalination [210]. The effect of electroconvection is relevant also in a number of nano-and microfluidic devices [194,242] such as electrokinetic micropumps [243,244]. Besides, intensive current regimes are applied in electrophoresis [242], electrodeposition [245], layering of colloid crystals on electrode surfaces [246] and other domains related to separation science [58,203].

Governing Equations and General Boundary Conditions
The first adequate 2D model of ED (named the hydrodynamic theory of desalination) was developed by Sonin and Probstein [247]. Later on other versions of this model was proposed, in particular, by Gnusin et al. [248,249] and Tedesco et al. [250], who reported a numerical solution of the problem formulated by Sonin and Probstein, when taking into account the difference in the diffusion coefficients of cation and anion in solution; Shaposhnik et al. developed an analytical solution for special cases [251] and accounted for the effect of ion-conducting spacer in numerical solution [252].
Solan et al. [253,254] proposed a "mesh step" model for the effect of a net type spacer-turbulence promoter-on mass transfer in an ED channel. The role of spacer in ED desalination was also studied by a number of authors [255][256][257][258][259][260][261]. Other aspects of modelling of ED are considered in Refs. [262,263] and other contributions are reviewed by Campione et al. [264].
The basic version of 2D mathematical model of ED involves [247][248][249] the extended Nernst-Planck equation (Equation (22)) in vector form, the equation expressing the current density as a sum of ion flux densities (Equation (23)), the continuity equation for the conservation of species i written for the steady state) and the LEN condition: The system of coordinates is shown in Figure 8. Equations (41)-(44) are written for a cation and an anion denoted by subscript i as 1 and 2, respectively.
It is assumed that the fluid velocity at the surface is zero (the no-slip condition) [80]: The fluid flow is considered laminar and steady. For such a flow between two parallel plates (without spacer), the solution of the Navier-Stokes equation with the boundary no-slip condition gives a parabolic velocity distribution along the transverse coordinate known as the Poiseuille distribution, the normal component of flow velocity is zero: where V is the average flow velocity between the plates. The continuity of the ion fluxes through the interface (Equation (37)), which is rewritten in the form ∂c ∂x x=0 is used as the boundary condition for both AEM and CEM. Note that the effective transport numbers of cations in the CEM (T c 1 ) and anions in the AEM (T a 2 ) can take into account the co-ion or H + (OH − ) transport across the membrane. For this, experimental [76,265] or theoretical [266] evaluation of T c 1 and T a 2 should be carried out. The feed solution concentration at the channel entrance is set as the initial condition: The condition specifying the parameters of the operating mode assumes that the PD ∆φ between two planes passing through the centre of two nearest concentration compartments (CC) is known. This PD can be found as ∆φ =(U − ∆φ el )/n, where U is the voltage measured between the polarizing electrodes of the ED cell, ∆φ el is the PD on the electrodes and n is the number of cell pairs. It is assumed that any plane passing through the centre of an CC is equipotential and the current density and ionic fluxes along this plane are zero.
The PD per a cell pair (considered as a known parameter) is determined by the equation: where R a , R c and R cc are the resistances of the AEM, CEM and the concentration compartment, respectively (all are assumed known), ∆ϕ D is the sum of (Donnan) PDs over the solution/membrane interfaces for AEM and CEM, E = − ∂ϕ ∂x is the electric field; h 0 Edx is the PD in the DC. The diffusion transport in the membranes and the CC is not taken into account, the → i and → E vectors are assumed directed normally to the membrane surface. According to Equation (27), Equations (41)-(49) form a 2D boundary value problem under the condition that E(x,y) is expressed as a function of electrolyte concentration c(x,y) and current density i(x,y).
Taking into account the condition V x = 0, the x-component of Equation (41) is written as the conventional Nernst-Planck equation. The electric field strength can be found from Equations (41), (42) and (43): (51) where c = z 1 , c 1 = −z 2 c 2 is the electrolyte concentration (in eq L −1 ). Substituting Equation (51) into Equation (41) yields the transport equation in the (reduced) form of Kedem-Katchalsky equation (Equation (5)): where D (the electrolyte diffusion coefficient) and t k (the transport number of ion k) are expressed as follows: In the vector form when taking into account the convective term, Equation (52) reads as follows: When considering the normal and longitudinal components of → j i , it is assumed that the normal velocity is zero and the longitudinal diffusion and migration are negligible compared to the convective transport: Then the application of Equation (43) to Equation (55) gives [80]: Here for the sake of simplicity, V(x) refers to the longitudinal velocity. The diffusion coefficient and the transport numbers are assumed constant and known quantities.

Concentration Distribution and Diffusion Layer Thickness
Numerical solution of the convection-diffusion problem presented by Equations (45)-(49) and (58) gives the distribution of concentration c(x,y), current density i(y) and electric potential φ(x,y) in the space between an AEM and CEM under the condition that the overall PD per cell pair ∆φ is known and the geometric parameters of the cell as well as the average flow velocity are set. In the version presented above, the variation of concentration for the sake of simplicity is taken into account only in the DC since namely this compartment determines the behaviour of the system in ED desalination: the concentration variations in the CC only slightly affect the overall system resistance.
The concentration profiles in the DC computed at different distances from the cell entrance are shown in Figure 9.
In the vector form when taking into account the convective term, Equation (52) reads as follows: When considering the normal and longitudinal components of i j  , it is assumed that the normal velocity is zero and the longitudinal diffusion and migration are negligible compared to the convective transport: Then the application of Equation (43) to Equation (55) gives [80]: Here for the sake of simplicity, V(x) refers to the longitudinal velocity. The diffusion coefficient and the transport numbers are assumed constant and known quantities.

Concentration Distribution and Diffusion Layer Thickness
Numerical solution of the convection-diffusion problem presented by Equations (45)-(49) and (58) gives the distribution of concentration c(x,y), current density i(y) and electric potential φ(x,y) in the space between an AEM and CEM under the condition that the overall PD per cell pair Δφ is known and the geometric parameters of the cell as well as the average flow velocity are set. In the version presented above, the variation of concentration for the sake of simplicity is taken into account only in the DC since namely this compartment determines the behaviour of the system in ED desalination: the concentration variations in the CC only slightly affect the overall system resistance.
The concentration profiles in the DC computed at different distances from the cell entrance are shown in Figure 9. It can be seen that the concentration profiles are not symmetric, which is due to the difference in the transport numbers of cation and anion leading to a difference in the concentration gradients at the AEM and CEM, according to Equation (47) used as the boundary condition for the problem. Despite the difference in the concentration gradients, the thickness of the Nernst DBLs (determined through the intersection point of tangents to the concentration profiles, Equation (39)) at a AEM and CEM are the same, and, in the first approximation, do not depend on the applied voltage [248].
When a sufficiently great PD is applied, the concentration at the surface of the membrane, for which the difference of counterion transport numbers (T i − t i ) is maximum, reaches a nearly zero value, the average current density for this membrane (the CEM in the case of NaCl) attains the limiting value, Equation (40). Subsequent growth of ∆φ does not lead to an increase in the current density. When the current density is limiting, the problem (45)- (49) and (58) can be simplified. Instead of the boundary condition of the second kind [Equation (47), expressing the concentration gradient at the membrane surface], the boundary condition of the first kind, the concentration equal to zero at this surface, can be used: This simplified boundary value problem can be solved approximately. Indeed, an asymptotic solution was obtained by Lévêque in the case of the diffusion-convection heat transfer [267]. Adapting this solution to diffusion-convection mass transfer gives [80]: where i lim loc and i lim av are the local and the average over the channel length (L) current densities, respectively; i lim av = 1 L L 0 i lim loc (y)dy, where y is the distance from the channel entrance.
Combining Equations (60) with the Peers equation, Equation (40), yields the following equations for the local, δ N loc and the average, δ N av , values of the Nernst DBL thickness, δ N , in the case where the second term ("0.2") in brackets on the right-hand side is neglected [248,268]: As it can be seen from Equations (60), the relative contribution of the "0.2" term on the right-hand side decreases with decreasing the channel length L. However, similar equations can be obtained by treating the values of i lim or δ N found numerically from a convection-diffusion problem where zero concentration at the wall is used as the boundary condition. In this case, the limiting current density is in excellent agreement with Equations (60), while the values of the local and average δ N are slightly (about 4%) higher than those presented by Equations (61) due to the effect of the omitted second term. Namely, instead of 0.68 and 1.02, the numerical solution gives 0.71 and 1.06, for the local and the average δ N , respectively [269].
Taking into account the definition of i lim av above and the fact that i lim and δ N are inversely proportional to each other, Equation (40) The mass transfer rate in an ED channel at the limiting current can be also expressed in terms of the Sherwood, Reynolds and Schmidt numbers [80,248]: where Sh av is the average value of Sh in a channel of length L, ν is the kinematic viscosity.
For an ED channel of a dimensionless length Y = LD h 2 V 0 = 0.01, the local and average DBL thicknesses reach 0.23h and 0.15h, respectively; the degree of desalination at the limiting current density is 10-13% (it depends on the type of electrolyte). These values are typical for "short" channels in cells commonly used for electrochemical studies of IEMs. For "long" channels, where the DBLs formed at the AEM and CEM are overlapped, it is possible to use numerically calculated Sh vs. Y curves [248]. For rough estimations of the local limiting current density using Equation (40), the effective average value of δ N can be taken equal to h/3.
Equations (60) agree well with the experimentally determined limiting current density in a cell specially designed [270] to satisfy the conditions required for obtaining this equation. The cell has the input and output comb-like devices, which assure the laminar flow for the average velocities from 10 −2 cm/s to 2 cm/s. The intermembrane distance (h = 5−6 mm) and the channel length (20 mm) give the dimensionless length Y in the interval 10 −5 -10 −2 , where the Lévêque approximations (Equations (60)) are well satisfied. The possibility of knowing the theoretical limiting current density and the DBL thickness is very important, as it allows evaluation the coupled effects of concentration polarization. Indeed, the theoretical values of i lim and δ are determined for the case of an "ideal" IEM, which does not generate current-induced convection and does not split water molecules. Hence, the comparison of the experimental and theoretical i lim can give some ideas, for example about the intensity of electroconvection [217]. The role of the Lévêque equations for the slot cells with laminar flow is similar to the Levich equation obtained for the rotating disk electrode [271] or rotating membrane disc [225,272,273].

ED with Heterogeneous Membranes
Competition between homogeneous and heterogeneous IEMs lasts for decades. Although the conductivity and permselectivity of a heterogeneous IEM is generally lower than these parameters of a homogeneous membrane, the great advantage of heterogeneous IEMs is that they are less costly.
The surface of heterogeneous membranes includes conductive and non-conductive areas of the size from several µm to several tens of µm. The so called homogeneous membranes (usually prepared the paste method via copolymerization [9]) also contain non-conductive inclusions (e.g., PVC binder), however, the size of these inclusions generally less than 1 µm (except the threads of backing cloth). For instance, the Neosepta ® AMX and CMX membranes have PVC particles with a diameter of about 100 nm. The presence of non-conductive areas leads to the accumulation of the electric current lines on the conductive regions and to the curvature of these lines. On one hand, this "funnel effect" [274] increases membrane concentration polarization and water splitting rate, while reducing the limiting current density [274]. On the other hand, tangential electric field enhances electroconvection and mass transfer [232,236,275,276]. From what has been said before, the electrical heterogeneity of the membrane surface is subject to optimization. Too low surface fraction of membrane conductive area (γ) will lead to very high concentration polarization, which could be remediated by electroconvection. However, one can imagine that a certain heterogeneity of the surface could effectively enhance electroconvection, while only slightly increase concentration polarization. Simulations made by Davidson et al. [233] and Zabolotsky et al. [272] predict that the rate of electroconvection and, as a consequence, the rate of mass transfer pass through a maximum with increasing γ at a fixed voltage. According to [233], the maximum is close to γ = 0.5; according to [272], it is close to γ = 0.8. The both models are based on the same governing equations (the Nernst-Planck-Poisson-Navier-Stokes equations). Perhaps, the difference is due to the fact that Davidson et al. [233] do not take into account the forced flow between the membranes, while Zabolotsky et al. [272] do.
Another important lever for enhancing electroconvection and mass transfer is geometric heterogeneity (surface undulation, waviness, formation of special patterns). However, this option is still poorly studied, which can be understood in view of the wide variety of possible geometric shapes and sizes. A deeper understanding of the peculiarities of ion and water transport at electrically and geometrically heterogeneous surface is needed for targeted optimization of membrane surface parameters. Later on we will consider some aspects of mathematical modelling of ion transport at heterogeneous membrane surface.

Effect of Electrical Surface Heterogeneity on the Limiting Current Density and Diffusion Layer Thickness
Consider an ED channel, which is formed by a homogeneous AEM (which may have areas on their surface that differ in conductivity but the size of these areas is less than 1 µm) and a heterogeneous CEM whose surface contains alternating conductive and non-conductive areas, whose size is of the order of 10 µm or more. Figure 10 gives a 2D schematic presentation of such a channel. The length of the conductive region on the surface is denoted as 2b, the distance between the centres of two nearest conductive regions (the spatial period of surface heterogeneity) is 2l. As earlier (Section 5.2), we assume that a binary electrolyte solution with initial concentration c 0 flows between the membrane, the flow being laminar and there is Poiseuille velocity distribution, Equation (46). A PD ∆ϕ is applied between two planes passing through the centres of the concentration compartments adjacent to the central DC.
Poisson-Navier-Stokes equations). Perhaps, the difference is due to the fact that Davidson et al. [233] do not take into account the forced flow between the membranes, while Zabolotsky et al. [272] do.
Another important lever for enhancing electroconvection and mass transfer is geometric heterogeneity (surface undulation, waviness, formation of special patterns). However, this option is still poorly studied, which can be understood in view of the wide variety of possible geometric shapes and sizes. A deeper understanding of the peculiarities of ion and water transport at electrically and geometrically heterogeneous surface is needed for targeted optimization of membrane surface parameters. Later on we will consider some aspects of mathematical modelling of ion transport at heterogeneous membrane surface.

Effect of Electrical Surface Heterogeneity on the Limiting Current Density and Diffusion Layer Thickness
Consider an ED channel, which is formed by a homogeneous AEM (which may have areas on their surface that differ in conductivity but the size of these areas is less than 1 μm) and a heterogeneous CEM whose surface contains alternating conductive and non-conductive areas, whose size is of the order of 10 μm or more. Figure 10 gives a 2D schematic presentation of such a channel. The length of the conductive region on the surface is denoted as 2b, the distance between the centres of two nearest conductive regions (the spatial period of surface heterogeneity) is 2l. As earlier (Section 5.2), we assume that a binary electrolyte solution with initial concentration c0 flows between the membrane, the flow being laminar and there is Poiseuille velocity distribution, Equation (46). A PD ϕ Δ is applied between two planes passing through the centres of the concentration compartments adjacent to the central DC. The governing equations used here are the same as in the case of homogeneous membranes: the Nernst-Planck with the convective term, Equation (41) and the continuity equation for the electrolyte species, Equation (43). However, instead of the LEN assumption, Equation (44), we use here the Poisson equation: The governing equations used here are the same as in the case of homogeneous membranes: the Nernst-Planck with the convective term, Equation (41) and the continuity equation for the electrolyte species, Equation (43). However, instead of the LEN assumption, Equation (44), we use here the Poisson equation: where ε 0 and ε are the vacuum permittivity and the relative permittivity of the solution, respectively. As for the boundary conditions, as in Section 5.2, we assume the continuity of the fluxes across both conductive and non-conductive regions. Both CEM and AEM are considered as "ideally selective": the current is carried solely by counterions, hence, at the solution/membrane interface, the counterion partial current density is equal to the total current density and the co-ion normal flux is zero. In addition, following Rubinstein and Shtilman [239], we set the counterion concentration at the conductive surface regions (c 1m at the CEM and c 2m at the AEM, their values are considered as parameters). This condition may be called the Dirichlet-Rubinstein boundary condition [214]. It is applied also by Pham et al. [277] and Kwak et al. [230]. Hence, for the conductive surface regions of the CEM and the whole surface of the AEM, we have: (66) where N a and N c show how many times c 2m and c 1m are higher than the inlet electrolyte concentration (c 0 ), respectively. Subscript "cond" refers to a conductive surface region. The boundary conditions for the non-conductive regions of the CEM set zero normal flux across the surface: When writing ∂c i ∂x + F RT z i c i ∂ϕ ∂x (0, y non−cond ) = 0 for both ions and taking into account that the LEN condition holds in the solution at the insulating surface, one finds that condition (66) leads to zero normal components of both concentration and potential gradients: It is assumed that the electric current does not flow through the channel entrance and the channel exit. At the channel entrance, the electrolyte concentration is given and the LEN condition holds. At the exit, the tangential diffusion is negligible: ions can leave the channel only with the convective flow. This gives the following conditions: As in the case of homogeneous membranes (Section 5.2), the electric mode operating parameter is the PD between two planes in the centres of CC. However, when applying the Poisson equation, the electric potential is set at the surfaces of both IEMs. Similar to Equation (49), the total PD in an ED cell pair (∆φ) is presented as the sum of three contributions: where ∆ϕ DC is the PD within the desalination compartment limited by x = 0 at the left and x = h at the right, ∆ϕ a is the PD across the AEM and a half of the concentration compartment, ∆ϕ c is a similar quantity for the CEM half-cell ( Figure 11).   To find ϕ(0, y) and ϕ(h, y) as functions of the normal local current density (i x ) and ∆φ, it is assumed [214] that between the point referred to the membrane surface facing the CC and the point on the surface facing the DC where the counterion concentration is set (x = 0 for the CEM and x = h for the AEM, Figure 11), there is the membrane bulk and two interfacial quasi-equilibrium EDLs [83,214]. The PD between these two points is the sum of the ohmic PD and two Donnan PDs. For the sake of simplicity, concentration polarization at the membrane interfaces with enriched solution is not taken into account, hence, everywhere in the concentrating compartment the concentration is set equal to c 0 . Under these assumptions, the PD between the plane passing through the centre of the CC and x = 0, ∆ϕ c and the PD between the plane passing through the centre of the CC and x = h, ∆ϕ a , read as follows [214]: where R ms is the total ohmic (membrane+solution) resistance of the system involving an AEM and a CEM as well as a concentration compartment. Equations (71) and (72) are written only for the conductive regions on the CEM surface. The second terms on the left-hand side of Equations (71) and (72) involve the total interfacial PD. Thus, for the CEM, we have Setting φ = 0 in the middle of the left-hand concentration compartment and taking into account that ϕ(0, y) = ∆ϕ c and ∆ϕ = ϕ(h, y) + ∆ϕ a , we find from Equations (71) and (72): The positions of the cross-sections for which the profiles are calculated are shown in Figure 12; the first cross-section is at 550 µm from the entrance; the intermembrane distance h = 1000 µm; the average flow velocity of solution in the channel between the membranes is 0.08 cm·s −1 ; the length of the conductive surface regions is fixed: a = 10 µm.  The part of the concentration profile close to the surface changes essentially when passing from one cross-section to another along the y-axis (Figure 13). At the centre of the conductive region, the normal component of the concentration gradient (the slope of the concentration profile) is maximum, while it is zero at the centre of the non-conductive region. However, at some distance from the surface the changes of the concentration gradient slope with y are rather small: this slope gradually decreases as the distance from the entrance increases. Thus, it is possible to define the effective DBL thickness δ N from the concentration profile slope at a distance d from the surface, where the changes in the slope caused by the surface being conductive or not are negligible: For a homogeneous membrane, d = 0 and Equation (75) gives the quantity, which is identical to the Nernst DBL thickness, as c(d) is the concentration at the membrane surface. Together with the Nernst DBL, the total DBL thickness, δ tot , is often used. δ tot can be defined [268], as the distance from the surface to the point in solution where c = 0.99c 0 . The results of the computations of δ N and δ tot for the cross section at a distance of 550 µm from the ED channel entrance are shown in Table 2.   When a current of density i is applied to an IEM initially in an equilibrium state with the bathing solution, the concentration profiles develop with time [12,31,278]. Electrolyte diffusion from the solution bulk to the membrane surface arise and the diffusion flux of counterions is added to their migration in solution, which slows down the rate of concentration changes at the interface. If the current density is not high (lower than the steady-state limiting current density, i lim ), the electrolyte diffusion reaches such a rate, that the difference between the counterion migration fluxes in the membrane and solution gets compensated. The system reaches a steady state whereas the main mechanism of current transfer remains electrodiffusion. However, if the current density is higher than i lim , at a certain time called "transition time" there is a change in current transfer mechanism across the membrane. At i > i lim , the electrolyte diffusion is not sufficient to stop the decrease of concentration at the membrane interface. When the electrolyte concentration at this interface reaches a rather small critical value [216,277], current-induced mechanisms of overlimiting current conductance such as coupled convection (gravitational convection and electroconvection), water splitting [24,220,279] and others [213,234] occur additionally to the electrodiffusion. Current-induced convection enhances the delivery of electrolyte to the membrane surface [217,230,232]. Since the rate of electrolyte delivery via current-induced convection increases as the current density increases, the system reaches a steady-state where its resistance is higher than that at i < i lim but remains finite.
The onset of these additional transport mechanisms occurs in the vicinity of the transition time τ, which corresponds to the inflection point on the chronopotentiogram. The rate of increase of the PD across the membrane attains its maximum. After reaching the transition state (t = τ), the rate of PD increase slows down. Then the system reaches a steady state where the PD does not change (or slowly changes) with time. The current-induced mechanisms of ion transport result in reducing concentration polarization. Even when the current density is higher than i lim , the ion concentration at the interface never reaches zero and PD is not tends to infinity [239,280].
Chronopotentiometry (CP) is a method of studying the transient behaviour of electrode and membrane systems, when a (constant) current density is applied and the response of time-dependent PD is measured. It is a powerful technique for studying interfacial ion transport and chemical reactions [15,278,281,282].
The transition time is an important characteristic of time-dependent ion transport [15,282,283]. In 1901 Sand [284] deduced a simple equation for the transition time, τ, when considering diffusion in a semi-infinite 1D domain with a constant-flux boundary condition: The adaptation of the Sand equation to IEMs with non-ideal permselectivity expressed by T k was first made apparently by Lerche and Wolf [285] and later by Krol et al. [282].
In Sand's theory, the DBL thickness is assumed infinite; the forced or natural convection occurring in the membrane/electrode system is not taken into account. The transition time is defined as the moment where the interfacial concentration attains zero and the PD tends to infinity.
The fact that the DBL thickness is assumed infinite leads to essential deviation between theoretical and experimental values of in certain conditions. The situation is similar to that in the theory of impedance. Warburg [286] proposed a simple theory for computing the impedance of solution/electrode system as a function of frequency assuming that the DBL thickness is infinitely large. However, taking into account the finite length of the DBL changes cardinally the shape of impedance spectrum [287,288]. The finite length of the DBL at membrane surface was accounted by Sistat and Pourcelly [289] and Larchet et al. [46], when basing on the Nernst-Planck equation and the LEN assumption. The value of the DBL thickness, δ, is considered as a fitting parameter. In Reference [46], the time dependence of δ is implicitly taken into account in conditions of current-induced gravitational convection. Recently, van Soestbergen et al. [290] presented a theory based on the Nernst-Planck-Poisson equations. The value of δ was fixed as the distance between two parallel electrodes forming the cell; electron-transfer reactions across the Stern monolayer at the electrode surface were taken into account [290]. Mareev et al. [269] took into account the decrease of δ, occurring with time, due to the development of current-induced convection. The value of δ is considered as a function of the Donnan PD (∆φ Don ) at the membrane/depleted solution interface. When an overlimiting current density is applied and the time approaches the transition time, current-induced convection (namely, electroconvection) arises near the depleted interface and causes a decrease of δ. The δ(∆φ Don ) function is assumed a linear one containing two adjustable parameters: a threshold value of ∆φ Don , which relates to the onset of electroconvection and a value of ∆φ Don related to the steady state.
As it is shown in Reference [269], the model describes well experimental chronopotentiograms of a CMX homogeneous CEM at all applied overlimiting current densities. The comparison of the calculation using the Sand theory and the model [269] with experimental transition times found for different current densities is shown in Figure 14. As it can be seen, the Sand transition time is less than the experimental value; this difference increases with decreasing (overlimiting) current density. The explanation is that the assumption of the semi-infinite DBL does not allow accounting for the contribution of forced (or natural) convection. The electrolyte delivery by convection slows down the growth of PD with time. However, when the current density is sufficiently high (>1.5 i lim in Figure 14), the concentration profile develops so rapidly that the DBL thickness does not reach the values, where the contribution of forced convection is significant; therefore, at i > 1.5 i lim , the model of Mareev et al. [269] and the Sand theory give the same results and both agree with the experiment. The same conclusion concerning the agreement between the Sand theory and a model taking into account the finite length of δ was obtained by van Soestbergen et al. [290] in the case of electrode systems.   (76); τModel refers to the finite-length timedependent diffusion layer model [269]; the DBL thickness and the limiting current density calculated using the Lévêque Equations (60) and (61) are δLev = 260 μm, ilim = 2.0 mA/cm 2 , respectively. Reproduced with permission from [269], Elsevier, 2016.

The Use of Electric Current Stream Function. Comparison with the Sand and Choi and Moon Theories
Electrochemical responses of materials having heterogeneous surface mainly depend on the parameters of surface heterogeneity [15,274,291]: on the size and shape [292,293], as well as the surface fraction [274,292,[294][295][296] of conductive and nonconductive areas. Most of mathematical descriptions of ionic transport in electrochemical systems are restricted to 1D geometry [296], not relevant to a heterogeneous permselective surface.
The passage from 1D to 2D and especially 3D geometry sharply increases mathematical difficulties. Nevertheless, there are a large number of systems with permselective surface for which multidimensional effects are important: mass transfer at partially blocked electrodes [292,295], microelectrodes [297,298], proton-exchange membrane fuel cells and batteries [299]; flow electrochemical devices, such as electrolysers and electrodialysers [80,300], electrochemically heterogeneous graphene [293] or carbon nanofilaments [274], heterogeneous IEMs [15,274,291] and so forth. Development of approaches for handling multidimensional modelling of interface phenomena represents an important issue.
2D/3D mathematical description of the PD response of membrane/electrode systems is complicated since an integral condition of current intensity is needed at the permselective surface. This constitutes a problem as the current density distribution over the conductive area is not known a priori. Rubinstein et al. [274], Green and Yossifon [301,302] used the condition of uniform current density distribution on the conductive regions of membrane surface. However, this condition was not verified using another approach. The use of the electric current stream function [303] seems promising for the solution of this kind of problems, since this approach does not require any assumption about the character of current density distribution. Stream functions are widely used in fluid mechanics [304]. The main feature of this approach is that it allows calculation of the fluid velocity distribution when the volumetric flow rate is known in heterogeneous systems [305]. Another advantage is that any solution obtained with the use of stream function satisfies automatically the law of mass conservation. Electrochemical responses of materials having heterogeneous surface mainly depend on the parameters of surface heterogeneity [15,274,291]: on the size and shape [292,293], as well as the surface fraction [274,292,[294][295][296] of conductive and nonconductive areas. Most of mathematical descriptions of ionic transport in electrochemical systems are restricted to 1D geometry [296], not relevant to a heterogeneous permselective surface.
The passage from 1D to 2D and especially 3D geometry sharply increases mathematical difficulties. Nevertheless, there are a large number of systems with permselective surface for which multidimensional effects are important: mass transfer at partially blocked electrodes [292,295], microelectrodes [297,298], proton-exchange membrane fuel cells and batteries [299]; flow electrochemical devices, such as electrolysers and electrodialysers [80,300], electrochemically heterogeneous graphene [293] or carbon nanofilaments [274], heterogeneous IEMs [15,274,291] and so forth. Development of approaches for handling multidimensional modelling of interface phenomena represents an important issue.
2D/3D mathematical description of the PD response of membrane/electrode systems is complicated since an integral condition of current intensity is needed at the permselective surface. This constitutes a problem as the current density distribution over the conductive area is not known a priori. Rubinstein et al. [274], Green and Yossifon [301,302] used the condition of uniform current density distribution on the conductive regions of membrane surface. However, this condition was not verified using another approach. The use of the electric current stream function [303] seems promising for the solution of this kind of problems, since this approach does not require any assumption about the character of current density distribution. Stream functions are widely used in fluid mechanics [304]. The main feature of this approach is that it allows calculation of the fluid velocity distribution when the volumetric flow rate is known in heterogeneous systems [305]. Another advantage is that any solution obtained with the use of stream function satisfies automatically the law of mass conservation.
Apparently, G. Taylor et al. [306] were the first who has introduced the stream function for modelling electric current across heterogeneous surface. They used the mathematical analogy between the flow of fluid and that of electric current. Later on, this approach was applied for simulation of electric current density distribution in semiconductors [307], in superconductors [308], in thin circuits [309], on the earth's surface [310], in electromagnetic devices [311] and other systems. In membrane systems, this method was applied by Pismensky et al. [312] and Mareev et al. [303]. 3D modelling in Ref. [303] was used for describing the chronopotentiograms of a heterogeneous IEM; the results of computations were compared with experiment and with other known theories, namely, that by Sand [284], Choi and Moon [291] and Rubinstein et al. [274].
The geometry used in Ref. [303] is shown in Figure 15. A hexagonal distribution of conductive areas (the open circles in Figure 15) over the membrane surface is assumed. Apparently, G. Taylor et al. [306] were the first who has introduced the stream function for modelling electric current across heterogeneous surface. They used the mathematical analogy between the flow of fluid and that of electric current. Later on, this approach was applied for simulation of electric current density distribution in semiconductors [307], in superconductors [308], in thin circuits [309], on the earth's surface [310], in electromagnetic devices [311] and other systems. In membrane systems, this method was applied by Pismensky et al. [312] and Mareev et al. [303]. 3D modelling in Ref. [303] was used for describing the chronopotentiograms of a heterogeneous IEM; the results of computations were compared with experiment and with other known theories, namely, that by Sand [284], Choi and Moon [291] and Rubinstein et al. [274].
The geometry used in Ref. [303] is shown in Figure 15. A hexagonal distribution of conductive areas (the open circles in Figure 15) over the membrane surface is assumed.
The three-layer system 'membrane + depleted DBL + enriched DBL' adjacent to the membrane from its both sides is considered as a system of cylindrical cells (Figure 15). The membrane upper surface layer consists of a conductive disk of radius R1 within a nonconductive ring of external radius R2. The DBL thickness, δ, is assumed identical from both membrane sides. Using cylindrical symmetry, the system under study may be presented in two coordinates r and z. The governing equations are the same as in the case of ED channel with heterogeneous membrane described in Section 5.3: the Nernst-Planck transport equation (without convective term), the continuity equation for ion transfer and the LEN condition. However, the time-dependent problem is considered here, that is, . As in Section 5.3, the flux continuity condition at the conductive surface/solution boundaries is applied. At the non-conductive boundaries, the zero flux condition like Equation (65) No condition on the distribution of normal current density over the conductive regions is imposed. Instead, the electric current stream function [306,312] η is introduced as follows: Figure 15. Distribution of round conductive areas over the membrane surface in hexagonal geometry (a). The DBL of thickness δ is assumed independent of lateral position (b). Conversion of the 3D system (b) into a 2D simulation problem (c). Reproduced with permission from [303], American Chemical Society, 2016.
The three-layer system 'membrane + depleted DBL + enriched DBL' adjacent to the membrane from its both sides is considered as a system of cylindrical cells ( Figure 15). The membrane upper surface layer consists of a conductive disk of radius R 1 within a nonconductive ring of external radius R 2 . The DBL thickness, δ, is assumed identical from both membrane sides. Using cylindrical symmetry, the system under study may be presented in two coordinates r and z.
The governing equations are the same as in the case of ED channel with heterogeneous membrane described in Section 5.3: the Nernst-Planck transport equation (without convective term), the continuity equation for ion transfer and the LEN condition. However, the time-dependent problem is considered here, that is, ∂c i /∂t = −div → j i = 0. As in Section 5.3, the flux continuity condition at the conductive surface/solution boundaries is applied. At the non-conductive boundaries, the zero flux condition like Equation (65) is used. In galvanostatic regime, the average current density across the membrane surface, i av , is a known parameter. It is linked with the unknown local current density as [303]: No condition on the distribution of normal current density over the conductive regions is imposed. Instead, the electric current stream function [306,312] η is introduced as follows: The charge conservation (continuity) equation in the cylindrical geometry is written in the form: Note that Equation (79) holds, if the rate of the space charge density (ρ e ) variation with time is negligible.
Equation (78) assures that Equation (79) is satisfied at any function η, if only the condition for equality of mixed derivatives of η is met.
The condition of uniform current density distribution over the conductive areas according to Refs. [274,313] is expressed as: However, as it is shown in Ref. [303], when integral condition, Equation (77), is applied, the current lines are concentrated within conductive areas at the boundary with the neighbour non-conductive region. The application of the uniform current density distribution, Equation (80), results in too fast decrease with time of the interfacial electrolyte concentration in the centre of conductive areas, which leads to a shorter transition time in comparison with the model [303].
The model [303] formulated above allows computation of chronopotentiograms and the transition times. As the LEN assumption and a fixed DBL thickness are used in this model, the transition time corresponds to the state where the PD tends to infinity, as in the model [269] described in Section 6.1, where the current-induced transport is implicitly taken into account.
It is known, that the transition time for heterogeneous membranes is lower than that for homogeneous ones [15,291]. However, the only mathematical description of the transition time for an electrically non-homogeneous permselective interface was made by Choi and Moon [291]. They had modified the Sand equation, Equation (76), replacing the average current density, i av , by the local current density, i local = i av /γ, across the conductive areas: The relationship between the transition times expressed according to Sand and to Choi and Moon [291] is expressed as τ Choi−Moon = τ Sand γ 2 .
The transition time calculated using the 3D model [303] is in a good agreement with the experimental data (Figure 16a) obtained for a specially designed membrane. The substrate is a polyethylene terephthalate track-etched membrane with cylindrical randomly distributed pores with a radius of 12.93 ± 0.22 µm. The substrate was coated on one side with a Nafion ® perfluorinated resin solution of thickness d 2 about 3 µm, the Nafion ® solution filled also the pores. Thus, a membrane with a smooth homogeneous permselective surface on one side and a heterogeneous surface on the other was obtained. The circular conductive areas filled with the Nafion ® material occupy about ≈7.7% of the surface, 92.3% of the surface was non-conductive. concentration distribution far from the heterogeneous surface is the same as in the case of homogeneous surface. The Choi-Moon theory [291] assumes the current lines normal to the conductive surface, hence, the funnel effect [274] is not taken into account and tangential transport is considered negligible. This condition is satisfied, when R1/δ and i are both sufficiently great ( Figure  16b,c). No fitting parameters in the model of Mareev et al. [303] are used. The curve computed for the lower bound of the range of possible radii, R 1 = 12.71 µm, is below the experimental points; that for the upper bound, R 1 = 13.15 µm, is above. The transition times calculated using the Sand equation, Equation (76), are essentially higher than the experimental ones and that calculated using the Choi-Moon equation are much lower (Figure 16a). The calculated values of τ using the 3D model [303] but under condition (80) are lower too. In fact, Equations (76) and (81) give two extreme asymptotes of the model [303]. As we mention in the previous section, the Sand Equation (76) holds, when the surface is homogeneous and the current density i > (1.5 ÷ 2)i av lim [269] (Figure 16b,c). The Sand equation holds also when γ < 1 and the ratio R 1 /δ is sufficiently small (Figure 16c). In this case, the concentration distribution far from the heterogeneous surface is the same as in the case of homogeneous surface. The Choi-Moon theory [291] assumes the current lines normal to the conductive surface, hence, the funnel effect [274] is not taken into account and tangential transport is considered negligible. This condition is satisfied, when R 1 /δ and i are both sufficiently great (Figure 16b,c).

Phenomenon of Two Transition Times of Chronopotentiograms for Membranes with Heterogeneous Surface
When an electric current is applied to an IEM with electrically heterogeneous surface, the current lines are concentrated within the conductive regions. With the passage of time, the electrolyte concentration decreases first near the conductive surface regions, then the concentration changes propagate to the solution domains near the non-conducting regions ( Figure 17). Tangential diffusion, which enhances the delivery of electrolyte to the surface of conductive regions occurs in this system.  (circles) are compared with the simulation using the 3D model [303] with integral current condition (77) (curve 1) or condition [274] (80) (curve 2) and the Choi-Moon theory [291], Equation (81)

Phenomenon of Two Transition Times of Chronopotentiograms for Membranes with Heterogeneous Surface
When an electric current is applied to an IEM with electrically heterogeneous surface, the current lines are concentrated within the conductive regions. With the passage of time, the electrolyte concentration decreases first near the conductive surface regions, then the concentration changes propagate to the solution domains near the non-conducting regions ( Figure 17). Tangential diffusion, which enhances the delivery of electrolyte to the surface of conductive regions occurs in this system. Figure 17. Distribution of electrolyte concentration in the DBL at a heterogeneous membrane surface, calculated using the 3D model [303] at i = 1.25ilim in moments of time 0.2, 0.8, 1.4 … 5.6 s. The parameters used in simulation relate to a membrane with randomly distributed conductive areas with the radius of 12.95 μm, the fraction of the conductive surface area γ = 0.077, the DBL thickness δ = 260 μm. The centre of the conductive area refers to r = 0, the boundary between the conductive and nonconductive regions is shown with the dashed line. Reproduced with permission from [314], Elsevier, 2018.
The contribution of tangential diffusion is detected by appearance of two transition times in chronopotentiometry, as it was observed by Butylskii et al. [314].
The first transition time, τ1, is attained, when the electrolyte concentration at the conductive regions of the surface becomes sufficiently low to give rise current-induced convection developed locally at the boundaries between the conductive and non-conductive regions. The second transition time, τ2, is reached when a critically small concentration is attained at the whole membrane surface, including the conductive and non-conductive regions [314]. As the 3D model [303] predicts, when there is no current-induced convection, the electrolyte concentration remains rather elevated at the non-conductive regions, while it reaches very small values at the conductive regions ( Figure 17). According to this model, only one transition time τ1 < Sand τ can occur in a heterogeneous membrane system if the ions in the depleted DBL are transferred by electrodiffusion solely. The decrease of concentration near the non-conductive regions may be caused by the current-induced convection, apparently electroconvection occurring as electroosmosis of the first kind [314]. Electroconvection The contribution of tangential diffusion is detected by appearance of two transition times in chronopotentiometry, as it was observed by Butylskii et al. [314].
The first transition time, τ 1 , is attained, when the electrolyte concentration at the conductive regions of the surface becomes sufficiently low to give rise current-induced convection developed locally at the boundaries between the conductive and non-conductive regions. The second transition time, τ 2 , is reached when a critically small concentration is attained at the whole membrane surface, including the conductive and non-conductive regions [314]. As the 3D model [303] predicts, when there is no current-induced convection, the electrolyte concentration remains rather elevated at the non-conductive regions, while it reaches very small values at the conductive regions ( Figure 17). According to this model, only one transition time τ 1 < τ Sand can occur in a heterogeneous membrane system if the ions in the depleted DBL are transferred by electrodiffusion solely. The decrease of concentration near the non-conductive regions may be caused by the current-induced convection, apparently electroconvection occurring as electroosmosis of the first kind [314]. Electroconvection enhances the exchange between the solution domains, which are located near the conductive and non-conductive surface regions. The onset of this enhancement (occurring at t = τ 1 < τ Sand ) reduces the rate of the PD growth, which gives the first inflection point on the chronopotentiogram. The appearance of larger electroconvective vortices, when the concentration reaches a critically small value on the entire membrane surface (t = τ 2 ≈ τ Sand ), gives the second inflection point. τ 2 is close to τ Sand , since the latter event occurs also at homogeneous membrane surface. However, electroconvection occurring when the near-surface concentration is small everywhere can enhance mass transfer and slows down the rate of concentration decrease resulting in increasing τ 2 higher than τ Sand . This case is shown in Figure 18 for the M1 and M2 membranes. , since the latter event occurs also at homogeneous membrane surface. However, electroconvection occurring when the near-surface concentration is small everywhere can enhance mass transfer and slows down the rate of concentration decrease resulting in increasing τ2 higher than Sand τ . This case is shown in Figure 18 for the M1 and M2 membranes.  [303]: curve 1 relates to a hypothetical homogeneous membrane, curves 2 and 3, to the M1 and M2 heterogeneous membranes, respectively. Reproduced with permission from [314], Elsevier, 2018.
The experimental values of τ2 are essentially higher than the results of calculations for these membranes using the 3D model [303], which takes into account the electrical heterogeneity of the surface but not account for electroconvection. Note that the increase in the values of τ2 with decreasing the lim / i i ratio is due to the effect of the finite length of DBL described in Section 6. lim / i i curve passes lower than the similar curves for the M2 membrane, which is apparently due to a low fraction of the conductive surface area.

Conclusions
In conclusion, we once again emphasize that the heterogeneity of both membrane bulk and surface plays an important role in conditioning membrane properties and enhancing their performance in electroseparation and energy production applications. The recent research shows that design of the membrane bulk affects such membrane properties as conductivity, diffusion permeability and permselectivity (mainly as related to the selective transport of ions of a certain sign of charge). An important example is the immobilization of inorganic nanoparticles in the membrane bulk, which increases membrane conductivity and permselectivity, especially at low relative humidity. Modification of membrane surface exerts a strong impact on the properties related to mass exchange between the solution and membrane surface. The use of membranes with tailored surface τ using the 3D model [303]: curve 1 relates to a hypothetical homogeneous membrane, curves 2 and 3, to the M1 and M2 heterogeneous membranes, respectively. Reproduced with permission from [314], Elsevier, 2018.
The experimental values of τ 2 are essentially higher than the results of calculations for these membranes using the 3D model [303], which takes into account the electrical heterogeneity of the surface but not account for electroconvection. Note that the increase in the values of τ 2 with decreasing the i/i lim ratio is due to the effect of the finite length of DBL described in Section 6.1. The behaviour of the MA-41 membrane is similar: at i/i lim ≤ 1.25, there are two transition times. However, the τ/τ Sand vs. i/i lim curve passes lower than the similar curves for the M2 membrane, which is apparently due to a low fraction of the conductive surface area.

Conclusions
In conclusion, we once again emphasize that the heterogeneity of both membrane bulk and surface plays an important role in conditioning membrane properties and enhancing their performance in electroseparation and energy production applications. The recent research shows that design of the membrane bulk affects such membrane properties as conductivity, diffusion permeability and permselectivity (mainly as related to the selective transport of ions of a certain sign of charge). An important example is the immobilization of inorganic nanoparticles in the membrane bulk, which increases membrane conductivity and permselectivity, especially at low relative humidity. Modification of membrane surface exerts a strong impact on the properties related to mass exchange between the solution and membrane surface. The use of membranes with tailored surface heterogeneity allows a significant increase of the limiting current density and overlimiting ion transport rate. It possible to talk about the formation of surfaces, self-inducing the mixing of the adjacent solution. Mathematical modelling brings a new insight into complicated transport phenomena occurring in heterogeneous systems. This deeper understanding may be very useful for the development of high-performance membranes and membrane devices. In particular, the modelling suggests that the very important factors affecting the intensity of electroconvection mixing are the size and surface fraction of the conductive/non-conductive surface regions. Numerical computations can give optimum values of these parameters. However, there is no good accordance between the computations made in different laboratories. Besides, the role of charge density and degree of hydrophilicity of surface is not yet well understood, while these parameters are also important.
Although there is a strong progress in understanding of the response of membrane performance on certain modifications of their bulk and surface, some important knowledge gaps remain. Most experiments and theoretical calculations are made for NaCl solutions and in conditions proper for laboratory tests. There are no studies on the effects of these modifications on the performance of real ED, no estimations of how these effects depend on solution concentration and forced flow velocity, what is the influence of the cell geometry, intermembrane distance, geometric shape of surface patterns, the presence of intermembrane spacers.
Even more questions remain as to the effect of modification on the membrane antifouling behaviour. The understanding of the relationships between the surface properties and antifouling resistance is rather poor. No mathematical modelling of these effects is available. The problem here is complicated by complexity of coupling between mass transport and chemical reactions, especially when the rate of these reactions is quite low. A similar situation occurs with respect to the relationship between membrane properties and the characteristics of competing transfer of ions having the same charge sign but different valence. However, given the importance of these problems for the sustainable development of society, we can expect rapid progress in this area. . The authors thank the Core Facility "Environmental Analytical Center" at the KubSU(unique identifier RFMEFI59317X0008) and and "Interfaces, Physical Chemistry and Polymers" at the IEM, University of Montpellier, for providing their equipment. V.N. thanks IEM, Montpellier, for giving possibility to realise this study.