**1. Introduction**

Ion exchange membranes play an important role for a wide range of engineering applications [1]. These applications can be divided in two main groups: separation processes and energy production and conversion. The main processes gathered in the separation field are electrodialysis (water desalination and salt pre-concentration), diffusion dialysis (acid and base recovery from industrial waste water), bipolar membrane electrodialysis (production of acids and bases), continuous electro-deionization (production of ultrapure water), and capacitive deionization (water desalination and water softening). Energy production and storage include reverse electrodialysis, fuel cells and redox flow batteries [2,3]. The recent development of microscale devices, such as microfuel cells, lab-on-a-chip biomolecule sensing and micropumps, use electro-kinetic phenomena in the presence of ion exchange membranes [4–6]. These membranes are mainly polymeric matrices composed of backbones on which cationic or anionic groups are attached. Owing to the key role of the ion exchange membrane, the role of the functional groups, polymer architecture and chemical procedure of functionalization have been extensively studied and have led to a very wide variety of membranes [1,7]. In the case of microscale devices, a new class of membrane seems to be promising. The structure of these materials is inspired by biological processes, such as enzymatic catalysis or transport regulation of species through membrane proteins [8–11].

Electro-convective instability involving an ion exchange membrane was predicted for the first time by Rubinstein and Zaltzman [12–16]. The authors showed by means of an asymptotic expansion that the deformation of the non-equilibrium electric double layer (EDL) induces a longitudinal force (along the membrane) initiating fluid motion [13,14]. This new kind of instability, confirmed by experimental visualizations [17–23], enhances the ion flux through the membrane system. This phenomenon is all the more important because electro-convective instability occurs at the microscale as it is not bound to the Reynolds number but to the non-linear coupling between the ionic charge density and the electric field.

Three regimes are commonly observed in a current–voltage response when a potential drop (ΔΦ) is applied normal to the membrane system (see Figure 1). For a small potential drop (region (1)), the ionic conductivity is constant (Ohmic regime). The polarization layer constitutes an electroneutral diffusion layer, over which the concentration gradient is constant, and a quasi-equilibrium EDL. The concentration gradient increases with ΔΦ, which means that the width of the diffusion layer decreases. It follows that beyond a first potential threshold ΔΦ\* (region (2)), the conductivity cancels and the current density reaches a plateau (limiting current regime) [24–26]. A new layer, called the extended space charge (ESC) region, appears between the diffusion layer and the quasi-equilibrium EDL [12]. The ESC region and the quasi-equilibrium EDL form the non-equilibrium EDL. In the ESC region, the electro-neutrality is not fulfilled. The ESC width increases with ΔΦ. Beyond a second threshold ΔΦ\*\* (region (3)), the non-equilibrium EDL becomes unstable. A longitudinal electric charge gradient, located in this layer, induces an electric and a pressure force along the membrane surface. The enhanced ionic transport leads to an increase in current density [27–30]. This last regime is called the over-limiting current regime. The diffusion layer is moved away from the membrane and a new layer, called the mixing layer, takes place between the diffusion one and the non-equilibrium EDL. The width of the mixing layer is much higher than the vortex size owing to a plume of small ion concentration in the region of the rising fluid [21,28,31].

**Figure 1.** Schematic drawing of the current density (*I*) vs. potential drop (ΔΦ) over a membrane system. The regions (1), (2) and (3) localize the Ohmic regime, the limiting current regime and the overlimiting current regime respectively.

Direct numerical simulations are necessary to reach a complete understanding of the electro-kinetic instability mechanism. With this aim, the Navier–Stokes equation, coupled with the Poisson–Nernst–Planck systems, has been solved in previous works for the case of one or two membranes immersed in a stationary reservoir. In this context, instability modes, transition to chaotic motion, 3D electro-convective pattern, bifurcation, non-ideal selectivity and chaotic mixing were investigated [28,31–37]. Recently, the influence of surface heterogeneity on the electro-instability has also been investigated by direct simulations and experimental visualizations [38,39]. In particular, the authors [38] have studied how the electric field, bended by an array of charged P2VP microgels coated on Nafion membranes, enhances the ionic transfer.

Few studies have been published in the case of electro-kinetic instabilities in the presence of an imposed inflow. Kwak et al. [18] found a scaling law governing the vortex thickness. The authors also observed that the vortex advection velocity is dependent on the imposed mean velocity. In a 3D approach, experiments and simulations reveal parallel helical vortices and a merging process along the channel [40]. Nikonenko's group [28,41] described the instability mode appearing successively during a linear potential sweep. The authors observed first a periodic mode with a time period about 0.9 s, when ΔΦ is close to the instability threshold. For a higher potential drop, they observed regularly spaced clusters of counter-rotating vortices. They noticed that vortices can emerge in these clusters. This research group also analyzed the width of different zones of the diffusion layer as a function of the current density. In particular, they found an agreement between experiments and simulations in spite of the difference of inlet concentration.

Despite the quality of the published results, the mechanisms controlling the instability modes remain only partially understood. The goal of the current work is to intend to explain some processes like vortex dissociation/association, counter-clockwise vortex birth and also time evolution of the cluster structure by considering the spatial distribution of the electric and pressure force. As in previous publications [28,33,40,41], the fluid flow is visualized by means of streamlines computed at a given time. Therefore, the time evolution of the vortex layer is a snapshot series of these streamlines named 'static' (SSL) in the article. However, considering the unsteady fluid flow, the fully Lagrangian approach provides access to the real hydrodynamics. This is visualized by trajectories called 'dynamic' (DSL) streamlines in the present work. So, in a second step, the close relation and complementarity between these 'dynamic' streamlines and the vortices computed with the SSL are clearly identified and analyzed with the spatial distribution of force and kinetic energy. This gives a new insight into the role of the vortices in the fluid motion and of the relation between the diffusion layer and the real instability layer. This approach is also extended to ion transfer in the periodic mode. Finally, the results are summarized in the conclusion.

#### **2. Theory**

#### *2.1. Model and Governing Equations*

