Simulation and Experimental Study of Ion Concentration Polarization Induced Electroconvective Vortex and Particle Movement

Ion concentration polarization (ICP) has been widely applied in microfluidic systems in pre-concentration, particle separation, and desalination applications. General ICP microfluidic systems have three components (i.e., source, ion-exchange, and buffer), which allow selective ion transport. Recently developed trials to eliminate one of the three components to simplify the system have suffered from decreased performance by the accumulation of unwanted ions. In this paper, we presented a new ICP microfluidic system with only an ion-exchange membrane-coated channel. Numerical investigation on hydrodynamic flow and electric fields with a series of coupled governing equations enabled a strong correlation to experimental investigations on electroconvective vortices and the trajectory of charged particles. This study has significant implications for the development and optimization of ICP microfluidic and electrochemical systems for biomarker concentration and separation to improve sensing reliability and detection limits in analytic chemistry.


Introduction
Ion concentration polarization (ICP) is a fundamental electrokinetic phenomenon resulting from selective ion transport through ion-exchange films, such as cation-exchange membranes (CEM) or anion-exchange membranes (AEM), which allow only cations or anions to pass [1][2][3]. Conventional electromembrane processes (e.g., electrodialysis) follow a typical system configuration in which spacers (i.e., channels) and ion-exchange membranes are alternatively stacked between electrodes to operate under two vertically applied electric and hydrodynamic flow fields [4,5]. After the electric potential is applied across CEM/AEM, the ion-depleted/enriched diffusion layer on the anodic side of the membrane is triggered by selective ion transport accompanying the ICP. A thicker ion depletion layer that develops beyond a diffusion-limited thickness leads to an electroconvective instability. The thicker layer also results in the reduction of system performance owing to the limitation of ion diffusivity and the degradation of membrane selectivity by highly depleted or accumulated ions on the membrane surface [3,6,7]. Recently, microfluidic systems have facilitated a well-defined steady-state electroconvection, owing to the advantages of spatial controllability in microfluidic systems, the application of in micro-and nanofluidic-scale pre-concentration bio-agents [8,9], the separation of charged particles [10,11], and fluidic mixing [12]. ICP microfluidic systems have various configurations according to device compartments and forces applied to control target particles for corresponding applications ( Figure 1 and Table 1).
A traditional ICP microfluidic system which comprises three components (ICP-3C) was depicted as an equivalent circuit model (Figure 1a). The first component is a source, and the second is a buffer where an anode and a cathode are generally placed. The third is an ion-exchange channel compartment (e.g., CEM or negatively charged nanochannel junction), which allows the movement of counter-ions (i.e., cations in the case of CEM) between the two compartments [2,8,13,14]. The source and buffer compartments can be described as resistances R S and R B , respectively, and the diode represents the ionexchange compartment with an internal resistance of R I (R I R S ). Target particles (i.e., ions, deoxyribonucleic acids (DNAs), proteins, or cells) are placed adjacent to the anodic side of the ion-exchange compartment alongside the source compartment. The buffer compartment collects the relocated counter-ion after it passes through the ion-exchange compartment. The ICP-3C system is classified into two applied categories according to external force fields: a single field (electric field (EF) only, ICP-3C-EF, Figure 1b) [8,13] and a coupled field (EF and hydrodynamic flow fields (HF), ICP-3C-EF/HF, Figure 1c) [14,15]. The ICP-3C-EF system develops an ion depletion layer next to the anodic side of the ionexchange compartment and generates a spatially redistributed electric field (EF) [8,9,16]. The redistributed EF forms an electro-osmotic flow alongside the source compartment and an electroconvective vortex (EV) on the anodic side of the ion-exchange compartment. The equilibrium of electrophoretic and drag forces on particles results in the enrichment or separation of bio-agents to detect concentrated or separated bio-agents in specific areas of the channel [8,13]. The ICP-3C-EF/HF system, in which an additional HF allows continuous control of target particles, interrupts the flow of nano-and micro-sized particles and lets them be collected or detected from the channel for further processing [10,14,15]. In spite of their usability, ICP-3C systems suffer from an accumulation of ionic species (i.e., transported counter-ions and newly generated chemicals by the electrode) in the buffer compartment, causing deterioration of the perm-selectivity of the ion-exchange compartment. Additional washing is therefore required to maintain stable operation.  Ion concentration polarization (ICP)-microfluidic system. (a) Equivalent circuit of ICP-3C system consisting of source, buffer, and ion-exchange compartments. ICP-3C with an application of (b) electric field (EF) (ICP-3C-EF) [8,13] and (c) electric and hydrodynamic flow fields (HF) (ICP-3C-EF/HF) [10,14,15]. (d) Equivalent circuit of ICP-2C system, consisting of source and ion-exchange compartment. ICP-2C with an application of (e) electric field (ICP-2C-EF) [17,18] and (f) electric and hydrodynamic flow fields (ICP-2C-EF/HF) [11,12,19]. S , B , and I indicate resistances of source, buffer, and ion-exchange compartment, respectively.  Ion concentration polarization (ICP)-microfluidic system. (a) Equivalent circuit of ICP-3C system consisting of source, buffer, and ion-exchange compartments. ICP-3C with an application of (b) electric field (EF) (ICP-3C-EF) [8,13] and (c) electric and hydrodynamic flow fields (HF) (ICP-3C-EF/HF) [10,14,15]. (d) Equivalent circuit of ICP-2C system, consisting of source and ion-exchange compartment. ICP-2C with an application of (e) electric field (ICP-2C-EF) [17,18] and (f) electric and hydrodynamic flow fields (ICP-2C-EF/HF) [11,12,19]. R S , R B , and R I indicate resistances of source, buffer, and ion-exchange compartment, respectively.
Instead of ICP-3C systems, a two-compartment configuration (ICP-2C) having only source and ion-exchange channel compartments was suggested ( Figure 1d) [12,[17][18][19]. The new configuration directly connects the anode and cathode via the source compartment. Ions in the source compartment can freely move toward their corresponding electrodes, and the counter-ions move through the ion-exchange compartment rather than the source compartment, owing to the high conductivity of the ion-exchange compartment toward the counter-ion (R S R I ). The counter-ion movement creates a higher electric field near the anodic side of the ion-exchange compartment, enabling protein concentration [17,18]. When integrated with additional HF, the ICP-2C-EF/HF system increased throughput to extract the target particle in the preconcentration process [11,12,19]. Compared with the ICP-3C systems which suffer from ion accumulation in the buffer compartment, the ICP-2C system has distinct merits in terms of its continuous discharge of the waste buffer and simple structures [8,17]. However, the ICP-2C system has not been fully optimized with numerical analyses by applying the three coupled governing equations of Nernst-Plank, Navier-Stokes, and Poisson (Table 2). This study develops an axial two-dimensional (2D) symmetry model of the ICP-2C-EF/HF system with coupled governing equations. Numerical predictions of the electroconvective vortex and movement of charged particles were confirmed experimentally under various hydrodynamic and electric fields. Significant implications are expected for the development of ICP microfluidic systems to control biomarkers by improving sensing reliability and detection limits in analytic chemistry. Table 2. Summary of current numerical study of ICP-microfluidic system.

Type
Force Field Boundary Condition for Ion-Exchange Compartment Ref.

ICP-3C
EF Fixed volumetric charge [30] Fixed volumetric charge [31] Fixed surface charge and electric potential [32,33] EF/HF Fixed surface charge [22] Fixed volumetric charge [34] Fixed surface charge and fixed electric potential [35] ICP-2C EF Fixed volumetric charge [31] EF/HF Fixed surface charge This work All studies applied a two-dimensional (2D) model. The volumetric and surface boundary conditions are constraints for area and line in the 2D model.

Experiment Setup
The experimental data were generated using a microfluidic device for the ICP-2C-EF/HF system, whose configuration, fabrication, and operation have been described and demonstrated in previous works [11,19]. Figure 2a and Figure S1 show the operation and the schematic depiction of the movement of ions (0.1× phosphate-buffered saline) and particles (0.2 mL −1 , 1-µm diameter; Invitrogen, Carlsbad, CA, USA). The device was designed to contain two reservoirs connected via an ion-exchange compartment, with Nafion-coated (Sigma-Aldrich, St. Louis, MO, USA) ion-concentration channel 500 µm wide and 120 µm high (hydraulic radius ≈ 100 µm). The platinum electrodes were connected to a source meter (Model 2400, Keithley Instruments, OH, USA) to regulate the electric potential at an inter-electrode spacing of 1 cm.

Experiment Setup
The experimental data were generated using a microfluidic device for the ICP-2C-EF/HF system, whose configuration, fabrication, and operation have been described and demonstrated in previous works [11,19]. Figures 2a and S1 show the operation and the schematic depiction of the movement of ions (0.1× phosphate-buffered saline) and particles (0.2 mL −1 , 1-µ m diameter; Invitrogen, Carlsbad, CA, USA). The device was designed to contain two reservoirs connected via an ion-exchange compartment, with Nafioncoated (Sigma-Aldrich, St. Louis, MO, USA) ion-concentration channel 500 μm wide and 120 μm high (hydraulic radius ≈ 100 μm). The platinum electrodes were connected to a source meter (Model 2400, Keithley Instruments, OH, USA) to regulate the electric potential at an inter-electrode spacing of 1 cm.

Theoretical Model and Simulation
Details of the simulation settings (e.g., governing equations, simulation models, and simulation methods) are described in the supporting information ( Figure S2 and Table  S1). Three coupled governing equations of Nernst-Plank, Navier-Stokes, and Poisson

Theoretical Model and Simulation
Details of the simulation settings (e.g., governing equations, simulation models, and simulation methods) are described in the supporting information ( Figure S2 and Table  S1). Three coupled governing equations of Nernst-Plank, Navier-Stokes, and Poisson were applied to an axial 2D symmetry model of the ICP-2C-EF/HF system. The governing equations were solved numerically on finite elements using COMSOL Multiphysics (v.5.2, COMSOL, Inc., Stockholm, Sweden). For a simplified expression, the cylindrical geometries, ion concentration, electric potential, and flow velocity are normalized as described: (r, z) = r r 0 , 10 z z 0 where r 0 (=100 µm) is the equivalent radius of the ion-permselective coated channel, and z 0 (=1000 µm) is the length scale of the channel length. C 0 (=0.1 mM) is the bulk concentration in the electrolyte solution with a diffusivity of 1.6 × 10 −9 m 2 /s for both cations and anions. ψ 0 (=1 V) is the reference electric potential to achieve an electric field of the same order as the experimental condition (=10-30 V/cm). U 0 (=0.106 mm/s, Q 0 = 2 µL/min) is the reference average velocity.

ICP and Electroconvective Vortex in ICP-2C-EF/HF System
In the ICP-2C-EF/HF, ICP was initiated at the entrance of the ion-exchange compartment, the ion concentration channel (IC channel) forming an electric double layer (EDL), owing to the lower electric resistance at 25 • C of nafion on the surface (3.65 × 10 5 /S) than that of the bulk solution (5.4 × 10 10 /S) [36] (Figure 2a). The electric field (−∇ψ) can be considered as two; one −∇ψ through the IC-channel (−∇ψ C ) and the other −∇ψ near the IC-channel surface (−∇ψ S ) formed by the space charge described by Poisson's equation. A higher ψ generates electro-osmotic instability because of the limitation of the ionic diffusion transport, which forms a plateau IC profile that leads to an electroconvective vortex (Figure 2b). This profile regulates a local electric field profile (i.e., −∇ψ S ) that generates an electrophoretic and dielectrophoretic force on a charged particle. The development of an electroconvective vortex is also regulated by the hydrodynamic flow (U), which works as shear flow, hydrodynamically suppressing the growth of the vortex and reducing ICP by its continuous ion supply. The U additionally acts as a drag force on a particle [10,37]. The fully coupled interactions redistribute the flow and electric fields and regulate particle movement. Figure 3 shows experimental streamlines as a function of U in and ψ in ( U in , ψ in ). A circular depletion layer (i.e., electroconvective vortex) is observed at the entrance of the IC channel. Lower U in ((1,3) Figure 3b) facilitates the growth of large-sized electroconvective vortices, while a higher U in ((3,3) Figure 3c) suppresses growth. Compared to the conventional electromembrane process in which electroconvective vortices are sequentially developed along the channel [38], in the ICP-2C-EF/HF, particles follow a stabilized path after experiencing a dramatic detour at the entrance by a spatially controlled electroconvective vortex.   Figure 4 shows IC distribution in the IC channel, polarized along the IC cha face (Figure 4a,b). The typical EDL and diffusion layer ( ̅ = 1.01, Figure 4b,c) noticed at the tight entrance of the IC channel ( ̅ = 1.01). The IC gradient then g decreased and stabilized, forming a flattened concentration profile ( ̅ = 1.11, Fig  similar to that where electroconvective vortices were observed (i.e., EV zone in F [38]. Two distinct ionic distributions of circular and parabolic profiles (Figure found, similar to the velocity profiles local at the IC channel entrance and global a channel. A circular profile of ions similar to the velocity profile developed by th convective vortex at the channel entrance will be discussed later. A parabolic p ions is observed along the channel where the Hagen-Poiseuille flow was develo profile similarity indicates that the ionic distribution and flow velocity were coup ure 4c,d shows the changes in the IC profile at the entrance ( ̅ = 1.01) and stabiliz ( ̅ = 1.11) along the ̅ -axis, respectively, as a function of flow velocity ( ̅ in ) an potential ( ̅ in ). Under a fixed ̅ in , the increase in ̅ in led to a thicker depletion la a steeper concentration gradient, because more ions migrated through the IC-cha face as ̅ in increased. Along the same principle, the increase in ̅ in increased flux through the channel, resulting in a thinner depletion layer with a lower conc gradient.  Figure 4 shows IC distribution in the IC channel, polarized along the IC channel surface (Figure 4a,b). The typical EDL and diffusion layer (z = 1.01, Figure 4b,c) could be noticed at the tight entrance of the IC channel (z = 1.01). The IC gradient then gradually decreased and stabilized, forming a flattened concentration profile (z = 1.11, Figure 4b,d) similar to that where electroconvective vortices were observed (i.e., EV zone in Figure 1b) [38]. Two distinct ionic distributions of circular and parabolic profiles (Figure 4a) were found, similar to the velocity profiles local at the IC channel entrance and global along the channel. A circular profile of ions similar to the velocity profile developed by the electroconvective vortex at the channel entrance will be discussed later. A parabolic profile of ions is observed along the channel where the Hagen-Poiseuille flow was developed. The profile similarity indicates that the ionic distribution and flow velocity were coupled. Figure 4c,d shows the changes in the IC profile at the entrance (z = 1.01) and stabilized point (z = 1.11) along the r-axis, respectively, as a function of flow velocity (U in ) and electric potential (ψ in ). Under a fixed U in , the increase in ψ in led to a thicker depletion layer with a steeper concentration gradient, because more ions migrated through the IC-channel surface as ψ in increased. Along the same principle, the increase in U in increased the mass flux through the channel, resulting in a thinner depletion layer with a lower concentration gradient. Micromachines 2021, 12, x FOR PEER REVIEW 7 of 12  Figure 5 shows a normalized electric field (−∇ ̅ ) and normalized flow velocity ( ̅ ) along the ̅ -axis at the entrance of the IC channel ( ̅ = 1.01). −∇ ̅ is maximized at the IC channel surface (−∇ ̅ S ) due to the field-focusing effect [39] and decreases sharply towards the center of the IC channel (−∇ ̅ C ). Notably, −∇ ̅ C is slightly higher than the electric field formed only from ion-free electric potential (−∇ ̅ ion−free ), because −∇ ̅ C is determined by both the application of electric potential ( ̅ in ) and the positive net charge. A higher ̅ leads −∇ ̅ C to approach −∇ ̅ ion−free (Figure 5a-enlarged) due to the reduced electric field by a higher convective flow (i.e., higher mass flux). Figure 5b presents velocities of developing flow near the entrance region. As simulated, flow field in the IC channel was determined by the hydrodynamic and electro-osmotic flows. ̅ near the surface of the IC channel is dominated by −∇ . A higher −∇ ̅ S leads to high electro-osmotic instability near the channel surface and the development of an electroconvective vortex, which will be discussed in the next section. ̅ of the central region of the IC channel (near ̅ = 0) is regulated by −∇ , and ̅ has a relatively flat profile under a higher −∇ ̅ in at ̅ = 1, which is known to be the profile of electro-osmotic flow.  Figure 5 shows a normalized electric field (−∇ψ) and normalized flow velocity (U) along the r-axis at the entrance of the IC channel (z = 1.01). −∇ψ is maximized at the IC channel surface (−∇ψ S ) due to the field-focusing effect [39] and decreases sharply towards the center of the IC channel (−∇ψ C ). Notably, −∇ψ C is slightly higher than the electric field formed only from ion-free electric potential (−∇ψ ion−free ), because −∇ψ C is determined by both the application of electric potential (ψ in ) and the positive net charge. A higher U in leads −∇ψ C to approach −∇ψ ion−free (Figure 5a-enlarged) due to the reduced electric field by a higher convective flow (i.e., higher mass flux). Figure 5b presents velocities of developing flow near the entrance region. As simulated, flow field in the IC channel was determined by the hydrodynamic and electro-osmotic flows. U near the surface of the IC channel is dominated by −∇ψ S . A higher −∇ψ S leads to high electro-osmotic instability near the channel surface and the development of an electroconvective vortex, which will be discussed in the next section. U of the central region of the IC channel (near r = 0) is regulated by −∇ψ C , and U has a relatively flat profile under a higher −∇ψ in at U in = 1, which is known to be the profile of electro-osmotic flow. Micromachines 2021, 12, x FOR PEER REVIEW 8 of 12

Development of Electroconvective Vortex
The formation of an electroconvective vortex was simulated according to various ̅ and ̅ conditions ( Figure 6). The vortex mainly develops near the surface of the IC channel entrance where ICP is significantly developed (Figure 6a). The streamline was plotted to visualize the electroconvective vortex as a function of ̅ (Figure 6b), the deformed elliptical shape deflected to the cathode. As ̅ increases, the elliptical shape tends to be inclined toward the cathode and decreased in size, suppressed by higher ̅ . The electroconvective vortex can also be regarded as a forced vortex having a maximum vorticity in the vortex core. Its vorticity and coordinates are shown in Figure 6c,d. The increase in both ̅ and ̅ leads to stronger vorticity and the early onset of the vortex, moving the position of the vortex core to the channel surface and to the IC channel entrance. ̅ not only suppresses the growth of the electroconvective vortex but also accelerates its vorticity with high shear stress. A higher ̅ results in a higher electro-osmotic instability and increases the early onset of the vortex.

Development of Electroconvective Vortex
The formation of an electroconvective vortex was simulated according to various U in and ψ in conditions ( Figure 6). The vortex mainly develops near the surface of the IC channel entrance where ICP is significantly developed (Figure 6a). The streamline was plotted to visualize the electroconvective vortex as a function of U in (Figure 6b), the deformed elliptical shape deflected to the cathode. As U in increases, the elliptical shape tends to be inclined toward the cathode and decreased in size, suppressed by higher U in . The electroconvective vortex can also be regarded as a forced vortex having a maximum vorticity in the vortex core. Its vorticity and coordinates are shown in Figure 6c,d. The increase in both U in and ψ in leads to stronger vorticity and the early onset of the vortex, moving the position of the vortex core to the channel surface and to the IC channel entrance. U in not only suppresses the growth of the electroconvective vortex but also accelerates its vorticity with high shear stress. A higher ψ in results in a higher electro-osmotic instability and increases the early onset of the vortex.

Movement of Charged Particles
The motion of charged particles in flow and electric fields is generally regulated by hydrodynamic drag (F D ), electrophoretic force (F EP ), and dielectrophoretic force (F DEP ) (see supporting information for details) [10,11]. To understand how these forces manipulate the motion of charged particles, a simulation was performed to track their trajectory and applied forces (Figure 7). Flows formed a flattened profile along the channel under a higher ψ in , due to the electro-osmotic flow (color of the surface in Figure 7a). The motion of the particle was mainly regulated by the electric field in the channel (−∇ψ C ) rather than that near the channel surface (−∇ψ S ) at its entrance (Region A in Figure 7a). ∇ψ S is much higher than ∇ψ C , because F EP generated by −∇ψ C near the channel entrance delays the particle movement toward the channel to form the larger circular trajectory. It placed a particle beyond the influence of −∇ψ S and ∇ψ ∇ (i.e., dielectrophoretic force). A higher ψ in (i.e., higher F EP ) slows particle movement toward the channel entrance (z = 1), leading to a shift in the particle trajectory toward the center of the IC channel (r = 0). After particles formed a circular trajectory, they were pushed toward the IC channel surface until the particle trajectory was parallel to the IC channel (Figure 7a,b), where particles experienced no net force (Figure 7c).

Movement of Charged Particles
The motion of charged particles in flow and electric fields is generally regulated by hydrodynamic drag ( ), electrophoretic force ( ), and dielectrophoretic force ( ) (see supporting information for details) [10,11]. To understand how these forces manipulate the motion of charged particles, a simulation was performed to track their trajectory and applied forces (Figure 7). Flows formed a flattened profile along the channel under a higher ̅ , due to the electro-osmotic flow (color of the surface in Figure 7a). The motion of the particle was mainly regulated by the electric field in the channel (−∇ ̅ ) rather than that near the channel surface (−∇ ̅ ) at its entrance (Region A in Figure 7a). ∇ is much higher than ∇ ̅ , because generated by −∇ ̅ near the channel entrance delays the particle movement toward the channel to form the larger circular trajectory. It placed a particle beyond the influence of −∇ ̅ and ∇ ̅ • ∇ (i.e., dielectrophoretic force). A higher ̅ (i.e., higher ) slows particle movement toward the channel entrance ( ̅ = 1), leading to a shift in the particle trajectory toward the center of the IC channel ( ̅ = 0). After particles formed a circular trajectory, they were pushed toward the IC channel surface until the particle trajectory was parallel to the IC channel (Figure 7a,b), where particles experienced no net force (Figure 7c). Position and color of sphere indicate where particle travels for corresponding time and magnitude of particle velocity, respectively. Path and color of line indicate trajectory and velocity throughout its entire course. The particle trajectory simulation started at (1, 0.95) to avoid a particle trapping in the vortex, because particles started at (1,1) are trapped in the electroconvective vortex (data not shown). (b) Simulation and experimental results of the charged particle trajectory and (c) forces experienced by the particle as a function of ( ̅ , ̅ ).

Conclusions
In conclusion, we performed numerical simulations on ICP phenomena in the ICP-2C-EF/HF device under the application of coupled hydrodynamic flow and electric fields. Particle trajectories are stimulated by the interaction of redistributed hydrodynamic flow and electric fields and negatively charged particles with one micrometer of diameter, representing bacteria. The simulated charged particle trajectory under ICP by the charged IC channel surface was validated with experimental results. The generation of an electroconvective vortex at the entrance could be analyzed by its position and velocity, enabling precise tracking of charged particle trajectory under hydrodynamic and electric force fields. The numerical analysis is not only tightly matched with experimental results but Position and color of sphere indicate where particle travels for corresponding time and magnitude of particle velocity, respectively. Path and color of line indicate trajectory and velocity throughout its entire course. The particle trajectory simulation started at (1, 0.95) to avoid a particle trapping in the vortex, because particles started at (1,1) are trapped in the electroconvective vortex (data not shown). (b) Simulation and experimental results of the charged particle trajectory and (c) forces experienced by the particle as a function of U in , ψ in .

Conclusions
In conclusion, we performed numerical simulations on ICP phenomena in the ICP-2C-EF/HF device under the application of coupled hydrodynamic flow and electric fields. Particle trajectories are stimulated by the interaction of redistributed hydrodynamic flow and electric fields and negatively charged particles with one micrometer of diameter, representing bacteria. The simulated charged particle trajectory under ICP by the charged IC channel surface was validated with experimental results. The generation of an electroconvective vortex at the entrance could be analyzed by its position and velocity, enabling precise tracking of charged particle trajectory under hydrodynamic and electric force fields. The numerical analysis is not only tightly matched with experimental results but also able to provide physical insight into ICP phenomena under parallel force fields to simulate experimental conditions for various applications, such as biomarker concentration and separation to improve sensing reliability and detection limits in analytic chemistry. However, the change in the particle movement according to various particle properties, such as size, charge, and density, within the same force field should be further studied.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/mi12080903/s1, Figure S1: The schematic illustration of the experimental set-up of ICP-2C-EF/HF system. To avoid a fluctuation of a generated bubble by electrolysis, flow generated by hydrostatic pressure difference is implemented with two open reservoirs. Figure S2: Schematic illustration of an ion-permselective material coated channel model and boundary conditions for numerical simulation. Table S1: Parameters used for simulation.

Conflicts of Interest:
There are no conflicts of interest to declare.