Transitions and Instabilities in Imperfect Ion-Selective Membranes

Numerical investigation of the underlimiting, limiting, and overlimiting current modes and their transitions in imperfect ion-selective membranes with fluid flow through permitted through the membrane is presented. The system is treated as a three layer composite system of electrolyte-porous membrane-electrolyte where the Nernst–Planck–Poisson–Stokes system of equations is used in the electrolyte, and the Darcy–Brinkman approach is employed in the nanoporous membrane. In order to resolve thin Debye and Darcy layers, quasi-spectral methods are applied using Chebyshev polynomials for their accumulation of zeros and, hence, best resolution in the layers. The boundary between underlimiting and overlimiting current regimes is subject of linear stability analysis, where the transition to overlimiting current is assumed due to the electrokinetic instability of the one-dimensional quiescent state. However, the well-developed overlimiting current is inherently a problem of nonlinear stability and is subject of the direct numerical simulation of the full system of equations. Both high and low fixed charge density membranes (low- and high concentration electrolyte solutions), acting respectively as (nearly) perfect or imperfect membranes, are considered. The perfect membrane is adequately described by a one-layer model while the imperfect membrane has a more sophisticated response. In particular, the direct transition from underlimiting to overlimiting currents, bypassing the limiting currents, is found to be possible for imperfect membranes (high-concentration electrolyte). The transition to the overlimiting currents for the low-concentration electrolyte solutions is monotonic, while for the high-concentration solutions it is oscillatory. Despite the fact that velocities in the porous membrane are much smaller than in the electrolyte region, it is further demonstrated that they can dramatically influence the nature and transition to the overlimiting regimes. A map of the bifurcations, transitions, and regimes is constructed in coordinates of the fixed membrane charge and the Darcy number.


Introduction
Problems of electrokinetics and micro-and nanofluidics have recently attracted a great deal of attention due to rapid developments in micro-, nano-, and biotechnology. Among the numerous modern micro-and nanofluidic applications of electrokinetics are micropumps, micromixers, currents [22]. Hence, the VC-curve can consist of two segments instead of three: the underlimiting and the overlimiting currents.
To date, there has not been a single study that takes into account the effect of nonzero hydrodynamic velocities within a porous membrane on the electrohydrodynamic response without. The reason for this is velocities in the membrane are considerably smaller than the velocities in the adjacent fluid; hence, the hydrodynamics inside the porous membrane seems reasonably negligible (see monograph by Zabolotsky and Nikonenko [27]). Moreover, such accounting greatly complicates the formulation and solution of the problem.
We point out that in addition to ion-selective membrane systems, composite systems consisting of a porous medium and an adjacent fluid layer are used in various engineering applications: drying processes, solid-matrix heat exchangers, electronics cooling, thermal insulation, heat pipes, nuclear reactors, and porous journal bearings. Mathematical models of such systems are well-developed and applied to numerous flows. An important for us result of these investigations is that even small velocities in the porous medium can have a significant effect on the loss of stability of the whole system, see for example [28]. Hence, it stands to reason to expect the same kind of influence in the membrane composite system; this fact inspired the present work.
For correct analysis of any composite flow, imposition of appropriate boundary conditions at the interface liquid-porous medium is crucial. For that reason, many investigators have proposed different types of interfacial conditions between the porous medium and the adjacent fluid layer, as summarized and compared in the work of Alazmi and Vafai [29]. There are two commonly recognized approaches to composite systems. The first one to match the momentum equations in the fluid-porous medium inter-region is to make use of an empirical expression proposed by Beavers and Joseph [30]. This expression relates the rate of strain in the fluid phase to the difference between the tangential velocities in the fluid and the porous media adjacent to the inter region. These boundary conditions require determination of an empirical parameter.
In the second approach, Ochoa-Tapia and Whitaker [31,32] presented a mathematical model for the jump of the shear stresses that involves Brinkman's correction [33] to Darcy's law in the porous layer. The method of volume averaging is used to derive a stress jump boundary condition. In this case, the differential equations that govern the momentum transport in both the fluid and porous layers are rendered in the same order so that it is possible to match the rates of strain in both regions. The stress-jump condition has an inherent problem; it involves an unknown coefficient β, which needs to be fitted either experimentally or numerically; it is shown to depend on the membrane texture: the porosity, the Darcy number and the pore characteristic size. Comparison of these two approaches (see, for example [34]) shows that they gives qualitatively rather similar results. The objective of the present work is to apply the Ochoa-Tapia and Whitaker [31,32] approach to the ion-selective membrane problem, following Kuznetsov [35] and take β as an independent parameter. This formulation and details on the numerical procedure are presented in the next section, followed by a discussion of the main results and conclusions.