The elementary electrodialytic cell consists of a channel between two parallel ion exchange membranes. The electrolyte solution enters one end of the channel and exits from the other. The transverse electric potential drop, imposed at the two outer sides of the membranes, induces two polarization layers at their inner sides [41]. In the current work (Figure 2), the membrane system consists of one polarization layer close to a cation exchange membrane of length *Lz*. The bulk is located at *Lx* from the membrane. At the inlet (*z* = 0), the ionic concentration is *co*, a constant transverse potential gradient ΔΦ/*Lx* is imposed and a part of the Poiseuille profile is prescribed depending on the distance between the two membranes *Lm*(> *Lx*). At the membrane surface (*x* = 0), the cation concentration is *c* + interf, the anion flux and the electric potential are zero. These two last conditions mean that the membrane is ideally selective and its electrical resistivity is zero as well. At the bulk limit (*x* = *Lx*), the potential and the ion concentration are imposed constant and equal to ΔΦ and *co*, respectively. The fluid velocity is longitudinal and equal to *U*slip. Figure 3 displays a snapshot of cationic profiles and streamlines for two values of ΔΦ in the case of the monovalent electrolyte (Na+, Cl−).

**Figure 2.** Schematic representation of the numerical model: an electric potential drop ΔΦ is imposed over the width of the electrolyte layer (polarization layer) lying between a bulk, containing an electrolyte solution of concentration *co* and a cation exchange membrane at the surface of which the cation concentration *c* + interf and a non-flux of co-ions are imposed. A Poiseuille flow is imposed at the inlet.

**Figure 3.** Current density vs. potential drop. (**A**): mean current density through the domain (*I*); (**B**): current density at the membrane surface (*I*m). The red vertical bars represent the variance of the current density signal.

*Int. J. Mol. Sci.* **2019**, *20*, 2393

The incompressible fluid motion is described by the Stokes equation (Equation (1)) coupled with the continuity equation (Equation (2)). The ionic transport is governed by the Nernst–Planck equation (Equation (3)) and the electric potential is related to the local ionic concentration by the Poisson equation (Equation (4)):

$$0 = -\overrightarrow{\nabla}P + \mu\Delta\overrightarrow{U} - \rho\_f \left(z^+c^+ + z^-c^-\right)\overrightarrow{\nabla}\Phi\tag{1}$$

$$
\overrightarrow{\nabla} \cdot \overrightarrow{\mathcal{U}} = 0 \tag{2}
$$

$$\frac{d}{dt}\mathbf{c}^{\pm} = \frac{\partial}{\partial t}\mathbf{c}^{\pm} + \overrightarrow{\mathbf{U}} \cdot \overrightarrow{\nabla}\mathbf{c}^{\pm} = D^{\pm}\overrightarrow{\nabla} \cdot \left(\overrightarrow{\nabla}\mathbf{c}^{\pm} + \frac{z^{\pm}\mathbf{c}^{\pm}}{k\_B T}\overrightarrow{\nabla}\Phi\right) \tag{3}$$

$$
\varepsilon\_0 - \varepsilon\_0 \varepsilon\_r \Delta \Phi = Q = z^+ c^+ + z^- c^- \tag{4}
$$

Here *<sup>P</sup>* is the pressure, <sup>μ</sup> is the dynamic viscosity, <sup>→</sup> *U* is the fluid velocity, ρ*<sup>f</sup>* is the fluid density, *z*<sup>±</sup> is the charge of the cation and the anion, *c*± is the concentration of the cation and of the anion, Φ is the electric potential, *t* is the time, *D*± is the diffusion coefficient of the cation and of the anion, *Q* is the electric charge density, ε*<sup>o</sup>* and ε*<sup>r</sup>* are the permittivity of the vacuum and the relative permittivity of the solution, respectively, and *kB* and *T* are the Boltzmann constant and the temperature, respectively.

These equations are scaled as in Druzgalski et al. [31]. The characteristic parameters are the transverse length *Lx*, the diffusional time *tdi f f* = *L*<sup>2</sup> *<sup>x</sup>*/*Dre f* , the diffusional velocity *Udi f f* = *Dre f* /*Lx*, the thermodynamic potential Φ*<sup>T</sup>* = *kBT*/*e*, the 'osmotic' pressure *Po* = μ*Dre f* /*L*<sup>2</sup> *<sup>x</sup>* and the bulk concentration *c*o. Here *e* is the elementary charge, *Dre f* is a reference diffusion coefficient whose value is arbitrary. For a symmetric univalent electrolyte, Equations (1)–(4), in their dimensionless form, become:

$$0 = -\overrightarrow{\nabla}P + \Delta\overrightarrow{U} - \frac{Pe}{2\nu^2}(\mathfrak{c}^+ - \mathfrak{c}^-)\overrightarrow{\nabla}\Phi\tag{5}$$

$$
\overrightarrow{\nabla} \cdot \overrightarrow{\mathcal{U}} = 0 \tag{6}
$$

$$\frac{d}{dt}\mathbf{c}^{\pm} = \frac{\partial}{\partial t}\mathbf{c}^{\pm} + \overrightarrow{\mathbf{U}} \cdot \overrightarrow{\nabla}\mathbf{c}^{\pm} = D^{\pm}\overrightarrow{\nabla} \cdot \left(\overrightarrow{\nabla}\mathbf{c}^{\pm} \pm \mathbf{c}^{\pm}\overrightarrow{\nabla}\Phi\right) \tag{7}$$

$$-2\nu^2 \Delta \Phi = Q = \mathfrak{c}^+ - \mathfrak{c}^- \tag{8}$$

Here *Pe* and ν are the Peclet number and the dimensionless Debye length, respectively:

→

$$Pe = \frac{\varepsilon\_o \varepsilon\_r}{\mu D\_{ref}} \left(\frac{k\_B T}{c}\right)^2 \tag{9}$$

$$\nu = \frac{\lambda\_d}{L\_x} \text{ with } \lambda\_d = \left(\frac{\varepsilon\_o \varepsilon\_r k\_B T}{2e^2 c\_o}\right)^{1/2} \tag{10}$$

For simple ions, the diffusion coefficient is roughly 10−<sup>9</sup> m2/s. Hence, this value will be used to define *Dre f* and the value of Pe is imposed equal to 0.5. Therefore, if *Lx* = 0.1 mm and ν = 10−3, then λ*<sup>d</sup>* = 100 nm, *co* = 10−<sup>5</sup> M (*T* = 300 K) and the diffusional time is 10 s. In what follows, the dimensionless variables will be used. Therefore, the value of dimensionless transverse length *Lx* and of the dimensionless bulk concentration *co* are equal to 1.

The boundary conditions used at the fluid/membrane interface are the no-slip condition, and the imposed value of the cationic concentration *c* + interf is assumed to be equal to *co*. At the inlet, the imposed velocity profile has the following expression:

$$
\Delta I\_z = \Theta(\mathcal{U}\_{\bar{z}}) \frac{\chi}{L\_m} \left(\frac{\chi}{L\_m} - 1\right) \tag{11}
$$

In this expression, *Uz* is the longitudinal velocity averaged over the distance between the two membranes and x ranges from 0 to *Lx*. Therefore, the slip velocity imposed at the bulk boundary is *Uslip* = *Uz*(*z* = 0, *x* = *Lx*) = 6*UzLx*/*Lm*(*Lx*/*Lm* − 1). At the outlet (*z* = *Lz*), the equality of the velocity, of the ion concentration and of the electric potential are imposed on the two sides of the boundary [18]. As the Stokes equation is solved with the finite volume method, the pressure is linearly extrapolated at all the domain boundaries.

At the membrane surface, the no-slip boundary condition is usually applied. However, the slip condition can be justified by the hydrophobic property of membranes. The velocity slip is modeled by the Navier equation, which is characterized by the slip length. The value of this parameter depends on the degree of hydrophobicity and heterogeneity. It is admitted that it is lower than 100 nm for a homogeneous surface and may reach a few micrometers for a grooved surface (superhydrophobicity) [42,43]. Experiments on cation exchange membranes have shown that surface modification by hydrophobization enhances the electro-convection [44,45]. The instability threshold is displaced to lower values of ΔΦ, and at a given value of ΔΦ, the current density increases. The consequence is that water splitting and salt precipitation are reduced. Numerical simulations [41,46] confirm the experimental observations. Shelistov et al. [46] have shown that moderate hydrophobicity increases the critical wavelength and the maximal growth rate, and the unsteady mode appears at lower values of ΔΦ as the hydrophobicity increases.

In this study, *Lz* = 6*Lx* and *Lm* = (10/3)*Lx*, *Uz* = 240 (8 mm/s), and ΔΦ varies up to 35 (0.9 V). Two values of <sup>ν</sup> were used (10−<sup>3</sup> and 1/3 <sup>×</sup> <sup>10</sup>−<sup>3</sup> ) and two cases were also considered for *D*± (*D*<sup>+</sup> = *D*<sup>−</sup> = 1 and *D*<sup>+</sup> = *DNa*<sup>+</sup> = 1.33, *D*<sup>−</sup> = *DCl*<sup>−</sup> = 2.06). In the following, the ionic flux and the current density are defined in their dimensionless form:

$$\overrightarrow{J}^{\pm} = c^{\pm}\overrightarrow{U} - D^{\pm} \left( \overrightarrow{\nabla c^{\pm}} \pm c^{\pm} \overrightarrow{\nabla} \Phi \right) \text{and} \ I = \left\langle \left( \overrightarrow{I} \right)\_{\mathbf{x}} \right\rangle\_{\Omega} = \frac{1}{A} \int\_{\Omega} \left( \left( \overrightarrow{J}^{+} \right)\_{\mathbf{x}} - \left( \overrightarrow{J}^{-} \right)\_{\mathbf{x}} \right) dA \tag{12}$$

In this equation, A is the area of the computational domain Ω.

#### *2.2. Numerical Method*

#### 2.2.1. Stokes and Poisson–Nernst–Planck Equations

The same discretization method and structured grid refinement as the previous article [37] were used. The number of nodes is 180 and 384 in the transverse and longitudinal direction, respectively. In the transverse direction, the common ratio of the geometric series is 0.0234. In the longitudinal direction, the refinement threshold is 2.5 Equations (5)–(8) are integrated in time with a time step δ*t* of 10<sup>−</sup>4.

#### 2.2.2. Streamlines

In References [18,28,30–32,40,41] the hydrodynamic streamlines were computed with the velocity field at a given time. In the present work, these streamlines are named static streamlines (SSL) and are computed with Matlab. In Section 3.2, the instability structure is analyzed with the SSL. They are visualized in black in Figure 4 for example. In Section 3.3, the real trajectories of the fluid elements and ions are computed with the numerical method described below. These trajectories will be named dynamic streamlines (DSL) owing to the use of the unsteady fluid and ion transfer. The trajectories are displayed in red in Section 3.3.

**Figure 4.** Color plot of the cationic concentration. (**A**): ΔΦ = 22, (**B**): ΔΦ = 30. The black lines represent the static streamlines.

Two kinds of dynamic streamlines are computed in this work: the hydrodynamic and the ionic one. As regards to the ions, their local velocity is computed by means of the ratio between the flux density and the concentration: <sup>→</sup> *J* ± /*c*±. As the velocity field and the spatial distribution of concentrations and potential are unsteady, the time integration is performed by means of the fourth order Runge Kutta method.

Consider an elementary volume of fluid located at the point (*zo*, *xo*) at time *to*. Four Eulerian values of velocity are needed to compute the Lagrangian one:

$$
\overrightarrow{\boldsymbol{U}}\_{1} = \overrightarrow{\boldsymbol{U}}(\boldsymbol{z}\_{o\prime}, \boldsymbol{x}\_{o\prime}, t\_{o})
$$

$$
\overrightarrow{\boldsymbol{U}}\_{2} = \overrightarrow{\boldsymbol{U}}(\boldsymbol{z}\_{o} + \left(\overrightarrow{\boldsymbol{U}}\_{1}\right)\_{\boldsymbol{z}} \times \delta t / 2, \,\,\mathbf{x}\_{o} + \left(\overrightarrow{\boldsymbol{U}}\_{1}\right)\_{\boldsymbol{x}} \times \delta t / 2, \,\,t\_{o} + \delta t / 2\right)
$$

$$
\overrightarrow{\boldsymbol{U}}\_{3} = \overrightarrow{\boldsymbol{U}}\Big{(}\boldsymbol{z}\_{o} + \left(\overrightarrow{\boldsymbol{U}}\_{2}\right)\_{\boldsymbol{z}} \times \delta t / 2, \,\,\mathbf{x}\_{o} + \left(\overrightarrow{\boldsymbol{U}}\_{2}\right)\_{\boldsymbol{x}} \times \delta t / 2, \,\,t\_{o} + \delta t / 2\Big{)}\tag{13}
$$

$$
\overrightarrow{\boldsymbol{U}}\_{4} = \overrightarrow{\boldsymbol{U}}\Big{(}\boldsymbol{z}\_{o} + \left(\overrightarrow{\boldsymbol{U}}\_{3}\right)\_{\boldsymbol{z}} \times \delta t, \,\,\boldsymbol{x}\_{o} + \left(\overrightarrow{\boldsymbol{U}}\_{3}\right)\_{\boldsymbol{x}} \times \delta t, \,\,t\_{o} + \delta t\Big{)}
$$

The Eulerian fluid velocity <sup>→</sup> *U*<sup>1</sup> at a point (*zo*, *xo*) at time *to* is computed by a spatial interpolation with an upstream biquadratic polynomial. <sup>→</sup> *U*<sup>2</sup> is computed in two steps. First a quadratic time interpolation of the nine velocities, used for the spatial interpolation, must be carried out because <sup>→</sup> *U*<sup>2</sup> is determined at the intermediate time *to* + δ*t*/2. Then the spatial upstream biquadratic interpolation is performed with the velocities interpolated in time. The time interpolation is performed with the velocity fields at time *to* <sup>−</sup> <sup>δ</sup>*t*, *to*, and *to* <sup>+</sup> <sup>δ</sup>*t*. The velocities <sup>→</sup> *<sup>U</sup>*<sup>3</sup> and <sup>→</sup> *<sup>U</sup>*<sup>4</sup> are computed in the same way as <sup>→</sup> *<sup>U</sup>*<sup>2</sup> and <sup>→</sup> *U*<sup>1</sup> , respectively. The computation of the ionic Lagrangian velocity is carried out in the same manner with the three flux contributions (<sup>→</sup> *J* ± *Co* <sup>=</sup> *<sup>c</sup>*±<sup>→</sup> *U*, → *J* ± *Fi* <sup>=</sup> <sup>−</sup>*D*±<sup>→</sup> ∇*c*±, → *J* ± *Em* <sup>=</sup> <sup>∓</sup>*D*±*c*±<sup>→</sup> ∇Φ) and with the ion concentration; here *Co*, *Fi* and *Em* mean convective, Fickian and electro-migration contribution, respectively. The interpolated flux is divided by the interpolated concentration to determine the four <sup>→</sup> *U* ± .

#### **3. Numerical Results**

#### *3.1. Comparison with Published Results*

To validate the numerical method, the results were compared to those published in References [18,28,41]. To my knowledge, no analytical solution is available in the literature. Several analyses of the polarization layer structure in an over-limiting regime have been carried out by assuming 1D ion transfer at steady state (see Zaltzman et al. [14] and references therein). With these assumptions, the electric field is a solution of the 2nd kind Painlevé equation which does not admit a simple analytical solution.

Figure 3 shows the mean current density through the boundary layer (*I*) and through the membrane surface (*Im*) as a function of the potential drop. Three series of computations were carried out. Each point of the curves corresponds to the value of the current density averaged in time at a prescribed potential drop. In the Ohmic and limiting regimes, the potential drop increases by a step of 2.5 until it reaches a value of 21. In the over-limiting regime, a first series is carried out with *D*± = 1, <sup>ν</sup> = 1/3 <sup>×</sup> 10<sup>−</sup>3. Beyond the value of 22, the initial condition is the unsteady solution at the previous value of ΔΦ. Concerning the other series, the initial condition at a given potential drop is the unsteady solution of the first series with the same value of ΔΦ in order to reduce the transitory period. The series 1 and 2 (*D*<sup>+</sup> = *DNa*<sup>+</sup>, *<sup>D</sup>*<sup>−</sup> = *DCl*−, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup>−3) will be used for the comparison with published numerical results and theoretical ones. The 3rd series is computed with *D*<sup>±</sup> = 1 and ν = 10−3. The time average of the current density is computed after a stabilization period of the current signal. By stabilization it is not meant that steady state is recovered. The transitory period can take several hundred time steps, especially after a jump of the potential. The current density is averaged over one to two thousand time steps. In Sections 3.2 and 3.3, the instability structure and the dynamic streamlines are analyzed after the transitory period.

These series show that the over-limiting threshold values lie between ΔΦ = 21 and 22 (0.57 V). This value, slightly sensitive to the diffusion coefficient and the length ratio ν, is higher than the value found in the case of periodic conditions [31,37]. In Nikonenko et al. [28], an expression of the local limiting current density, in the case of a laminar flow between two plates, was proposed. Averaged over the domain length, the limiting current density *I*lim has the following expression (Equation (17.56) in Newmann et al. [47]):

$$I\_{\rm lim} = \frac{2}{L\_{\rm x}} \left[ 1.47 \left( \frac{L\_{\rm x}^2 \langle \ell L\_{\rm z} \rangle}{L\_{\rm z} D} \right)^{1/3} \right] \tag{14}$$

where *D* = 2*D*+*D*−/(*D*<sup>+</sup> + *D*−). The expression (14) is equivalent to Equation (3b) in Nikonenko et al. [28] because the cation transfer numbers *T*<sup>1</sup> and *t*<sup>1</sup> are equal to 1 and 0.5, respectively. In Equation (14) *Lx* is used instead of *Lm* because the potential drop is imposed over *Lx*. In the same idea, *Uz* is the inlet velocity averaged over the computational domain width. This relation gives a value of 9 and 12.4 for the series 1 and 2, respectively. Therefore, it seems from Figure 3 that *I* is the relevant quantity in the present model.

Let us compare the cation profiles obtained in the present work and those published in References [28,41]. The authors investigated the ion transfer between an anion exchange membrane and a cation one. In Urtenov et al. [41], the two membranes were non-ideal, i.e., co-ions could cross the membranes. The potential drop is imposed at the outer membrane surfaces. In this publication, Figure 5 shows the static streamlines and the map of the cation concentration with the parameter values *Lx* = *Lm* = <sup>1</sup> mm, *Lz* = <sup>2</sup> mm, <sup>ν</sup> = <sup>10</sup>−4, *Uz* = <sup>8</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> ms<sup>−</sup>1, *<sup>D</sup>*+/<sup>−</sup> = *DNa*<sup>+</sup>/*Cl*−. Therefore, if a part of the domain is modeled, i.e., *Lx* = 0.3 mm, one obtains <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, which corresponds to series 2 if *D*+/<sup>−</sup> = *DNa*+/*Cl*−. However, the value of the potential threshold is close to 1 V instead of 0.57 V. Two reasons may explain this difference in value: (1) The potential gradient is located mainly in high resistance regions, i.e., where the ion concentration is low. So, the transverse gradient takes place in the two polarization layers [28,41]. Assuming the problem is symmetric, the potential drop

through the cationic polarization layer would be equal to about half of the overall potential drop. (2) The membranes have a non-zero electrical resistivity and the potential drop is imposed outside the electro-dialytic cell. These two conditions lead to an increase in the potential threshold. However, it also seems that the potential threshold depends on *Lm* [28,41]. At the threshold, the authors [41] found that the outlet cation concentration was equal to 0.6 at *x* = *xc*+=0.6 = *Lm*/10 (100 μm) from the cation exchange membrane, which corresponds to *Lx*/3 in the present work at ΔΦ = 22 (see Figure 4A). Their values ΔΦ = 1.8 V, *<sup>c</sup>*<sup>+</sup> = 0.6 at *xc*+=0.6 = *Lm*/7 (140 <sup>μ</sup>m) correspond to ΔΦ = <sup>30</sup> = 1.4 <sup>×</sup> 21.5, *c*<sup>+</sup> = 0.6 at *xc*+=0.6 = 0.5*Lx* (150 μm) found here. The counterion concentration color map shown in Figure A2 (Nikonenko et al. [28]) come from simulations carried out with the same type of membrane cell (*Lx* <sup>=</sup> *Lm* <sup>=</sup> 0.5 mm, *Lz* <sup>=</sup> <sup>2</sup> mm, *Uz* <sup>=</sup> <sup>8</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> ms<sup>−</sup>1, *<sup>D</sup>*+/<sup>−</sup> <sup>=</sup> *<sup>D</sup>Na*+/*Cl*<sup>−</sup> ) and with several inlet ion concentrations. Two cases may be used for comparison: <sup>ν</sup> = 2 <sup>×</sup> 10−<sup>4</sup> (Figure A2a) and <sup>ν</sup> = 6 <sup>×</sup> 10−<sup>5</sup> (Figure A2b). In these Figures ΔΦ = 1.43ΔΦ∗∗ with ΔΦ∗∗ = 0.7V. In the present work, *Lm* = 0.5 mm and ΔΦ = 1.43ΔΦ∗∗ correspond to *Lx* = 0.15 mm and ΔΦ = 31, respectively. In Figure A2a,b, *c*+/*co* = 0.6 at *xc*+=0.6 = 0.6*Lx* (90 μm), which is close to *xc*+=0.6 = 0.5*Lx* (75 μm) obtained with <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup>−<sup>3</sup> (Figure 4B). The color bar range in terms of dimensionless concentration is assumed to be identical to the range in Figure 5 (Urtenov et al. [40]). The values of *xc*+=0.6 are given in Table 1.

**Table 1.** Width cation boundary layer defined by c<sup>+</sup> = 0.6 at the outlet of the domain. 1st row: published values (Nikonenko et al. [28] and Urtenov et al. [41]). 2nd row: this work.


Kwak et al. [18] found a linear correlation between the vortex size δ<sup>v</sup> and the potential drop ΔΦ<sup>v</sup> over δv: δv/*Lm* = a ΔΦ<sup>2</sup> v/*Uz* 1/3 + b where *Uz* is the mean inlet velocity (see Table 2). With their numerical model, Urtenov et al. [41] verified this linear behavior. In the present work, as ΔΦ increases from 23 to 30 (Figure 4), at the outlet ΔΦ<sup>v</sup> and δv/*Lm* increase from 19.4 to 25.7 and from 0.04 to 0.07, respectively. If the velocity is averaged over *Lx*, the values of a and b are very close to those found by Urtenov et al. [41] (Table 2). Similar values are found if the velocity is averaged over *Lm* (a = 0.0215, b = 0.107).

**Table 2.** Parameter values of the correlation δv/*Lm* = a ΔΦ<sup>2</sup> v/*Uz* 1/3 + b. 1st and 2nd columns: published values. 3rd columns: this work.


### *3.2. Hydrodynamic Instability*

The Nikonenko research group has mainly investigated the structure of the polarization layer. However, the instability modes are briefly studied, and the cluster dynamics and structure are not addressed. On the contrary, in this section, the key role of the force spot upstream of the vortex front in the vortex shape and in the birth of counter-clockwise vortices and clusters is shown as an example. In Section 3.3, the role of SSL bending, induced by the layer ( → *<sup>F</sup>Tot*)*<sup>z</sup>* <sup>&</sup>lt; 0 at the vortex rear, in the fluid flux through the outer edge of the unsteady vortex layer will be also pointed out. The investigation is performed with *D*<sup>±</sup> = 1 and ν = 1/310−3. Another mode of instability found with *D*<sup>±</sup> = 1 and ν = 10−<sup>3</sup> is briefly discussed in the Appendix A.

#### 3.2.1. Periodic Mode Flow (ΔΦ = 22)

Before going further, the main results described in the previous article [37] about the spatial distribution of the electric and pressure forces and of the fluid kinetic energy are recalled briefly. In forced flow, the periodic mode close to the instability threshold is very similar to the marginal instability because the rolls are stable during their travel to the outlet. In periodic systems, the steady instability is characterized by two counter-rotating vortices. The non-equilibrium EDL has a minimum width in the downward flow region and presents a hump at the upward one. The hump extends between the two vortex centers. *Q* decreases mainly in the hump because its transverse mean value remains constant. It is recalled that the non-equilibrium EDL gathers the ESC region and the quasi-equilibrium EDL. Two force layers along the membrane induce a longitudinal flow. In the first layer placed in the quasi-equilibrium EDL, the pressure force is dominant, whereas the electric one is dominant in the second layer located around the local maximum of electric charge density (*Q*max) in the ESC region. Two spots of kinetic energy are located in the second layer and close to the ESC hump. A third spot is located between the two vortices in the upward region due to the convergence of the fluid streamlines. The pressure and potential disequilibrium on both sides of the vortices is explained first by the transverse profile of charge density, which controls the variation of the potential gradient, and also by the equilibrium between the electric force and the pressure one in the transverse direction leading to the quadratic dependence of the pressure on the transverse potential gradient.

Figure 5 displays the vortices (SSL) and the color map of the two components of the total force ( → *FTot*), which is the sum of the electric ( → *FEl*) and pressure force ( → *FPr*). All vortices rotate clockwise. The region of upward fluid motion defines the front of the vortex. Inversely, the region of downward fluid motion defines the rear of the vortex. The space between two successive vortices, i.e., the inter-vortex region, replaces the counter-clockwise vortex in the periodic model. The circumstances in which a counter-clockwise vortex occurs this space will be discussed. The minimum width of the non-equilibrium EDL is located at the rear of the vortex and the hump takes place around the vortex front (upward flow). The two spots of force (Figure 5A) are equivalent to the spots shown in Reference [37] but with a much higher intensity. In these spots, the pressure force is dominant.

**Figure 5.** Color plot of: (**A**) ( → *<sup>F</sup>Tot*)*z*; (**B**) ( → *<sup>F</sup>Tot*)*x*. Black lines: static streamlines starting from the inlet. ΔΦ <sup>=</sup> 22, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, *<sup>D</sup>*<sup>±</sup> = 1.

A layer and a spot of the same force sign, located at two different transverse positions, may be linked as shown in Figure 5A. This means that a continuous path of the same force sign may go through the front or the rear of the vortex. This controls the inclination of the static streamlines. For this reason, the streamlines are oriented downstream of the front. When there is no link, the front is vertical. At the

rear, in the vortex located at *z* ≈ 5.3 for example, the spot ( → *<sup>F</sup>Tot*)*<sup>z</sup>* <sup>&</sup>gt; 0 is linked to the layer of positive sign, whereas, in the upstream vortices rear (*z* ≈ 5 and *z* ≈ 4.7), the pair spot/layer ( → *<sup>F</sup>Tot*)*<sup>z</sup>* <sup>&</sup>lt; 0 are linked. This explains the small difference in orientation of the streamlines at the rear of the two vortices and the small deviation from the exact periodicity. Below, the effect of these spot/layer links on the form of the vortices and their influence on the vortex dynamics will be more clearly shown and described.

Figure 5B shows the transverse component of the total force ( → *<sup>F</sup>Tot*)*<sup>x</sup>* . ( → *<sup>F</sup>El*)*<sup>x</sup>* is always negative as a consequence of the positive sign of *Q*. The transverse pressure force is opposed to the electric one in the front region, whereas it is negative in the rear region. Three spots are located around the vortex front. At the vortex front, the pressure force dominates (( → *<sup>F</sup>Tot*)*<sup>x</sup>* <sup>&</sup>gt; 0). The spot is located inside the ESC region. On both sides of the front, the electric force is dominant (( → *<sup>F</sup>Tot*)*<sup>x</sup>* <sup>&</sup>lt; 0) suppressing the upward fluid motion. The bottom of these two spots are located at the outer edge of the ESC region [37].

In Figure 6, the fluid kinetic energy (*Ec*) also presents three spots of high magnitude. However, the kinetic energy at the front is much smaller than when there is no forced flow [37], and it induces a smaller roll size. This also induces a periodic trajectory of the anions (see Section 3.3.3). A disequilibrium of *Ec* is also visible on both sides of the front. As when there is no forced flow, the two *Ec* spots are located at the transverse position of *Q*max: the fluid is accelerated in the longitudinal direction by the electric force and reaches a maximum kinetic energy when *Q*max and the ESC width begin to decrease and increase, respectively.

**Figure 6.** Color plot of the hydrodynamic kinetic energy *E*c. Black lines: static streamlines starting from the inlet. ΔΦ <sup>=</sup> 22, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, *<sup>D</sup>*<sup>±</sup> = 1.

Contrary to [28,41], another periodic mode is also observed at ΔΦ = 23. The period of the vortex structure is composed of a vortex pair formed from a large and a small vortex (Figure 7). The vortex pair size is seen to increase with time. Two consecutive vortices begin to form a pair when the upstream vortex (V1) reaches the position *z* = 3.5. The pair formation consists of three steps. During the first one, the downward vortex (V2) comes closer to the downstream pair so that its size decreases. At the two ends of V2, the streamlines are vertical close to the membrane because the negative layer and the negative spots ( → *<sup>F</sup>Tot*)*<sup>z</sup>* are linked at these ends (for an example, see a similar situation at *<sup>z</sup>* <sup>=</sup> 4.7 in Figure 8D). In the second step, when V2 reaches its minimum size, V1 begins to move closer to V2. The negative layer and the negative spot ( → *<sup>F</sup>Tot*)*<sup>z</sup>* get closer at the rear of V1, whose size decreases during its displacement. In the third step, the negative spot ( → *<sup>F</sup>Tot*)*<sup>z</sup>* between V1 and V2 flattens and the two layers ( → *<sup>F</sup>Tot*)*<sup>z</sup>* <sup>&</sup>lt; 0 at the rear of V1 join leading to vertical streamlines near the interface. At the same time, at the rear of V2, the positive spot and the positive layer ( → *<sup>F</sup>Tot*)*<sup>z</sup>* also join. The length of V2 increases and its height decreases accordingly. Figure 7 shows the third and final step.

**Figure 7.** Color plot of ( → *<sup>F</sup>Tot*)*z*. Black lines: static streamlines. ΔΦ <sup>=</sup> 23, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, *<sup>D</sup>*<sup>±</sup> = 1.

3.2.2. Non-Periodic Mode Flow (ΔΦ = 27)

As the potential drop increases, the periodicity disappears gradually. At ΔΦ = 27 and higher, an unsteady instability mode is clearly established and the rolls associate to form clusters. The simulations with *<sup>D</sup>*+/<sup>−</sup> = *DNa*<sup>+</sup>/*Cl*−, ΔΦ <sup>=</sup> 27, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, show vortex groups independent from one another as observed in Nikonenko et al. [28] and no roll exchange occurs. Inside the groups, dissociation and association are observed. With *D*± = 1, the spatial periodicity of the vortex group is not observed, and the clusters are a transitory structure. They become one vortex (Figure 8) or break into several rolls (Figure 9). It must be noticed that the definition of envelope and cluster is not the same as in Nikonenko et al. [28]. In the latter article, the envelope is the black SSL around a vortex group and a cluster is a counter-rotating vortex pair. In the current work, a cluster is a vortex group wrapped in an envelope (e.g., red SSL in Figure 8C). It does not contain counter-clockwise rolls but several clockwise ones.

In Figures 8 and 9, the static streamlines (SSL) are associated with the color map of ( → *<sup>F</sup>Tot*)*<sup>z</sup>* . The blue, the green and the red SSL represent clockwise rolls, counter-clockwise rolls and the cluster envelopes, respectively. The time evolution of the cluster begins at δ*t* = 1600δ*t* (with δ*t* = 10<sup>−</sup>4). For the sake of simplicity, this time is used as the initial one and the time is expressed in terms of δ*t* in the comments (with *t*' = *t*/δ*t* − 1600).

First, let us consider Figure 8. In Figure 8A (*t*' = 0), a counter-clockwise roll (V1/−) lies between the two clockwise vortices V1/+ and V2/+ located at *z* = 3.4 and 3.8, respectively. The signs + and − refer to the rotation direction. In Figure 8A at the rear of V2/+, the positive layer ( → *<sup>F</sup>Tot*)*<sup>z</sup>* is the pressure force layer. It can be recalled that the negative layer ( → *<sup>F</sup>Tot*)*<sup>z</sup>* is composed of the pressure layer next to the membrane and the electric one just above. Owing to the negative pressure force spot between V1/+ and V2/+, the black SSL at the rear of V1/+ and at the front region of V2/+ are vertical. V1/− is surrounded by the fluid film separating the vortex from V1/+, V2/+ and the membrane. V1/− rotates counter-clockwise due to the friction with the fluid film, despite the negative sign of ( → *<sup>F</sup>Tot*)*<sup>z</sup>* at the roll bottom. That is the pressure force acts against the fluid motion at the bottom of V1/−. The inter-vortex region can be viewed as a cavity with a slip condition at the walls. This condition is induced by the fluid film.

In Figure 8B (*t*' = 20), the volume of V1/− increases and its bottom moves to the membrane, so that it is clearly located inside the positive layer ( → *<sup>F</sup>Tot*)*<sup>z</sup>* . However, the thin fluid film continues to separate V1/− from the membrane. The vortex V1/+ front gets closer to V2/+ which does not move. The negative pressure spot ( → *<sup>F</sup>Tot*)*<sup>z</sup>* begins to join the layer of the same sign through the front of V2/+. A new vortex appears in the tail of V2/+.

**Figure 8.** Temporal evolution of the cluster structure. (**A**) *t*' = 0; (**B**) *t*' = 20; (**C**) *t*' =40; (**D**) *t*' = 60: (**E**) *t*' = 80; (**F**) *t*' = 100. Color plot: ( → *<sup>F</sup>Tot*)*z*. Black lines: static streamlines starting from the inlet. Blue streamlines: clockwise rotating rolls, green streamlines: counter-clockwise rotating rolls, red streamlines: cluster envelope. ΔΦ <sup>=</sup> 27, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, *<sup>D</sup>*<sup>±</sup> = 1.

**Figure 9.** Temporal evolution of the cluster structure. (**A**) *t*' = 0; (**B**) *t*' = 20; (**C**) *t*' =40; (**D**) *t*' = 60: (**E**) *t*' = 80; (**F**) *t*' = 100. Color plot: ( → *<sup>F</sup>Tot*)*z*. Black lines: static streamlines starting from the inlet. Blue streamlines: clockwise rotating rolls, green streamlines: counter-clockwise rotating rolls, red streamlines: cluster envelope. ΔΦ <sup>=</sup> 27, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, *<sup>D</sup>*<sup>±</sup> = 1.

In Figure 8C (*t*' = 40), the negative layer ( → *<sup>F</sup>Tot*)*<sup>z</sup>* inside V1/+ connects with the negative spot ( → *<sup>F</sup>Tot*)*<sup>z</sup>* inside V1/−. Therefore, a continuous path occurs along the membrane joining V1/+ and V2/+. This coincides with the apparition of an envelope around V1/+ and V2/+. Now, V2/+ clearly constitutes two rolls: V2a/+ (head) and V2b/+ (tail). During the merging process (Figure 8A–E), V2/+ remains at the same place, so that the distance from the downstream vortex increases. The very long non-equilibrium EDL becomes unstable and a hump of charge appears at the front of V2b/+. The SSL forming the rear of V1/+ and the front of V2/+ respectively are clearly tilted. Moving downward, V1/− gets into contact with the membrane as explained in Figure 9. However, V1/− remains isolated from V1/+ and V2/+ by the fluid film forming the cluster.

In Figure 8D (*t*' = 60), V1/− and the front of V2/+ disappear owing to the decrease in the charge hump in the non-equilibrium EDL. The spot ( → *<sup>F</sup>Tot*)*<sup>z</sup>* <sup>&</sup>gt; 0 and the trace of the negative spot ( → *<sup>F</sup>Tot*)*<sup>z</sup>* on both sides of the previous front of V2/+ remain. As soon as the charge hump disappears, the former spot moves upstream and joins the positive spot ( → *<sup>F</sup>Tot*)*<sup>z</sup>* of V1/+. In Figure 8E (*t*' <sup>=</sup> 80), V2/+, i.e., the charge hump moves downstream this time. The tail is composed of two vortices rotating clockwise separated by a secondary envelope. In Figure 8F (*t*' = 100), the tail, constituting one roll, is separated completely from the head of V2/+ and the cluster disappears.

The concept of thin film can be extended to the cluster envelope. This film separates the clockwise vortices from the membrane and the surrounding electrolyte. Along the membrane surface, it spreads over a width smaller than the non-equilibrium EDL one. The presence of this film means that a kind of upstream electro-osmotic velocity is imposed between the vortex group and the membrane. However, the immobility of the cluster is not related to this velocity but to the stability of the non-equilibrium EDL.

In the second series (Figure 9), V1/− is initially located between the two rolls V1/+ (*z* = 2.1) and V2/+ (*z* = 2.35). Contrary to Figure 8A, a continuous path of negative ( → *<sup>F</sup>Tot*)*<sup>z</sup>* extends along the membrane over the vortex group. Despite this path, no SSL envelops the vortex group because the tops of V1/+ and V2/+ are too separated, allowing V1/− to spread far from the membrane interface. V1/+ and V2/+ get closer (Figure 9B) and the SSL between the vortices are more and more tilted owing to the profile of ( → *<sup>F</sup>Tot*)*<sup>z</sup>* (see comments in Figures <sup>7</sup> and 8).

At the top of V1/+ and V2/+, the imposed main flow drives the fluid located close to the vortices in the downstream direction by means of the viscous force. At the same time, in the inter-vortex region, the longitudinal pressure force pushes the fluid upstream. Therefore, at the outer edge of the inter-vortex region, the SSL changes direction and the two points of the same SSL, on both sides of the neck, are close to each other. This behavior is all the more pronounced given that the two vortices are close to each other and the magnitude of the negative pressure spot ( → *<sup>F</sup>Tot*)*<sup>z</sup>* increases. Therefore, the SSL join at the upper region of the roll group at *t*' = 40 (Figure 9C). This induces a displacement of V1/− to the membrane and the presence of an envelope around the vortex group which becomes a cluster. At the same time, the roll (V0/+ located at *z* = 1.6) upstream of the cluster moves closer with a decreasing volume. From *t*' = 40 to 80 (Figure 9C–E), the tail of V2/+ extends as the downstream vortex goes away and a charge hump appears at *x* = 2.5 as explained in Figure 8C. Later, the cluster breaks, V0/+ fuses with V1/+, V2/+ and its tail becomes independent, and the negative layer ( → *<sup>F</sup>Tot*)*<sup>z</sup>* in the inter-vortex region is transformed into a spot. It seems that a cluster including three clockwise vortices is not stable or may only exist during a very short time.

### *3.3. Dynamic Streamlines of Fluid Flow and Ionic Transfer*

A fully Lagrangian approach to the fluid flow and ion transfer is necessary to investigate the mixing mechanisms in unsteady flow. In this section, the dynamic streamlines (DSL) display the real trajectories of the fluid and ions, and point out that the vortices visualized by the SSL are fictitious. In this way, for instance, the DSLs show that, in a periodic regime, the vortices and the region far from the membrane, i.e., the diffusion layer containing high ion concentration, are linked by a mechanism of fluid flux through the vortices. However, in an aperiodic regime, the unsteady size of the inter-vortex regions and the vortices promotes fluid exchange with the diffusion layer and enhances the efficiency of the mixing. The same parameter values as in Section 3.2 are used.

#### 3.3.1. Fluid DSL in Periodic Regime (ΔΦ = 22)

Figure 10 shows four dynamic streamlines at five different times. The line source is located inside a vortex at *zo* = 4.7. The vortices move with a longitudinal velocity (*Uroll <sup>z</sup>* ) of 36. As in Section 3.2, the time value will be expressed in terms of time step and one assumes that initially *t*' = 0.

**Figure 10.** Fluid trajectories. Black lines: static streamlines starting from the inlet. Red lines: dynamic streamlines starting from *z*o = 4.7 (*x*o = 0.005, 0.045, 0.075, 0.115). Color plot: cation concentration. The red arrow represents the rotation direction of the vortices and the position of the vortex at time *t*' (the trajectories start at *t*' = 0 inside this vortex). (**A**) *t*' = 5; (**B**) *t*' = 70; (**C**) *t*' =95; (**D**) *t*' = 195; (**E**) *t*' = 275. ΔΦ <sup>=</sup> 22, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, *<sup>D</sup>*<sup>±</sup> = 1.

Let us consider the starting point *xo* = 0.005. This point is located between the membrane surface and the *Ec* maximum inside the *Ec* spot. At this point, the transverse equilibrium between the electric force and the pressure force is fulfilled. Initially, *Uz* = −43 and *Ux* = −1. At *t*' = 10, the fluid begins to move upward as it comes closer to the vortex front. Then *Ux* and *Uz* increases and decreases, respectively. Consequently, at *t*' = 20 (*x* = 0.033), *Ux* reaches a maximum value of 49 and *Uz* = 24 < *Uroll <sup>z</sup>* . Then *Uz* increases and becomes greater than *Uroll <sup>z</sup>* . *Uz* decreases first because the fluid element moves away from the *Ec* spot (see Figure 6) and then increases owing to the velocity inlet boundary condition. At *t*' = 95 (Figure 10C), the fluid is just above the vortex (*x* = 0.17, *z* = 5) with *Ux* ≈ 0 and *Uz* = 74.5. During the upward motion, the DSL deviates from the vortex front for two reasons: first *Uz* < *Uroll z* for a period of time and then the angle θ between the membrane surface and the front is lower than 90o. With regards to the 2nd reason, if the fluid element moves during a time interval Δ*t* upward over the distance Δ*x* = *Ux*Δ*t*, the distance to the front becomes Δ*z*<sup>2</sup> = *Ux*Δ*t* × tan(π/2 − θ). Finally, the distance from the fluid element to the front is <sup>Δ</sup>*<sup>z</sup>* <sup>=</sup> <sup>Δ</sup>*z*<sup>1</sup> <sup>+</sup> <sup>Δ</sup>*z*<sup>2</sup> with <sup>Δ</sup>*z*<sup>1</sup> <sup>=</sup> *Uroll <sup>z</sup>* − *Uz* Δ*t* and Δ*z* cancels if *Uz* = *Uroll <sup>z</sup>* + *Ux* × tan(π/2 − θ). If *Ux* > 0, *Uz* must be much greater than *Uroll <sup>z</sup>* if θ is small. As *Q* is positive, ( → *<sup>F</sup>El*)*<sup>x</sup>* is always negative. In the front region ( → *<sup>F</sup>Pr*)*<sup>x</sup>* <sup>&</sup>gt; 0 (see Figure 5B). During the fluid rising, ( → *<sup>F</sup>Pr*)*<sup>x</sup>* and ( → *<sup>F</sup>El*)*<sup>x</sup>* decrease but the pressure contribution is higher than the electric one. Above the vortex, the force balance occurs. However, the force magnitude is very small because *Q* is very small. On the other side of the vortex, ( → *<sup>F</sup>Pr*)*<sup>x</sup>* is negative. Therefore, the downward motion of the fluid is accelerated due to the two negative force contributions. The fluid enters the vortex at *t'* = 190 far from the center of rotation and close to the membrane (*x* = 0.03, *z* = 5.47—Figure 10D) (the fluid entrance can be explained by the above relation Δ*z* = Δ*z*<sup>1</sup> + Δ*z*<sup>2</sup> with *Ux* < 0 and θ>π/2). The fluid element continues to move to the membrane surface until the transverse force equilibrium is fulfilled (*Ux* = 0). Owing to the vortex displacement in the downstream direction, the *Ec* spot (with *Uz* < 0) moves towards the fluid element. This induces the displacement of the fluid to the front with *Ux* ≈ 0, so that at *t*' = 230/240, the fluid element leaves the vortex and then moves upward.

Figure 10 shows that the fluid motion looks like a trochoid around the same vortex. The time and the spatial period are approximately 210 and 0.8, respectively. As *xo* increases to 0.035, the fluid element is less and less accelerated in the ejection region at the vortex front and spends less time in it. Therefore, the maximal position in the transverse direction decreases, the loop of the trochoid appears more and more early and its size decreases. So that the period decreases to 70 and 0.39 in time and space, respectively.

If the starting point is located at (*xo* = 0.045 , *zo* = 4.7), the streamline shape shows another kind of trochoid because the starting point is located above the center of rotation. Initially *Ux* ≈ −5 and *Uz* = −9. The fluid moves upward then downward as explained above. This time, the DSL enters the vortex at *t*' = 70 at the position (*x* = 0.1, *z* = 4.97) close to the top of the vortex (Figure 10B). During the downward displacement, |*Ux*| and *Uz* decrease because the fluid element moves closer to the center of rotation. *Uz* becomes lower than *Uroll <sup>z</sup>* and at *t*' = 115 the center of rotation catches up with the fluid element. Then the fluid is pushed downstream with an increasing value of *Uz*. So, the fluid element stays above the center of rotation. As the vortex front moves closer to the position of the fluid element, the latter leaves the vortex at *t*' = 160 (*x* = 0.071, *z* = 5.13). The time and spatial periodicity are around 114 and 0.37, respectively.

As *xo* increases to 0.09, the transverse amplitude of the trochoid decreases and the periodicity remains equal to about 110 and 0.37 in time and space, respectively. This time periodicity corresponds to the time needed for the vortices to move over the distance between two successive vortex fronts. Outside the vortex (*xo* > 0.09), the trajectory is in fact a prolongation of another trajectory beginning upstream (*z <sup>o</sup>*, *x <sup>o</sup>*) with *x <sup>o</sup>* lower than the rotation center and at an earlier time given the periodic shape of the trajectory. As a last observation, Figure 10 shows that if *xo* = 0.075, the trajectory is linear (*Ux* ≈ 0). The fluid element remains close to the rotation center and moves with the same velocity as the roll.

Therefore, the figure shows that the real width of the "hydrodynamic instability" layer is much larger than the vortex one. Its outer edge is in contact with the diffusion layer. In the reference frame of the vortices, the trajectories form concentric circles. The rotation center, placed at *x* ≈ 0.75, is above the SSL vortex one. The fluid exchange between the vortex layer and the diffusion layer is controlled by an inward flux at the rear and an outward one at the front. The trajectories passing close to the membrane are the most efficient for the fluid exchange. This explains the presence of the mixing layer whose width is higher than the SSL vortex one.

Now if *zo* = 4.55, i.e., upstream of the vortex front, the DSL shape is not always a trochoid because the fluid element can remain inside the inter-vortex region as shown in Figure 11. As an example, let us consider the starting position (*xo* = 0.015, *zo* = 4.55). The point is located inside the *Ec* spot. Initially, *Ux* = −6 and *Uz* = 87 > *Uroll <sup>z</sup>* <sup>=</sup> *<sup>U</sup>spot <sup>z</sup>* . Therefore, as the DSL gets closer to the front vortex, the fluid element begins to move upward, and *Uz* decreases to a value smaller than *Uroll <sup>z</sup>* (see the discussion about Figure 10). Consequently, at *t*' = 35 (Figure 11B), the vertical line defined by *Ux* = 0 (black dashed line) and the DSL meet at *z* = 4.67 (the fluid element is at *x* = 0.063). This line moves with the same velocity as the vortices. Therefore, the fluid element passes on the other side of this line and moves towards the membrane with an increasing velocity. Then it enters the *Ec* spot, so that *Uz* reaches a value of 85 at *t*' = 75 (see Figure 11C). Therefore, the fluid element remains in the inter-vortex region (at *t*' = 275 *z* = 5.5, *x* = 0.06—Figure 11E). This behavior is observed until *xo* = 0.065. If 0.075 ≤ *xo* ≤ 0.095, i.e., at the outer edge of the inter-vortex region, the trochoid shape is observed but the DSL does not wind around the same roll. Most of the time, the fluid element moves with *Uz* << *Uroll <sup>z</sup>* and with |*Ux*| small enough for the fluid element to enter the rear of the upstream vortex close to the membrane. This is shown for example with *xo* = 0.095 in Figure 11C (*x* = 0.014, *z* = 4.76). The loop is located at *z* = 4.8 (Figure 11D). Figure 11E shows that the same fluid element enters the even upstream vortex a second time at *t*' = 255 (*x* = 0.024, *z* = 5.06). The loop at *z* = 5 belongs to the trajectory beginning at *xo* = 0.085. Therefore, the inter-vortex is not connected to the diffusion layer and it has a small role in the mixing process.

#### 3.3.2. Fluid DSL in Aperiodic Regime (ΔΦ = 27)

In Figure 12 the source is located inside the vortex V0. Initially, the front of V0 and of the downstream vortex V1 are at *z* = 3.7 and *z* = 4.25, respectively. Until *t*' = 25, the volume of V0 is unchanged and its height is around 0.18. Then its volume decreases during the merging process with V1. The front of V0 moves with a longitudinal velocity of 66.7 until *t*' = 85. Furthermore, *Uroll <sup>z</sup>* is smaller than the fluid inlet velocity at *x* = 0.18 (*Uz* = 74). The front of V1 remains at the same place in the time interval. Two starting points, (*xo* = 0.01, *zo* = 3.9) and (*xo* = 0.1, *zo* = 3.9), are chosen for the discussion because they give rise to two kinds of motion.

In the first case, the fluid element moves first upstream along the membrane to the vortex front, then moves upward in the front region at *t*' = 25 (*x* = 0.01, *z* = 3.8) and finally passes V0 at *t*' = 85, as explained in Section 3.3.1. The height difference between V0 and V1 makes the SSL tilted to the upward direction for a long time, such that at *t*' = 85 the fluid element is at a transverse position much greater than the height of V0. Reaching V1, the fluid leaves the instability region (Figure 12C). The electric force acts no more. The fluid moves mainly in the longitudinal direction with *Uz* ≈ <sup>120</sup> >> *Uroll <sup>z</sup>* . Therefore, the addition of the high jump at *t*' = 25 and of the increase of the height of V1 leads to the ejection of the fluid outside the vortex layer.

On the contrary, the fluid remains close to the instability layer if *x*<sup>o</sup> = 0.1. The source is located above the center of rotation of V0. Therefore, the fluid element leaves the vortex outside the fluid ejection region as explained in Section 3.3.1 (see Figure 12A). The fluid element is just above V0 (Figure 12B) and then above V1 at *t*' = 135 because *Uz* > *Uroll <sup>z</sup>* (Figure 12C). It reaches its maximum transverse position *x* = 0.25 at *z* = 4.77. From *t*' = 135 to 190, as the fluid element overtakes V1, the height of V1 increases to 0.18 and its shape becomes symmetric. The static streamlines are more and more directed downward near the rear region of V1 (see Section 3.2.1 ΔΦ = 23) and the fluid element is driven towards the membrane. For the same reason as explained in the case of ΔΦ = 22, the fluid element must enter the vortex V1 and the loop of the trochoid begins at *t*' = 210 (*x* = 0.03, *z* = 5.2).

**Figure 11.** Fluid trajectories. Black lines: static streamlines starting from the inlet. Red lines: dynamic streamlines starting from *z*o = 4.55 (*x*o = 0.015, 0.045, 0.085, 0.095). Black dashed line: curve defined by *Ux* = 0. Color plot: cation concentration. The red arrow represents the rotation direction of the vortices and the position of the vortex at time t' (this vortex is downstream the trajectories source at *t*' = 0). (**A**) *<sup>t</sup>*' <sup>=</sup> 10; (**B**) *<sup>t</sup>*' <sup>=</sup> 35; (**C**) *<sup>t</sup>*' <sup>=</sup> 80; (**D**) *<sup>t</sup>*' <sup>=</sup> 155; (**E**) *<sup>t</sup>*' <sup>=</sup> 275. ΔΦ <sup>=</sup> 22, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, *<sup>D</sup>*<sup>±</sup> = 1.

As *x*o decreases from 0.1, the DSL reaches a higher transverse position during the fluid excursion, i.e., the fluid explores the diffusion layer more deeply, and it enters the vortex layer at a higher distance from the source position. This means that the high transverse amplitude of the trajectory and the large dispersion of the DSL allow a fluid exchange between different distant vortices in contrast to the periodic regime.

**Figure 12.** Fluid trajectories. Black lines: static streamlines starting from the inlet. Red lines: dynamic streamlines starting from *z*o = 3.9 (*x*o = 0.01, 0.08, 0.09). Color plot: cation concentration. Red numbers: vortex numbering. Vortex n◦0 is the vortex in which the trajectories source is located at *t*' = 0. (**A**) *t*' = 55; (**B**) *<sup>t</sup>*' <sup>=</sup> 85; (**C**) *<sup>t</sup>*' <sup>=</sup> 135; (**D**) *<sup>t</sup>*' <sup>=</sup> 190; (**E**) *<sup>t</sup>*' <sup>=</sup> 245. ΔΦ <sup>=</sup> 27, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, *<sup>D</sup>*<sup>±</sup> = 1.

Figure 13 displays the fluid trajectories when the source is located at *z*<sup>o</sup> = 2.7 between V0/+ and V1/+ (see Figure 13A). Between these rolls, a counter-clockwise roll (V/−) is located at *x* ≈ 0.04. During the rolls' displacement, the inter-vortex region between V0/+ and V1/+ expands and V/− disappears at *t*' = 55. At *t*' = 135 (Figure 13C), V1/+ stops to move downstream (*z f ront* = 3.2), V0/+ moves closer to V1/+ and the roll V/− appears one more time between these two vortices, which merge at *t*' = 225 (Figure 13D).

**Figure 13.** Fluid trajectories. Black lines: static streamlines starting from the inlet. Red and green lines: dynamic streamlines starting from *z*o = 2.7 (*x*o = 0.01, 0.025, 0.05, 0.085). Color plot: cation concentration. Red numbers: vortex numbering. Vortex n◦0 is the vortex upstream the trajectories source location at *t*' = 0. Vortex n◦1 is the vortex downstream the source. (**A**) *t*' = 5; (**B**) *t*' = 45; (**C**) *t*' = 135; (**D**) *t*' = 220; (**E**) *<sup>t</sup>*' <sup>=</sup> 245. ΔΦ <sup>=</sup> 27, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, *<sup>D</sup>*<sup>±</sup> = 1.

Let us consider the first trajectory *x*<sup>o</sup> = 0.01 (red trajectory). This starting point is located inside the *E*<sup>c</sup> spot close to the membrane, i.e., in the fluid film below V/−. If *t*' < 80 (*z* < 3), the fluid trajectory is similar to the trajectory beginning at *xo* = 0.015, shown in Figure 11. The fluid element remains in the thin film surrounding V/− in the inter-vortex region. However, the fluid element passes above V/− during the first upward period and leaves for a very short time the fluid film at the outer edge of the inter-vortex region. This means that V/− is also isolated from the neck, a behavior which the SSL cannot not reveal (see Figures 8 and 9). At *t*' = 100 (*x* = 0.12, *z* = 3), the fluid element reaches the region where *Ux* is small. However, the fluid element does not cross the boundary *Ux* = 0 because the inter-vortex region expands. Therefore, the fluid element continues to move upward and leaves

the inter-vortex region with an increasing longitudinal velocity, which becomes higher than *Uroll <sup>z</sup>* (see Figure 13C). As the size of V1/+ increases when *t*' > 135, the DSL enters the V1/+ rear at *t*' = 175 as explained previously and then the fluid trajectory makes a loop inside V1/+.

For the source point *x*<sup>o</sup> = 0.025 (green trajectory), the fluid element is initially inside V/−. The vortex is isolated from the neighborhood. So, until *t*' = 55, the fluid element remains inside V/− and rotates counter-clockwise two times. At *t*' = 60, V/− disappears; therefore, until *t*' = 125 (see Figure 13C), the fluid trajectory oscillates in the transverse direction as explained in Section 3.3.1. At *t*' = 135, V/− appears, and the fluid element makes two counter-clockwise rotations until *t*' = 180 (*x* = 0.012, *z* = 3.2). The inter-vortex region between V0/+ and V1/+ is so small that the rear of V0/+ passes the fluid element during its upward motion at *t*' = 200 (*x* = 0.045, *z* = 3.25). The fluid consequently moves upstream inside V0/+ to the front (Figure 13D), then it moves upward during the merging of V0/+ with V1/+. The volume of the new vortex increases. Therefore, as the fluid element moves with *Uz* > *Uroll z* at *t*' = 245, it must enter the new vortex by the rear later on and must then begin a loop motion.

The two last Figures show that several kinds of trajectories follow one another owing to the vortex dynamics. Therefore, the trajectories playing a small role in the exchange between the vortex layer and the diffusion layer have a finite lifetime. The inter-vortex regions exchange fluid with the vortices. As in the periodic mode, the passage through the vortices plays a key role in this exchange and enhances the ion exchange. The vortices trap the fluid elements by means of their changing shape. The counter-clockwise vortices do not exchange fluid through its boundary.

#### 3.3.3. Ionic DSL in the Periodic Regime (ΔΦ = 22)

In this section, the ion transfer in the aperiodic regime is not investigated because the anions are simply ejected from the vortex region if ΔΦ > 24, and the cations follow the same kind of trajectories as in periodic regime.

Let us consider the anionic motion when the starting longitudinal position is *zo* = 4.7, i.e., inside a vortex as shown in Figure 14. It is recalled that the ion velocity (<sup>→</sup> *U*±) is computed by means of the ion flux density (<sup>→</sup> *J* ±) (see Section 2.2.2).

Consider first the case *xo* = 0.005. At this distance from the membrane, the negative Fickian contribution ( → *J* − *Fi*) *<sup>x</sup>* balances the positive electric one ( → *J* − *Em*) *x* . The convection ( → *J* − *Co*) *<sup>z</sup>* and the concentration gradient ( → *J* − *Fi*) *<sup>z</sup>* push the anions to the vortex front. At the beginning of the rising period (Figure 14A), the anions are pushed transversely by convection because ( → *J* − *Fi*) *<sup>x</sup>* <sup>=</sup> <sup>−</sup>( → *J* − *Em*) *x* . Then around the position (*x* = 0.1, *z* = 4.8), the electro-migration contribution becomes dominant with ( → *J* − *Fi*) *<sup>x</sup>* ≈ −( → *J* − *Co*) *x* . Finally at *t'* = 95 Figure 14B, the ions reach the maximal transverse position with the balance ( → *J* − *Fi*) *<sup>x</sup>* ≈ −( → *J* − *Em*) *x* ·( → *J* − *Co*) *<sup>x</sup>* <sup>≈</sup> 0) and with *<sup>U</sup>*<sup>−</sup> *<sup>z</sup>* > *Uroll <sup>z</sup>* . In the longitudinal direction, the ions move by convection. During the downward motion (see Figure 14C), the ions enter the same vortex at *t* ≈ 210 and get closer to the vortex rotation center. The ions are pushed to the membrane by convection and diffusion and with *U*− *<sup>z</sup>* < *Uroll <sup>z</sup>* . At *t* ≈ 210, the anions arrive at the same transverse position as the rotation center. ( → *J* − *Co*) becomes negligible and ( → *J* − *Fi*) ≈ −( → *J* − *Em*), i.e., the anions arrive at a stagnation point. Therefore, after a while the anions cross over to the other side of the rotation center, then they begin an upward motion and leave the vortex a second time. As *x*o increases from 0.035 to 0.085, the ions deviate from the ejection region as they leave the vortex. Consequently, the transverse amplitude of their trajectory decreases, and the cusps move upstream.

**Figure 14.** Fluid and anionic trajectories. Black lines: fluid static streamlines starting from the inlet. Red lines: anionic dynamic streamlines starting from *z*o = 4.7 (*x*o = 0.005, 0.035, 0.055, 0.085). Color plot: cation concentration. The red arrow represents the rotation direction of the vortices and the position of the vortex at time *t*' (the trajectories start at *t*' = 0 inside this vortex). (**A**) *t*' = 10; (**B**) *t*' = 95; (**C**) *t*' = 210; (**D**) *<sup>t</sup>*' <sup>=</sup> 250; (**E**) *<sup>t</sup>*' <sup>=</sup> 275. ΔΦ <sup>=</sup> 22, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3, *<sup>D</sup>*<sup>±</sup> = 1.

Figure 15A displays the anionic DSL when the anion source is located between two vortices (*z*<sup>o</sup> = 4.55). Two families of trajectory can be observed. The first one is identical to the trajectories shown in Figure 14 if *x*o lies between 0.035 and 0.055. However, the stagnation point is replaced by an elongated and tilted loop because the anions reach a transverse position lower than the rotation center one and move to the vortex front by diffusion ( → *J* − *Fi*) *<sup>z</sup>* and convection ( → *J* − *Co*) *z* . At higher values of *x*o, *U*− *<sup>z</sup>* is always positive because the DSL never enter the vortex owing to higher fluid velocity and to the positive contribution of ( → *J* − *Fi*) *z* . The sign of ( → *J* − *Fi*) *<sup>z</sup>* is positive because the anion position is always downstream of the location of the high ion concentration region at the vortex rear. The consequence is that the anions move faster than the vortices along the membrane if *xo* ≥ 0.075.

**Figure 15.** Fluid and ionic trajectories. Red lines: (**A**) anionic dynamic streamlines starting from *z*o = 4.55 (*x*o = 0.005, 0.045, 0.055, 0.065, 0.1). The red arrow represents the position of the vortex at time *t*' (this vortex is downstream the trajectories source at *t*' = 0); (**B**): cationic dynamic streamlines starting from *z*o = 3, *x*o = [0.23, 0.35]. Black lines: fluid static streamlines starting from the inlet. Color plot: cation concentration. *<sup>t</sup>*' <sup>=</sup> 275. ΔΦ <sup>=</sup> 22, <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup>−<sup>3</sup>

Therefore, the longitudinal Fickian contribution is of great importance to the Lagrangian anionic motion. The two last Figures also show that despite the electric force directed upward, the anions remain in the layer explored by the fluid DSL. At the outer-edge of this layer, the Fickian diffusion balances the electro-migration and keeps the anions leaving the region.

In regards to the cations (Figure 15B), the dynamic streamlines are computed with the vertical source (*z*<sup>o</sup> = 3, 0.23 ≤ *xo* ≤ 0.35). This time, the electro-migration flux is directed downward like the Fickian one. ( → *J* + *Em*) *<sup>x</sup>* is the main contribution to the cation motion towards the membrane, and Figure 15B shows that the cations move monotonically. At the initial time, cation motion is mostly longitudinal (*U*<sup>+</sup> *<sup>z</sup>* <sup>&</sup>gt; *<sup>U</sup>*<sup>+</sup> *<sup>x</sup>* ). As they get closer to the membrane, the DSL become vertical because ( → *J* + *Em*) *x* increases and becomes dominant. The convective contribution has the same order of magnitude as the electro-migration one in the region around the outer edge of the vortex layer and it leads to the increase of the cationic flux. The DSL form four groups at the vicinity of the membrane because they reach the membrane through the inter-vortex region and the rear of the rolls, i.e., regions where the transverse cationic flux is important.

#### **4. Conclusions and Summary**

In this work, the analysis of electro-convective instabilities in a laminar boundary layer along a cation exchange membrane is presented. The periodic and aperiodic modes were studied. In particular, for the latter mode, vortex association/dissociation and fusion, as well as time evolution of cluster and counter-clockwise vortex, were investigated. The influence of the instability structure and dynamics on the fluid and ion transfer was also investigated by means of a Lagrangian approach. The aperiodic mode is sensitive to the value of the parameters *D*± and ν. Despite this observation, the present study was performed mainly with *<sup>D</sup>*<sup>±</sup> = 1 and <sup>ν</sup> = 1/3 <sup>×</sup> <sup>10</sup><sup>−</sup>3.

The vortex structure and dynamics were analyzed by considering the spatial distribution of the longitudinal component of the electric and pressure force. This distribution is composed of a layer along the membrane inside the non-equilibrium EDL and of two spots on both sides of the vortex fronts and above the force layer. These two structures undergo a change of sign at the front and the rear of the vortices. The link between a spot and a layer of negative sign plays a key role in the vortex dynamics and in the birth of counter-clockwise vortices and clusters. This link bends the streamlines in the upward direction and induces a neck at the outer edge of the inter-vortex region where the fluid is pushed downward by the forced flow. The counter-clockwise vortices occur in the inter-vortex region over which the negative spot is located. This vortex is isolated from the neighboring clockwise vortices and the membrane by a thin fluid film coming from the neck. The clockwise vortices associate into a cluster when the neck closes, and a continuous path of negative longitudinal force occurs over the vortex group along the membrane. In this case, a thin film also wraps the vortex group, which is therefore separated from the membrane and the boundary layer. The counter-clockwise vortices are excluded from the clusters in contact with the membrane during the clustering process. The upstream fluid motion in the film along the membrane plays a similar role as the electro-osmotic boundary conditions.

The real motion of the fluid, as well as the ions, must be computed by considering the unsteadiness of the velocity and flux fields. In the periodic mode, the fluid trajectories passing through a vortex, computed with the SSL, form a trochoid. The SSL vortices are fictitious. In the vortex reference frame, these DSLs form rolls whose rotation center is above the SSL vortex center and whose size is much greater. The real vortex layer is therefore in contact with the diffusion layer containing high ion concentration and it controls the mixing layer expansion. This is not the case for the inter-vortex region. Surprisingly, the anions follow the same kind of trajectory with roughly the same transverse amplitude. The Fickian diffusion plays a key role in the loop formation and in the transverse amplitude. The outer edge of the region explored by the anions corresponds to the balance between the electro-migration and the Fickian diffusion fluxes in the transverse direction. Further away from the instability threshold, the electric field is high enough for the anions to be ejected from the vortex layer.

For the aperiodic mode, the real instability layer visualized by the DSL is also much larger than the SSL rolls. Its size increases with ΔΦ, imposing the upward displacement of the diffusion layer. In this case, fluid exchange is performed between distant SSL vortices. The transverse amplitude of the trajectories is high enough for the mean Lagrangian velocity to be much greater than the SSL roll one. The size increase allows the rolls to attract and trap the fluid elements. It is also observed that the inter-vortex regions contribute to the mixing process. The variation of the inter-vortex length destabilizes the fluid trajectory leading to the fluid ejection outside the SSL vortex layer and initiating a trochoid-like trajectory.

**Funding:** This research received no external funding.

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

#### **Abbreviations**



#### **Appendix A. Influence of the Length Ratio on the Instability Modes**

If (*D*<sup>±</sup> = 1, <sup>ν</sup> = <sup>10</sup><sup>−</sup>3), the periodic instability, without roll pair formation, occurs in the range 23 <sup>≤</sup> ΔΦ ≤ 25. Aperiodic mode appears at ΔΦ = 27. For this potential drop value, rolls of different sizes cross the domain. Vortex association and dissociation are also observed from time to time. However, a new kind of periodicity takes place at ΔΦ ≥ 30. Figure A1 shows the time evolution of the vortex structure at ΔΦ = 35.

Three steps have been observed. Every other roll is in the same phase and the rolls between are in another one, that is two phases coexist. The Figures show clearly that, at the vortex front and rear, the SSL orientation is related to the unsteady link between a spot and a layer of longitudinal force of the same sign. Let us consider the smallest rolls at *z* = 3.1 and 4.3 in Figure A1A. As a first step, these rolls, of initial size 0.2 × 0.11 in the *z* and *x* direction, respectively, move to *z* = 3.2 and 4.4 and increase in volume (Figure A1B). At the beginning, the negative spot ( → *<sup>F</sup>Tot*)*<sup>z</sup>* is transformed into a tilted layer in the downstream inter-vortex region. As the volume increases (Figure A1B), the SSL turn progressively to the upstream direction at the roll front, and a counter-clockwise vortex in the upstream inter-vortex region appears (*z* = 3.05 and 4.3). It disappears in the next step. This first step lasts 60δ*t*. In a second step (Figure A1C), the roll volume decreases and becomes elongated and symmetric (*z* = 3.4 and 4.6). The elongation occurs at the rear, the front remaining roughly at the same place. The roll dimension is about 0.6 × 0.16. This step lasts 30δ*t*. In the last step, the roll volume continues to decrease so that the rolls return to their initial form. This third step lasts 80δ*t*. As every second roll is in step N, the other rolls are in step N − 1 or N + 1. However, in the same roll family, a small phase shift in the shape evolution may occur from one roll to another and the growth amplitude may also differ. But the mechanism remains the same.

**Figure A1.** Temporal evolution of the vortices. (**A**) *t*' = 0; (**B**) *t*' = 60; (**C**) *t*' = 90. Color plot: ( → *<sup>F</sup>Tot*)*z*. Black lines: static streamlines starting from the inlet. ΔΦ = 35, ν = 10<sup>−</sup>3, *D*<sup>±</sup> = 1.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