Dimensional Equations
Let us consider a three layer system occupying the space −h <ỹ < 2h, in which there is fluid in the regions −h <ỹ < 0 andh <ỹ < 2h overlying a porous cation exchange membrane region 0 <ỹ <h saturated with the same fluid (see Figure 1). Notations with tilde are used for the dimensional variables, as opposed to their dimensionless counterparts without tilde. (x,ỹ) are the Cartesian coordinates, wherex is directed along the membrane interface andỹ is normal to it. The fluid is assumed to be a symmetric (valence or charge number, z + = −z − = 1), binary electrolyte with an equal diffusivity of cations and anionsD, dynamic viscosityμ, and electric permittivityε. Anode Cathode + -Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange Imperfect cation-exchange ∂Ũ ∂x wherec ± , (Ũ,Ṽ),Π andΦ are molar concentrations of ions, velocity, pressure and electric potential respectively;F is the Faraday number,R is the universal gas constant, andT is the absolute temperature. Let us consider the porous membrane flanked into the electrolyte fluid. For the sake of simplicity, the porosity of the membrane is taken equal to unity, the viscosity, permittivity and diffusivity of the porous medium are assumed the same as their analogues in the fluid layer. Then the ion transport and the Poisson equations in the porous layer 0 <ỹ <h can be written as, wherec ± , (Ũ,Ṽ),Π andΦ are molar concentrations of ions, velocity, pressure, and electric potential respectively, andÑ is the fixed charge density. As the membrane is defined to be a cation exchange membrane,Ñ is positive. The governing equations for the hydrodynamics are those of the Darcy-Brinkman model in the Stokes approximation, corrected by the Coulomb forces in the right hand side of the momentum equations, ∂Ũ ∂x whereκ is a permeability of the porous material. For the sake of simplicity the porosity of the membrane is taken equal to 1. This system must be complemented by the proper boundary conditions (BCs). With respect to the hydrodynamics boundary conditions, we assume the continuity of the velocity components and normal stresses, at the interfacesỹ =h andỹ = 0. According to the Darcy-Brinkman model, which is proposed by Ochoa-Tapia and Whitaker [31,32], the tangential stresses have a jump, crossingỹ =h andỹ = 0, µ ∂Ũ ∂ỹ Here β is an empirical constant which depends on the porosity, the Darcy number and the characteristic pore size. In the present work we adopt Kuznetsov's approach [35] and assume β is an independent parameter, −1 ≤ β ≤ 1.
Let us consider the BCs for the electrostatics and the ion concentration. The electric potential, the concentration of cations and anions and their the fluxes, the electric field normal to the interfaces y = 0 andỹ =h, are all continuous at the interfacesỹ = 0 andỹ =h. Under our assumptions these conditions turn into the following, At the outer boundary of the diffusion layers,ỹ = −h andỹ = 2h, the so called reservoir conditions [11] are applied: where ∆Ṽ is a fixed potential drop.
An important characteristic of the electrodialysis cell is the electric current density through the interfaceỹ =h, which is determined by the flux of cations and anions, Adding initial conditions for the cations and anions completes the system (1)- (19). These initial conditions arise from the following viewpoint: when there is no potential difference ∆Ṽ, the distribution of ions is homogeneous and neutral, The neutral conditions are corrupted by a small environmental random noise, so called "room disturbances".

Dimensionless Equations
To render the Equations (1)- (20)  Here,F is the Faraday constant,R is the universal gas constant, andT is the absolute temperature. The problem is described by the following non-dimensional equations, Here γ = 1 in the porous membrane, 0 < y < 1, and γ = 0 in the fluid region, −1 < y < 0 and 1 < y < 2.
Upon introducing of the stream function Ψ, U = ∂Ψ/∂y and V = −∂Ψ/∂x and excluding the pressure, the Stokes equations turn into the one biharmonic equation with the Darcy and Coulomb corrections, with the BCs at the upper boundary of the electrodialysis system, y = 2: the BCs on the lower boundary, y = −1: the BCs at the electrolyte-membrane interface, y = 1, Similar with (27)- (30), BCs must be set at the boundary y = 0. The fluxes of cations and anions at the interface y = 1 can be obtained from the following relations, giving the density of the electric current as the difference between these two fluxes as follows, Adding initial conditions for the cations and anions (20) in the dimensionless form completes the system, t = 0 : c ± = 1 + "room disturbances".
The problem is described by the following six dimensionless parameters: whereλ D = εΦ 0 /c ∞F is the Debye length, the coupling coefficient κ depends only on the physical properties of the electrolyte. Just to provide perspective, typical values of dimensional and dimensionless parameters are given below. The typical bulk concentration varies within the window,c ∞ = 1 − 10 3 mol/m 3 and the typical potential drop is about ∆Ṽ = 0 − 5V. The diffusivity is takenD = 10 −9 m 2 /s. κ is fixed for a given electrolyte but typically lies between 0.05 to 0.5. In this work it is taken, κ = 0.1. The fixed charge concentrationÑ inside the membrane depends on the type of membrane used and is usually written in the membrane passport; for regular membranes used for desalination, this value is of the order ofÑ = 10 3 to 2 × 10 3 mol/m 3 (see Ref. [27].) The membrane thicknessh for regular desalination devices is on the order of 0.5 mm, however, for some other devices in use, the thickness can be much smaller, micrometers and even down to 10 nm, see Ref. [21]. The Debye lengthλ D , depending on the concentrationc ∞ , lies within the window 1 to 100 nm; accepted regular typical value is 10 nm. The permeabilityκ varies withinκ = 10 −18 m 2 toκ = 10 −12 m 2 , hence, the porous characteristic size changes in the same window as the Debye length. It also means that the problem has two small parameters, the Debye number ν =λ D /h and the Darcy number Da = √κ /h. In this case, the dimensionless parameters change in the following windows: ∆V = 0 − 500, ν = 10 −5 − 10 −3 , Da = 0 − 5 × 10 −2 , N = 0.1 − 20, and β varies from −1 to +1. Results of calculations depend only weakly on the Debye number (see [12,[14][15][16]), and so it will be fixed in our calculations, ν = 5 × 10 −4 . The Brinkman coefficient β, properly speaking, is not an independent parameter, and it depends on the texture of the membrane, namely, on the porosity, the Darcy number, and the characteristic pore size. It must be found either experimentally for a given texture of the membrane or using special numerical calculations. In order to avoid these painstaking procedures, we adopted Kuznetsov's approach [35], where −1 < β < 1 was taken as an independent parameter. Fortunately, similar to the case in Ref. [35], we found that dependence on β is very weak and most calculations were made for β = 0. In subsequent calculations, three dimensionless parameters are varied, namely, ∆V, N and Da, and we explain this weak dependence on β.

Numerical Method and Its Justification
For a wide parameter range the problem (21)-(33) can be solved only numerically. The main difficulties in the numerical solution arise in the direction normal to the electrolyte-membrane surface, the y-direction. In the electrolyte, near the electrolyte-membrane interface, a thin Debye layer, O(ν), is formed. This layer is connected with a small coefficient in the Poisson equation. Moreover, there ia also a small coefficient in the Brinkman-Darcy equation, the Darcy number. Thus another thin layer, O(Da), arises in the membrane, near the interface. For perfect membranes, with a large fixed charge N, an effect of the charge jump upon crossing the interface must also be taken into account. Small layer thickness leads to rapid changing of the unknowns and large spatial derivatives. The expected distributions of the unknowns are shown in Figure 2a: The charge density ρ = c + − c − has a large jump when crossing the electrolyte-membrane interface, caused by large fixed charge N in the membrane. Moreover, for the depleted electrolyte region the singularity is strengthened by the space charge region with a fast variation of ρ. As it is seen from Figure 2b, the salt concentration K = c + + c − jump in the EDLs is nearly vertical and is close to zero in the depleted region; these facts also must be taken into account during choosing of the numerical method and the grid of the discretization. Variation of the electric potential Φ, Figure 2c, is much stronger in the depleted region and it is becoming more pronounced for larger ∆V. The presented behavior gives a hint to select a proper numerical method and proper discretization. The finite-difference methods show poor performance and efficiency for such thin boundary layers and such rapid jumps. Usually in such cases, mesh stretching near singularities causes stiffness of the system, unacceptably small time step, and even possible failure of numerical convergence. The best method to overcome these difficulties is to apply spectral and quasi-spectral methods for the space discretization (see Canuto et al. [36]), recommended for investigation of the laminar-turbulent transitions and even a well-developed turbulence for Navier-Stokes equations (see Spalart [37]). For the perfect electric membranes a quasi-spectral method was successfully applied in Ref. [13]. Chebyshev's orthogonal polynomials with accumulation of zeros and, hence, the best resolution near the boundaries, are chosen as the Galerkin basis. The τ-modification of the spectral method is employed to satisfy the boundary conditions.
Three kinds of problems are considered in this work, (a) one-dimensional (1D) steady solution for the underlimiting and limiting currents, (b) transition to the overlimiting currents, and (c) well-developed overlimiting regimes. Finding the one-dimensional solution and investigating its linear stability (problems (a) and (b)) require only discretization in y-direction. However, the direct numerical simulation to solve problem (c) also needs discretization along x-direction, which is tangential to the electrolyte-membrane surface. It stands to reason to also apply a quasi-spectral method in the x-direction by choosing the Fourier series as basis functions.

The Underlimiting and Limiting Currents
The one-dimensional (1D) steady-state, ∂/∂x = ∂/∂t = 0, solution of the system (21)-(33) describes both underlimiting and limiting current regimes. The tangential electric field, ∂Φ/∂x, for such solutions is absent and the hydrodynamic flow is not involved in the solution. Both regimes are the result of competition between electromigration. The solution does not depend on the coupling coefficient κ and the Darcy number Da; it is function of the potential drop ∆V, the fixed charge N and the Debye number ν.
The system (21)-(33) can be transformed into a system of ordinary differential equations. For the all three layers, Equations (21) can be integrated once. Taking into account boundary conditions (27)-(31), after long and tedious derivations this system was transformed into the following three equations, Here index k = 1, 2, 3 refers to one of three layers, the fixed space charge in the electrolyte layers is absent, N (1) = N (3) = 0, N = N (2) and y (k) m are constants of integration. The boundary conditions (24) and (30) turn into the ones, At the interfaces y = 0 and y = 1 function Φ and its first derivatives dΦ/dy and d 2 Φ/dy 2 are taken continuous; it follows from the original conditions (30), which, in turn, are valid when permittivity of the membraneε m and its diffusivityD m are equal to the permittivityε and diffusivitỹ D of the electrolyte, respectively.
Three different expansions in the Chebyshev series were exploited in the layers. Upon substitution of these series into Equation (34) and into boundary conditions of their continuity and boundary conditions (35), we obtained a system of nonlinear algebraic equations with respect to unknown coefficients of the expansion. Depending on the voltage ∆V from 20 to 150 terms of the expansion were taken, so that the largest order of the nonlinear system was about 500 equations.
Our system of equations, f = ( f 1 , f 2 , . . . , f n ), can be considered as a system of parameterized, nonlinear algebraic equations, where x = (x 1 , x 2 , . . . , x n ) is a vector of unknowns and the Debye number ν, the voltage ∆V, and the fixed charge N are the parameters of the system. Newton's method was used to solve Equation (36). When the solution x for some parameters ν, ∆V was found, this solution was extended to the solution with close parameters, etc. The previous solution was taken as the initial approximation for the search for the next one and eventually to find the whole family. When Φ(y) and its first derivatives were known, the ion concentrations c ± , the charge density ρ = c + − c − , the salt concentration K = c + + c − , and the ion fluxes j ± were readily found. We begin our analysis with an important quantity, the selectivity of the membrane, j + /j = j + /(j + − j − ), which shows deviation from the perfect membrane with j + /j = 1 and j − /j = 0. It was found that selectivity depends weakly on the voltage, but its dependence on the dimensionless fixed charge, N =Ñ/c ∞ , is strong. For different standard membranes used in industry the dimensional fixed charge is always about 1 mol/L. Hence, it is instructive to plot the selectivity versus the bulk concentrationc ∞ instead of N. Such dependence for the Nafion 120 membrane is shown in Figure 3. For low-concentration salt solutions,c ∞ < 10 −2 mol/L, the membrane with a good accuracy can be assumed perfect. Increasing the bulk salt concentration decreases the selectivity significantly; forc ∞ = 10 mol/L it is only 0.5! Note that our numerical model allows to avoid painstaking experiments and, instead, seek selectivity theoretically. Consider detailed results for three different membranes with the fixed dimensionless charge density N, N = 5, N = 0.5 and N = 0.1 thus having characteristic ion-selectivities which are, respectively, close to perfect, intermediate and imperfect. In Figure 4a-c the dependence of the salt concentration K, the charge density ρ and the potential Φ on the spatial coordinate y for N = 5 (perfect membranes or low-concentration electrolyte solutions)and four values of the voltage ∆V are presented. Pay attention to the jumps in salt concentration, charge density, and potential when crossing the electrolyte-membrane interface for all ∆V. This is characteristic of the Donnan enrichment/exclusion of membranes approaching perfectly ion-selective behavior. For the first values of ∆V, 10 and 50, corresponds to the underlimiting current regime, where ρ and K vary with distance practically linearly, excluding narrow O(ν) EDLs near phase interfaces. For the limiting current regime, ∆V = 100 and 300, neither the behavior in the salt enriched region nor the porous membrane changes qualitatively. However, in the salt depleted region, the behavior of all the unknowns changes dramatically. For K a zone of practically zero salt concentration appears near y = 1, which expands as the potential difference ∆V increases. In the ρ(y)-profile, a maximum appears in the space charge which increases and departs from the membrane with increasing ∆V. The potential distribution, Φ(y), along the membrane, shown in Figure 4c, has its own distinctive features. The potential drop inside the nanoporous membrane is almost linear, i.e., obeying Ohm's law. The potential change in the enriched solution region is relatively small due to its good electrical conductivity. In the case of limiting current regimes, a jump in potential both in the membrane and in the zone of an enriched solution can be neglected in comparison with a change of the potential in the zone of a depleted solution. This is due to the fact that in the last zone the electric conductivity tends to zero at K → 0. As a consequence of the reduced selectivity, at −1 < y < 1; the maximum of the charge in the ρ-distribution of the depleted region, 1 < y < 2, does not appear until higher voltages and becomes less profound. The zone of depleted solution in 1 < y < 2 becomes smaller and Φ(y)-distribution becomes smoother; specifically the transition between the membrane and the electrolyte near the interface y = 0 does not have a jump. The steep gradient at y = 1 remains, however, reduced in comparison to the perfect case due to the increased conductivity in the depleted region. In Figure 6a-c a smaller value of the fixed charge, N = 0.1 (high-concentration electrolyte solutions), is considered. The most astonishing difference between large and small N (or low-concentration and high-concentration solutions) lies in absence of the salt depleted region. In turn, the maximum of the space charge, typical for low-concentration solutions, disappears (compare Figures 4a,b and 6a,b). As a result, the gradient of Φ on the coordinate y in the depleted region decreases even further, see Figure 6c. Summarizing, we can say that for the small N (high-concentration electrolyte solutions) the difference between the underlimiting and limiting currents vanishes. Let us consider this fact from the view point of the voltage-current characteristic; the VC-dependence for the different N = 0.1, 0.5 and 5 is shown in Figure 7. We point to a rather long segment of the limiting currents for N = 5 (low-concentration electrolyte solutions) and its practical absence for N = 0.1; there is no difference between the underlimiting ang limiting regimes for the high-concentration electrolytes. This fact was theoretically pointed out in the work by Rubinstain and Zaltzman [22]. We will return to this VC-characteristic in the next chapters.

Transition to the Overlimiting Currents
As was first discovered theoretically [11,12], the main reason for the transition to the overlimiting regimes is a special kind of electrohydrodynanic instability, the electrokinetic instability. In this subchapter, the transition to the overlimiting currents is considered from the viewpoint of classical linear stability theory. Assume that at some ∆V the 1D quiescent steady-state solution, described by Equations (34) and (35), is disturbed by small sinusoidal perturbations, with the wave number k and the growth (decay) factor λ. These perturbations trigger a hydrodynamic flow so that now the velocity componentsÛ andV are nonzero. The subscript 0 is related to the mean 1D solution, the 'ˆ', to perturbations. Upon substitution into the system (21)-(33), then linearizing with respect to perturbations and omitting the subscript 0 in the mean solution, we get the following system, y = 2 : 2 < y < 1 : 1 < y < 0 : −1 < y < 0 : which is an eigenvalue problem for λ for a system of linear ordinary differential equations (ODEs). The coefficients of these equations depend on the solution of 1D problem (34) and (35). If real part of λ for all wave numbers k is negative, the 1D solution is stable; if real part of λ is positive for at least one k, the 1D solution is unstable and it gives the threshold of transition to the overlimiting currents. The system has two small parameters before the highest order derivatives, the Debye number, The length of each layer is 1. Since the Chebyshev polynomials are defined in the interval −1 < z < +1, the length of all three layers was doubled. The resolution of these polynomials increases as we approach the boundaries of the region; therefore, the Chebyshev polynomials ideally take into account the specifics of the problem. Individually, the Chebyshev polynomials do not satisfy the boundary conditions at the cathode and anode and the continuity conditions. Therefore, a τ-version of the Galerkin method was used.
The essence of the τ-version is that the relations (61) were substituted into the Equations (38)-(42), (46)-(53) and (54)-(59), and the condition of orthogonality was applied of the residual of the right-hand side of the equations. The last 22 conditions of orthogonality were excluded from the system and replaced by boundary conditions (37), (43)-(45), (51)-(52) and (60), where the Galerkin polynomials (61) were also substituted. Eventually, the eigenvalue problem for the ODE system was replaced by a generalized algebraic eigenvalue problem of matrices, as it follows, This problem was solved numerically by a standard QR-algorithm.The maximum dimension of the matrix reached 6000; we usually confined them to 2000.
The transition to the overlimiting currents is described by the following parameters: the voltage, ∆V, the fixed charged density, N, the Darcy number, Da, and the Brinkman coefficient, β (we recall, that the Debye number and the coupling coefficient were fixed, ν = 5 × 10 −4 and κ = 0.1). The Brinkman coefficient β is determined by the texture properties of the nanoporous membranes and the size of its pores; this parameter varies from −1 to 1. Kuznetsov [35] showed that the influence of this parameter is insignificant. However, in his formulation the electrodynamical effects are not involved. Following [35], we also varied the parameter β in the above range, and found that the difference between solutions with β = −1, 0 and +1 is less than 1% for all the studied modes. Physically, this effect can be explained as follows: (a) For large values of N, the membrane is close to a perfect and, hence, the instability is dominated by the surface galvanosmotic slip velocity and must be independent of the membrane inner properties (see Ref. [22,23,25]) and, in particular, of β. (b) In the case of small N, the influence of the inner properties of the membrane as a whole becomes important. However, since for small N the electroconvection is largely caused by the volumetric residual charge in the electrolyte layer rather than the electroosmotic slip velocity (see [25]), again, influence of the inner properties of the membrane layer remains negligible. Hypothetically, the parameter β may be important for the equilibrium electrokinetic instability caused by the electroosmotic velocity, but such a regime was not found in the acceptable ranges of the parameters. Taking into account all this arguments, the results with β = 0 are presented below.
The discrete spectrum of eigenvalues λ k consists of both real and complex-conjugate eigenvalues. We will enumerate the eigenvalues according to their real part, Re{λ 1 } > Re{λ 2 } > Re{λ 3 } > . . . > Re{λ n } > . . .. When the first real eigenvalue λ 1 or the first pair of the complex-conjugate eigenvalues Re{λ 1,2 } crosses the imaginary axis of the complex λ-plane, the 1D steady-state solution looses its stability and induces transition to the overlimiting currents. For large N, the membrane approaches the perfect case and the spectrum is real (see Ref. [17]). In this case, the instability is monotonic. However, small N corresponds to the complex-conjugate pairs of eigenvalues and, hence, instability is oscillatory (Figure 8a,b, respectively). So, transition to the overlimiting regimes is realized through one of two scenarios: Figure 8a the monotonic transition for a perfect membrane (low-concentration electrolyte solutions), or Figure 8b the oscillatory one for the imperfect membranes (high-concentration electrolyte solutions). The instability and transition are connected with the largest eigenvalue, λ = λ 1 . For stable ∆V = ∆V s < ∆V * λ(k) decays for all wave numbers; for the unstable ∆V = ∆V u > ∆V * there is a window of unstable wave numbers with λ > 0; the for critical case ∆V = ∆V * there is critical k = k * with λ(k * ) = 0, while all other wave numbers are stable, λ(k) < 0 for k = k * . The critical value ∆V * gives the electrokinetic instability threshold and critical wave number k * and the boundary between the one-dimensional underlimiting or limiting regimes and two-dimensional microvortex regime, see Figure 9a Let us return to the VC-characteristics for the underlimiting and limiting currents from the previous chapter (Figure 7), described by the 1D steady-state quiescent equilibria. The instability violates this state of equilibrium and leads to the ovelimiting regimes. While the one-dimensional equilibria do not depend on the hydromechanics of the process and thus on the Darcy number, the transition point ∆V * does in fact depend on it (see shaded areas in the Figure). The portion of the VC-characteristic for the supercritical voltage, ∆V > ∆V * , can be described only by the direct numerical simulation of the problem in the full nonlinear formulation (21)- (33). ∆V * is a function of the fixed charge N: (a) For sufficiently large N, N from 5 to 10, the membrane becomes practically perfect (low-concentration electrolyte solutions) and the critical voltage ceases to depend on N. (b) For a small N (high-concentration electrolyte solutions) ∆V * increases with decreasing of N. An interesting feature of high-concentration electrolyte solutions is that the flat portion of the VC-characteristic, responsible for the limiting currents, disappears and transition to the overlimiting regimes happens right from the underlimiting regimes, thus bypassing the portion of the current saturation in the VC curve. This fact was first elucidated in the works [22,23] for a simplified solution for small N.
Marginal stability curves depend on the parameters N and Da in a rather sophisticated way. In Figure 10a a typical marginal stability curve is shown for a small but nonzero Da. For a small subcriticality there is only one unstable mode, region II. However, with increasing voltage a thin region III with a second unstable mode appears near the lower branch of the marginal curve (the long-wave mode). The left boundary of this region is determined by the point "a". With increasing Da-number, this point moves towards the nose of the stability curve. At the same time a point "b" on the upper branch of the stability curve appears; it characterizes an additional unstable region III which arises for the large wave numbers (the short-wave mode). As the Darcy number increases, both points move towards each other and merge at a critical voltage ∆V * , see Figure 10b,c. With increasing of N region III disappears, and the instability is governed by only one unstable real eigenvalue.
Dependence of the critical voltage ∆V * on both N and Da is shown in Figure 11. The dependence of ∆V * on Da is not monotonic; it has a minimum: at N = 5 at Da = 5.1 × 10 −2 , at N = 0.5 at Da = 10 −2 and at N = 5 at Da = 8.1 × 10 −3 . Existence of this minimum is connected with the fact of interplay of the short-wave and long-wave modes. The results are summarized on the mode map, Figure 12, where a composite dependence of ∆V * on N and Da is shown. The critical voltage is highlighted by the background color, see to the left of the map. The nonequilibrium instability is characteristic for the large N (perfect membranes and low-concentration electrolyte solutions). The equilibrium instability begins to prevail with decreasing of N (imperfect membranes and high-concentration electrolyte solutions). The critical voltage, ∆V * , is increasing with decreasing fixed charge N. An increase in the Darcy number enhances the growth of critical voltage ∆V * . Instability is monotonic for perfect membranes and this fact is independent of the Darcy number. As the Darcy number increases, the transition region in N between monotonic and oscillatory instabilities narrows.

Direct Numerical Sumulation of Well-Developed Overlimiting Regimes
In order to find numerical solutions of the full nonlinear system (21)-(33) for the developed overlimiting currents, the quasi-spectral discretization in space was used. With respect to the independent variable x directed along the membrane, the unknown concentrations of anions and cations, the electric potential, and the velocity components in all three layers were represented as Fourier series: where k is a basic wave number for a characteristic length of the channel, 2π/k and nk, n = 2, 3, . . . are the overtones. Therefore, for each layer, M 1 one-dimensional unsteady problems for the each harmonic, n = 0, 1, . . . M 1 − 1 arise, The problem in the membrane layer is now formulated as, for zero harmonic, n = 0, and the Darcy-Brinkman equations turn into, Then, for each function F n (t, y), discretization by y was carried out using the Galerkin method with the choice of Chebyshev's polynomials T n (z) as basis functions. The y-coordinate was stretched twice to be mapped to −1 < z < 1, Chebyshev's polynomials are well suited to our problem, since they have an accumulation of zeros and therefore a high resolution near the Debye and Darcy layers. Substituting the finite Chebyshev series into the system (63)-(71) and using the Lanczos τ method (see [36,37]) to satisfy the boundary conditions (24)-(30) leads to a system of coupled ODEs for the unknown Galerkin coefficients c ± mn (t) and two systems of linear algebraic equations with respect to Φ m,n and Ψ m,n for the each layer. To obtain these systems, all nonlinear algebraic operations are executed in physical space, in the collocation points, while derivatives with respect to both spatial variables x and y are calculated in the space of the Galerkin coefficients. Derivatives of the Chebyshev polynomials are calculated by means of the collocation matrix method (see [36]). The connection between the collocation points and the Galerkin coefficients is performed by means of the fast Fourier transform.
A special method is developed to integrate the system in time. The numerical solution of the problem turned out to be very computationally costly. Therefore (a) the parallelization of the problem was applied on the Lomonosov supercomputer of the Moscow State University, (b) a special numerical scheme was used to discretize the problem in time. The second order Adams-Bashforth scheme for nonlinear terms along with the Crank-Nicholson scheme for linear terms were used in our work [13] for the perfect one-layer membrane. In order to significantly accelerate the calculations this method was changed. The system of ODE's with respect to dF j,n /dt was integrated by the special third-order Runge-Kutta semi-implicit scheme adapted from the work [38].
"Room disturbances" as conditions at t = 0 and initial stage of their evolution. At t = 0 electroneutrality conditions were used, c − = c + = 1 for the electrolyte layers and c − = c + − N = 1 in the membrane layer. These electroneutral conditions always contain small random (thermal) perturbations at a minimum. Sometimes they are called "room disturbances". To mimic the initial small-amplitude broad-banded white noise the initial conditions are presented in the form, which is consistent with the decomposition (62). The amplitudes of these harmonicsĉ ± are assumed to be small and their absolute value and phase θ ± m for the each harmonic mk set by a random number generator uniformly distributed in the region (0, 2π). The simulations were carried out for amplitude of the harmonic disturbance varying from 10 −5 to 10 −7 , with most at 10 −6 . However, changing the amplitude did not change the results, but only the saturation time. Two basis wave numbers were taken; the rough calculations were done at k = 1, but for more precise calculations k = 0.5 were taken.
In particular k = 0.5 were taken for the evaluation of critical ∆V * and k * . Only the ion transport equations contain time derivatives and, hence, require the initial conditions.
The initial evolution of the entire spectrum of small-amplitude harmonics is then described by the linear stability theory (see the previous Chapter). The amplitude of each harmonic mk, m = 1, 2, . . ., for the stable wavenumbers is exponentially decaying, while for the unstable wavenumbers it is growing. As it follows from the linear stability theory, evolution of each harmonic takes place independently from each other and is described by the eigenvalue problem (37)-(60). For the critical conditions ∆V = ∆V * there is a critical wave number k * = m * k, for which the amplitude does not change in time. Amplitudes of all other harmonics, m = m * , are exponentially decaying in time. Comparison between the linear stability theory and our direct numerical calculation is presented in Figure 13 for two values of the fixed charge N. According the linear theory ∆V * 1 ≈ 25.8, k * 1 ≈ 3.42 and ∆V * 2 ≈ 24.2, k * 2 ≈ 4.49, while according to the DNS ∆V * 1 ≈ 28.5 k * 1 ≈ 3.40 and ∆V * 2 ≈ 22.5, k * 2 ≈ 4.50. It shows reasonably good match between the linear stability theory and the direct numerical simulation of the full nonlinear system (21)- (33).
Typical evolution of the initial small-amplitude broad-banded white noise for a nonzero supercriticality ∆V − ∆V * is shown in Figure 14. k = 0.5 was taken as a basic wave number, so that the harmonics, involved in the calculation, were (0, k, 2 k, 3 k, . . . , M 1 k) with M 1 = 128. Evolution begins with a random distribution of the amplitudes of these harmonics versus the wave number, t = 0. At t = 0.003 the linearly stable harmonics decayed. Amplitude of the linearly unstable harmonics, in the wave number window from 5 to 20, increased about 200 times. The amplitude of the "surviving" harmonics was still small and they continued to develop according to the linear stability theory, exponentially growing in time. At t = 0.5017 the initial small-amplitude broad-banded white noise was filtered into practically monochromatic disturbances with a dominant wave number k max = 3, corresponded to the maximum growth rate λ max = λ(k max ) according to the linear theory.  Note, that the primary band focuses due to the filtering until the amplitudes of the disturbances are small. With amplitudes growing, the nonlinear effects began to corrupt this linear filtering. Overtone 3k max = 6 and and zero-frequency bands appeared, as seen in the Figure 14 at t = 0.5017. The exponential growth downstream was also saturated by this nonlinear interaction.

The Nonlinear Stage of Evolution
The next stages of evolution are nonlinear processes with eventual saturation of the disturbance amplitude. For a small supercriticality ∆V − ∆V * , the nonlinear saturation leads to steady periodic electroconvective vortices along the membrane surface. With increasing supercriticality the attractor can be described as a structure of periodically oscillating vortices. With a further increase in supercriticality, the behavior eventually becomes chaotic in time and space.
Established solutions are shown in Figures 15-17, where snapshots of the stream lines Ψ(x, y) are presented. All the calculations are for small supercriticalities ∆V − ∆V * , so that the attractor at t → ∞ is regular.
We begin the discussion with the simplest case of the perfect membranes (low-concentration electrolyte solutions) for different Darcy numbers and voltage, see Figure 15a-c. According to the results of the linear stability analysis, the attractor corresponds to a pair of steady vorticies for all the Darcy numbers. For very small Da, Figure 15a, despite of strong instability and hydrodynamic movement in the depleted electrolyte layer, the membrane and enriched electrolyte layers practically do not react to it and velocity field in the regions 1 < y < 2 and 0 < y < 1 is rather weak. With the Darcy number increasing, see Figure 15b,c, the character of the vortex pair in the depleted region does not change much, but even small increasing of Da leads to a rise in the velocity field in the layers 0 < y < 1 and −1 < y < 0, comparable with the one in the depleted region. With increasing of the supercriticality ∆V − ∆V * , the regime of regular microvortices, through several secondary instabilities and transitions, changes to a chaotic regime. These instabilities and transitions qualitatively are the same, as described in Ref. [17] for N = ∞. Eventually note, that the charge field, ρ = c + − c − , near the electrolyte -membrane interface y = 1 acquires a typical spike-like distribution with an angle at the apex of the spike about 120 • , see Ref. [17], Figure 10 (the charge distribution is not presented in this work). The results of the previous chapter on linear stability show the complex hydrodynamic behavior and, in particular, the appearance of several unstable modes and the change of the monotonic nature of the instability to the oscillatory one at decreasing of the fixed charge N, see Figure 11. The established regimes for N = 0.5, for a small supercriticality ∆V − ∆V * and different Darcy numbers are presented in Figure 16. Again, for very small Da the picture is similar with N = 5 and the attractor is a pair of steady vorticies, Figure 16a. With increasing of Da there is a possibility either of oscillatory or monotonic instabilities, see Figure 16b,c. The velocity field in the electrolyte depleted region initiates the hydrodynamic movement in the membrane and the electrolyte enriched region. The latter has the same order of magnitude as in the layer 1 < y < 2.
According to the linear analysis the stability for N = 0.1 is always oscillatory. The nonlinear direct numerical simulation confirms this fact. Its resalts are presented in Figure 17. For all Da the instability is caused by a complex-conjugate pair and is thus oscillatory. Several vortices appear in the depleted region, 1 < y < 2, the complex behavior begins already for rather small supercriticality. With further very small increasing of ∆V the dynamics becomes chaotic.
For sufficiently large Da, the hydrodynamics in all three regions (electrolyte/membrane /electrolyte) begin to be linked and the analysis presented above, which affects only the depleted zone of the electrolyte and the membrane region, while ignoring the enriched zone, ceases to adequately reflect the behavior of the system.

Conclusions
Using a three layer composite electrolyte-nanoporous membrane-electrolyte model system, three characteristic regimes of VC response for imperfect and perfect ion selective membranes were studied numerically: the underlimiting, the limiting and the overlimiting. Mathematically, the Nernst-Planck-Poisson-Stokes system was taken in the electrolyte layers and the Darcy-Brinkman approach was employed in the membrane layer. To avoid numerical difficulties connected with the Debye and Darcy singularities, the quasi-spectral Galerkin method was applied. The transition to the overlimiting currents was assumed to be connected with the electrokinetic instability. The threshold of the transition was studied by the linear stability theory and the well-developed overlimiting currents were a subject of the DNS the full nonlinear system of equations. It was found that for a large fixed charge density of the membrane, it operates as a perfect membrane, while for a small fixed charge density the behavior of the imperfect membrane is much more sophisticated. In particular, the direct transition from the underlimiting regimes to the overlimiting ones, bypassing the limiting currents, was found to be possible for imperfect membranes. The transition to the overlimiting currents for the perfect membranes is monotonic, while for the imperfect ones it is oscillatory. Despite the fact that velocities in the porous membrane are much smaller than in the electrolyte region, they could dramatically influence to the transition to the overlimiting regimes. Increasing Darcy number generally speaking introduces more oscillatory instability into the system and consequently results in more complex hydrodynamics. Moreover, the results of the preceding analysis can be experimentally verified in a rather straightforward manner by observation of the nature of the instability in a micro-nanochannel system. Such a system has already been employed to investigate the effects of imperfect ion transport on the transient regime [39]. The Darcy number in such systems can be controlled to reasonably observe both oscillatory and monotonic regimes.