Electroosmotic Flow Hysteresis for Fluids with Dissimilar pH and Ionic Species

Electroosmotic flow (EOF) involving displacement of multiple fluids is employed in micro-/nanofluidic applications. There are existing investigations on EOF hysteresis, i.e., flow direction-dependent behavior. However, none so far have studied the solution pair system of dissimilar ionic species with substantial pH difference. They exhibit complicated hysteretic phenomena. In this study, we investigate the EOF of sodium bicarbonate (NaHCO3, alkaline) and sodium chloride (NaCl, slightly acidic) solution pair via current monitoring technique. A developed slip velocity model with a modified wall condition is implemented with finite element simulations. Quantitative agreements between experimental and simulation results are obtained. Concentration evolutions of NaHCO3–NaCl follow the dissimilar anion species system. When NaCl displaces NaHCO3, EOF reduces due to the displacement of NaHCO3 with high pH (high absolute zeta potential). Consequently, NaCl is not fully displaced into the microchannel. When NaHCO3 displaces NaCl, NaHCO3 cannot displace into the microchannel as NaCl with low pH (low absolute zeta potential) produces slow EOF. These behaviors are independent of the applied electric field. However, complete displacement tends to be achieved by lowering the NaCl concentration, i.e., increasing its zeta potential. In contrast, the NaHCO3 concentration has little impact on the displacement process. These findings enhance the understanding of EOF involving solutions with dissimilar pH and ion species.


Introduction
An applied electric field induces electroosmotic flow (EOF) of fluids in a micro-/nanochannel. Upon contact with an aqueous solution, the negative charges developed on the channel wall result in the repulsion and attraction of the negative and positive ions in the electrolyte solution. Electrical double layer (EDL) of nanometer thickness is formed as a result of the net positively charged layer. When an electric field is applied, an electrical body force is experienced by the EDL, which drives its motion along the direction of the electric field. Through viscous drag, the bulk fluid acquires the flow momentum from the EDL. This generates a plug-like fluid flow. With comparatively thin EDL thickness to channel dimensions, the Helmholtz-Smoluchowski equation can be used to determine the EOF velocity: where ε r and ε o are the fluid relative permittivity and the free space permittivity, respectively, µ is the fluid viscosity, E is the electric field applied externally, and ζ is the electrostatic wall zeta potential. The EDL thickness (for a symmetric electrolyte) can be characterized by the Debye length: where k b is the Boltzmann constant, N a is the Avogadro constant, e is the electron charge, T is the temperature, c o is the electrolyte solution concentration, and z is the ion species charge number. Electroosmotic flow (EOF) involving displacement of multiple fluids is widely employed in micro-/nanofluidic applications. These include electrokinetically controlled DNA hybridization [1][2][3], analyte mixing [4,5], fluid pumping [6][7][8][9], chemical species/particle separation [10], and pre-concentration [11,12]. In micro-energy systems, liquid flooding is a typical problem in the Proton Exchange Membrane (PEM) fuel cell, which reduces the fuel cell power greatly. With the integration of planar electroosmotic pump (EOP), the excess liquid is actively removed through microflow channels [7]. Mixing of analytes can be accomplished with the electrokinetic instability (EKI) approach by rapid folding and stretching of fluids through proper designs of micro-mixing devices [4].
In numerous practical applications, multiple fluids with dissimilar properties, e.g., concentration and conductivity, are commonly involved. For instance, pre-concentration methods such as isotachophoresis (ITP) and field amplified sample stacking (FASS) [11,12] involve micro-mixing of fluids with dissimilar conductivities. With a constant current setting, ITP forms zones of concentrated/focused sample ions in order of the ionic mobility [13][14][15]. In contrast, FASS employs an applied electric field with a constant voltage to stack/accumulate the sample ions at the interface between the low conductivity sample solution and the high conductivity background electrolyte (BGE) [16,17]. Sample dispersion due to the non-uniform EOF velocities adversely affects the resolution and sensitivity of both ITP and FASS; this can be minimized by coating water-soluble polymers [18,19], and constructing nanostructures within a microchannel [20][21][22].
Many investigations have been conducted on EOF of fluids with different conductivities/concentrations [23][24][25][26]. Tang et al. [23] performed an investigation on electroosmotic displacement flow with two or three fluids of different conductivity ratios, with the velocity of the fluid interface and the current-time response theoretically derived. Ren et al. [24,25] conducted two-fluid displacement flow experiments for various solution concentrations. They related the flow processes with the electrical equivalent circuit model, whereby the fluid resistances were determined by their electrical conductivities. Mampallil et al. [26] adopted two-fluid electroosmotic displacement flow experiments to measure the surface charge of microchannel. The current-time response of the displacement process was curvefitted to a theoretical expression to acquire the surface charge measurement. In addition, there are also numerous applications that involve EOF of several fluids, such as ionic transistor and diode [27,28], and transdermal drug delivery [29]. These examples [23][24][25][26][27][28][29] have demonstrated the importance of EOF on multiple fluid systems in practical applications.
However, the EOF of two dissimilar fluids demonstrates the hysteresis effect experimentally, known as EOF hysteresis, as coined by our group [30][31][32][33], whereby the EOF velocity/flow rate for fluid 1 displacing fluid 2 is different from fluid 2 displacing fluid 1. The hysteretic behavior is exhibited as the EOF velocity/flow rate and the final content in a microchannel are influenced by the interplay between the electric field distribution, induced EOF flow and ion migration; these are functions of the direction of displacement flow. Understanding the effect of EOF hysteresis is essential as its neglection may result in experimental data to be interpreted inappropriately.
Our group have numerically investigated EOF hysteresis involving fluids with different concentrations [30]. Rather counter-intuitively, it was discovered that the hysteresis effect arises from the depletion and accumulation of the minority pH-governing ions caused by the electromigrative flux imbalance; the pH, zeta potential, and EOF velocity are altered differently according to the displacement flow directions. The mechanics have been affirmed by the investigation of Li et al. [34]. By including the proton transport, carbonate equilibrium, and surface complexation reaction equations, their semi-analytical model can capture the experimentally observed hysteretic effect. Further study had been performed by our group to validate the hypothesis [31]. With the application of pH-sensitive dye, the pH changes in the microchannel during EOF for different displacement flow directions were experimentally quantified.
We had also investigated EOF hysteresis for fluids with dissimilar ionic species [32,33]. Two mechanics were discovered to induce EOF hysteresis for fluids with dissimilar cation species [32], namely the interfacial sharpening/widening effect caused by the difference in solution conductivities, and the ion concentration evolutions, that lead to the variation of EOF velocities/flow rates for different flow directions. Contrary to conventional electroosmotic displacement flows, the EOF of fluids with dissimilar anion species [33] was found to illustrate complex behavior, i.e., the displacing fluids cannot fully displace the residing fluids. The observed EOF hysteresis is caused by the changing ion concentrations due to the upstream migrating anions, and the diffusive-interface-induced concentration changes at the inlet of microchannel.
Despite previous studies [30][31][32][33], no investigation has heretofore been carried out on EOF hysteresis for a solution pair system comprising of dissimilar ionic species and substantial pH difference. Their rather different characteristics tend to result in a complex hysteretic behavior, thus rendering a more complex analysis. Our previous simulation models [30][31][32][33] will require excessive computational effort for this investigation, and thus may not be able to accurately reflect the displacement flow conditions due to computational constraint. Therefore, a slip-velocity model with modified wall condition is developed, which could describe the phenomena exhibited with sufficient precision and accuracy as compared to our previous models [30][31][32][33], but with much less computational effort.
A comprehensive understanding of their hysteresis phenomena is necessary to provide accurate EOF interpretation for analytical systems [35][36][37] involving fluids with dissimilar pH and ionic species. The findings of this investigation will aid in the precise manipulation of EOF conditions for effective improvement in micro-/nanofluidic applications with inhomogeneous solutions.

Experimental Parameters and Conditions
Investigation of EOF hysteresis for a solution pair system comprising of dissimilar ionic species and substantial pH difference was conducted. Sodium bicarbonate (NaHCO 3 ) and sodium chloride (NaCl) solution pair was selected as a model for the examination of the hysteresis behavior. NaHCO 3 alkaline salt breaks down to form sodium ion Na + and bicarbonate ion HCO 3 − in water, resulting in an alkaline solution (high pH). While NaCl neutral salt breaks down to form Na + and chloride ion Cl − in water, which gives a slightly acidic solution (low pH) due to the dissolved carbon dioxide from the atmosphere [30,31]. The characteristics of the solution pair complicate the electroosmotic displacement flow process. The solution concentrations and applied electric field were varied to comprehend their effects on the hysteresis phenomena. The experimental conditions and parameters investigated are shown in Table 1. Table 1. Experimental conditions with parameter variations for electroosmotic displacement flow of sodium bicarbonate (NaHCO 3 ) and sodium chloride (NaCl) solution pair in two flow directions, i.e., NaCl → NaHCO 3 (arrow indicates EOF direction), and NaHCO 3 → NaCl.

Parameters Investigated Experimental Conditions with Parameter Variations (Two Flow Directions)
Concentration of solution pair 0.5 mM NaHCO 3 -0.5 mM NaCl Electric field = 125 V·cm −1 A numerical investigation was conducted to study the effects of varying the concentration of the solution pair and applied electric field (see Section 3). Knowledge gained from understanding these systems with reference to the basic fluid configuration (i.e., same concentration for NaHCO 3 and NaCl solutions, see Table 1) was employed in the exploration of the underlying mechanics that influence the displacement process, when the concentration of each individual solution was varied.

Experimental Setup and Materials
Current monitoring technique [30][31][32][33] was employed for the observation of EOF behavior during displacement flow process (see Section 2.3), and the measurement of the zeta potential (see Section 2.4). An illustration of the experiment setup is depicted in Figure 1. A high voltage power supply (CZE1000R, Spellman, Hauppauge, NY, USA) was used to supply the electric field for inducing EOF. To monitor the current change across the microchannel, a picoammeter (Keithley 6485, Tektronix, Singapore) was connected in series. The two devices were controlled by a customized Labview program to record the readings of current and voltage with a data acquisition card (PCI-6052E, National Instrument, Austin, TX, USA).

Experimental Setup and Materials
Current monitoring technique [30][31][32][33] was employed for the observation of EOF behavior during displacement flow process (see Section 2.3), and the measurement of the zeta potential (see Section 2.4). An illustration of the experiment setup is depicted in Figure 1. A high voltage power supply (CZE1000R, Spellman, Hauppauge, NY, USA) was used to supply the electric field for inducing EOF. To monitor the current change across the microchannel, a picoammeter (Keithley 6485, Tektronix, Singapore) was connected in series. The two devices were controlled by a customized Labview program to record the readings of current and voltage with a data acquisition card (PCI-6052E, National Instrument, Austin, TX, USA). By dissolving the NaHCO3 and NaCl salts (Sigma-Aldrich, Saint Louis, MO, USA) in deionized (DI) water, 0.01M stock solutions for NaHCO3 and NaCl were prepared. The prepared solutions were further diluted to the concentrations required for experimentation (see Section 2.1). Measurements of the solution properties, i.e., pH and conductivity, By dissolving the NaHCO 3 and NaCl salts (Sigma-Aldrich, Saint Louis, MO, USA) in deionized (DI) water, 0.01 M stock solutions for NaHCO 3 and NaCl were prepared. The prepared solutions were further diluted to the concentrations required for experimentation (see Section 2.1). Measurements of the solution properties, i.e., pH and conductivity, were carried out with a pH meter (Mettler Toledo, FiveEasy plus, Singapore) and conductivity meter (IONCheck 65, Radiometer Analytical, Loveland, CO, USA). Table 2 presents the measured pH values and conductivities. Microchannels, i.e., flexible fused silica microcapillaries (polyimide coated, Polymicro Technologies, Phoenix, AZ, USA), with 100 µm nominal internal diameter were cut into 8 cm length with a Shortix Column Cutter (SGT Ltd., Singapore). Acetone was used to flush the microcapillaries, followed by DI water. Lastly, the capillaries were filled with the electrolyte solutions. The two ends of the microchannel were connected to two reservoirs made of Telfon having both depth and diameter of 2 cm (see Figure 1). The zeta potential developed at the glass/silica surface is caused by the deprononation (proton removal) of silanol (SiOH) groups in contact with an aqueous solution: SiOH + H 2 O SiO − + H 3 O + [30]. For a decrease in pH, i.e., an increase in H 3 O + concentration, the equilibrium will shift to the left, resulting in less SiO − groups; this lower concentration of the negatively-charged groups on the glass surface will reduce the absolute zeta potential value. The reverse is true for an increase in pH.
Upon contact with the NaCl or NaHCO 3 electrolyte solutions, the negative charges developed on the channel wall result in the repulsion of the Cl − or HCO 3 − (negative) and attraction of the Na + (positive) ions. This forms a net positively charged EDL layer. When an electric field is applied, an electrical body force is experienced by the EDL, which drives its motion along the direction of the electric field. Through viscous drag, the bulk fluid acquires the flow momentum from the EDL, which generates a plug-like fluid flow. The ion mobilities of Na + , Cl − and HCO 3 − are 5.194 × 10 −8 m 2 ·V −1 ·s −1 , −7.919 × 10 −8 m 2 ·V −1 ·s −1 and −4.303 × 10 −8 m 2 ·V −1 ·s −1 , respectively. Based on the ion mobilities, the ion transport numbers, defined as the fractions of the total electrical current carried by the ion species, can be calculated. For NaCl solution, the ion transport numbers of Na + and Cl − are 0.396 and 0.604, respectively; for NaHCO 3 solution, the ion transport numbers of Na + and HCO 3 − are 0.547 and 0.453, respectively. As the size of the reservoirs was sufficiently large, the change in liquid level during the experiment was negligible. This minimized the back pressure from the liquid level difference in the reservoirs [38]. With small reservoirs, EOF might be affected by the local pH change induced by electrolysis at the electrodes due to the generation of H + and OH − ions [39]. The use of large-volume reservoirs significantly diluted the H + and OH − concentrations. The electrodes were also placed distant from the microchannel inlet and outlet (approximately 2 cm) to avoid unnecessary change in pH in the channel [40].

EOF of NaHCO 3 -NaCl Solution Pair System
Electroosmotic displacement flow experiments for the NaHCO 3 -NaCl solution pair system were performed via current monitoring technique (see Figure 1), according to the experimental parameters and conditions in Table 1. The experiments were conducted in two displacement flow directions, i.e., NaCl → NaHCO 3 (the arrow indicates EOF direction), and NaHCO 3 → NaCl. The displacing solution was placed in the anode reservoir for displacing the original residing solution in the microchannel; the channel and cathode reservoir were filled with the residing solution. For the revelation of the hysteresis behavior, comparisons of the current-time responses in two opposite flow directions were made.
Due to the low conductivities and concentrations of the electrolytes in this investigation, the effect of joule heating is negligible [41]. By considering the balance in energy between the energy storage ∆E st and energy generation E g in the fluid, the effect of joule heating can be conservatively estimated with the methodology detailed in the work of Arulanandam and Li [42]. The ion concentrations will be changing during the experiments. Hence, the worst scenario was having NaHCO 3 concentration = 5 mM, microchannel cross-sectional area = 7850 µm 2 , applied electric field = 125 V·cm −1 and experimental duration = 600 s; the estimated temperature increase is only 0.8 • C and is insignificant to affect the electroosmotic displacement flow process.

Zeta Potential Measurements
Zeta potential is a crucial parameter which determines the EOF velocity for electroosmotic displacement flow in a microchannel. As numerical simulations are required for the evaluation of the displacement flow behaviors and ion concentration distributions, it will be necessary for the prescription of exact wall zeta potential conditions to obtain accurate EOF simulations (see Section 3). Therefore, zeta potential measurements of the solutions required for simulations will be carried out.
Current monitoring technique (see Figure 1) was employed for the measurement of zeta potential. The cathode reservoir and microchannel were filled with the solution for measurement (residing solution), while the anode reservoir was filled with 95% concentration of the measurement solution (displacing solution), i.e., 5% difference in concentration. Electric field of 125 V·cm −1 was applied to generate EOF across the two reservoirs. To ensure reliability and consistency of the results, for each data point, the experiments were performed five times.
Through dividing the microchannel length with the displacement time, i.e., time for the current to attain the steady-state value of displacing solution (see Figure 2), the EOF velocity can be calculated with [30][31][32][33]: where L and T Displace represent the microchannel length and the displacement time, respectively.
reservoir were filled with the residing solution. For the revelation of the hysteresis behav-ior, comparisons of the current-time responses in two opposite flow directions were made. Due to the low conductivities and concentrations of the electrolytes in this investigation, the effect of joule heating is negligible [41]. By considering the balance in energy between the energy storage ∆Est and energy generation Eg in the fluid, the effect of joule heating can be conservatively estimated with the methodology detailed in the work of Arulanandam and Li [42]. The ion concentrations will be changing during the experiments. Hence, the worst scenario was having NaHCO3 concentration = 5 mM, microchannel cross-sectional area = 7850 µm 2 , applied electric field = 125 V·cm −1 and experimental duration = 600 s; the estimated temperature increase is only 0.8 °C and is insignificant to affect the electroosmotic displacement flow process.

Zeta Potential Measurements
Zeta potential is a crucial parameter which determines the EOF velocity for electroosmotic displacement flow in a microchannel. As numerical simulations are required for the evaluation of the displacement flow behaviors and ion concentration distributions, it will be necessary for the prescription of exact wall zeta potential conditions to obtain accurate EOF simulations (see Section 3). Therefore, zeta potential measurements of the solutions required for simulations will be carried out.
Current monitoring technique (see Figure 1) was employed for the measurement of zeta potential. The cathode reservoir and microchannel were filled with the solution for measurement (residing solution), while the anode reservoir was filled with 95% concentration of the measurement solution (displacing solution), i.e., 5% difference in concentration. Electric field of 125 V·cm −1 was applied to generate EOF across the two reservoirs. To ensure reliability and consistency of the results, for each data point, the experiments were performed five times.
Through dividing the microchannel length with the displacement time, i.e., time for the current to attain the steady-state value of displacing solution (see Figure 2), the EOF velocity can be calculated with [30][31][32][33]: where L and TDisplace represent the microchannel length and the displacement time, respectively.  Thereafter, by substituting the measured EOF velocity from Equation (3) into the Helmholtz-Smoluchowski slip velocity equation (see Equation (1)), the zeta potential may be written as [30][31][32][33]: where µ is the fluid viscosity, ε r and ε o are the fluid relative permittivity and the free space permittivity, respectively, and E is the applied electric field. The measured zeta potential values of the solutions required for the numerical simulations (see Section 3) are shown in Table 3. Table 3 indicates that by lowering the NaCl concentration from 1.24 mM to 0.5 mM, i.e., by approximately 60%, the absolute zeta potential value increased significantly by about 40%. By decreasing the NaCl concentration, there will be lesser Na + ions (positive ions) for the shielding of the negative wall charge (fused silica), which leads to the increase of absolute zeta potential value. While lowering the NaHCO 3 concentration from 1 mM to 0.375 mM, i.e., by approximately 60%, the absolute zeta potential value increased only by about 10%, i.e., much less. The absolute zeta potential value is supposed to similarly increase significantly when the NaHCO 3 concentration is decreased. However, decreasing the NaHCO 3 concentration causes the pH to decrease (see Table 2), which will cause a decrease in the absolute zeta potential value. The competing effect between the solution concentration and pH results in the slight increase of the absolute zeta potential value for the case of NaHCO 3 .
The Poisson-Nernst-Planck (PNP) model [32,33] considers the electric field variation (due to a difference in solution conductivity) and the transportation of main constituent ionic species during two-fluid displacement flow. With its ability to capture the transportation of the ionic species under diffusive, convective, and electromigrative effects, the variation of zeta potential is accounted for by prescribing a constant surface charge density at the microchannel wall. Surface charge regulation model [30,49] provides a better elucidation of the two-fluid displacement flow process. In addition to the features of PNP model, the surface charge regulation model includes the transportation of minority pH-governing ionic species with its associated reversible acid-based chemical reactions. The model can simulate the variation of zeta potential due to the concentration changes of main constituent ionic species, and the pH changes attributed by the accumulation/depletion of minority pHgoverning ionic species.
However, excessive computational effort will be required to implement the surface charge regulation model. Our investigation employs the NaHCO 3 -NaCl solution pair system, where dissimilar ionic species and substantial pH difference are involved, thus rendering a more complex analysis. However, the PNP model is not applicable for our study. This is because the large pH difference between the NaHCO 3 and NaCl solutions (see Table 2) yields a large surface charge difference; as such, the specification of an average constant charge density will not accurately reflect the displacement flow conditions. Therefore, a slip-velocity model with modified wall condition, which could describe the phenomena exhibited, will be adopted. This model, which has the same capabilities as the PNP model, considers the electric field variation (due to the difference in solution conductivity) and transportation of main constituent ionic species during two-fluid displacement flow. By prescribing a varying wall zeta potential boundary according to the flow condition, it can execute simulation runs with sufficient precision and accuracy as compared to the surface charge regulation model, but with much lesser computational effort.

Governing Equations
The modified slip-velocity model is developed based on several fundamental equations [30][31][32][33]. The applied electric potential ϕ is governed by the Laplace equation, which accounts for the electric field variation: where σ = F∑z i u m(i) c i is the solution conductivity, c i , z i , and u m(i) are the ion concentration, ion charge number, and ionic mobility, respectively, and F is the Faraday constant. Governed by the gradients of diffusive, convective, and electromigrative fluxes respectively, the change of ionic species concentration with time is described by the Nernst-Planck (NP) equation: where v is the velocity of the fluid, and D i the diffusion coefficient of the ion. The flow of an incompressible Newtonian fluid is defined by the Navier-Stokes (NS) and continuity equations: where p is the pressure, ρ the density, and µ the fluid viscosity. As the Reynolds number is usually below 1 for a microchannel, the inertial term is to be neglected by assuming Stokes flow. The flow field is obtained by simultaneously solving the Laplace (Equation (5)), NP (Equation (6)), NS and continuity equations (Equations (7) and (8)). The values of the constants employed for the simulations can be found in Table 4. Table 4. Symbols and values of constants for simulations, with ion mobility determined by (z i D i F)/(RT).

Constants Symbol (Unit) Value
Fluid relative permittivity ε r 80 Free space permittivity

Simulation Domain
Finite element numerical simulations were performed using the software COMSOL Multiphysics. The simulated domain (see Figure 3) is a 10 µm diameter microchannel of 0.5 cm in length. Axisymmetry is assumed for the fluid flow about the center axis. The strong coupling between the applied electric potential, concentrations of ions and fluid flow velocity in the governing equations (see Section 3.2) will result in significant computational workload. To ease this workload, the domain diameter and length were reduced from 100 µm (experimentally) to 10 µm (numerically), and from 8 cm (experimentally) to 0.5 cm (numerically), respectively. For non-overlapping EDL, EOF will be independent of the microchannel diameter. Thus, the reduction of the microchannel diameter in the simulations will not affect the EOF displacement process. In addition, with the application of the same electric fields for both experiments and simulations, the reduction of the channel length in the simulations should not affect its accuracy in representing the flow behaviors of the experimental runs after normalization.

Simulation Domain
Finite element numerical simulations were performed using the software COMSOL Multiphysics. The simulated domain (see Figure 3) is a 10 µm diameter microchannel of 0.5 cm in length. Axisymmetry is assumed for the fluid flow about the center axis. The strong coupling between the applied electric potential, concentrations of ions and fluid flow velocity in the governing equations (see Section 3.2) will result in significant computational workload. To ease this workload, the domain diameter and length were reduced from 100 µm (experimentally) to 10 µm (numerically), and from 8 cm (experimentally) to 0.5 cm (numerically), respectively. For non-overlapping EDL, EOF will be independent of the microchannel diameter. Thus, the reduction of the microchannel diameter in the simulations will not affect the EOF displacement process. In addition, with the application of the same electric fields for both experiments and simulations, the reduction of the channel length in the simulations should not affect its accuracy in representing the flow behaviors of the experimental runs after normalization. The simulated domain consisted of 24,000 quadrilateral elements (1000 and 24 edge elements, respectively, in the axial and radial directions). The fluid flow velocity and pressure were discretized with linear elements. While the ion concentrations and applied electric potential were discretized with quadratic elements. Through steady-state simulation, a convergence test was carried out with a higher element number, i.e., 30,000 elements (25% increment). The results were found to have a percentage change less than 1%, which is negligible. Between subsequent iterations, a relative tolerance of less than 0.01 (i.e., approximately 1%) was used as the convergence criterion.

Boundary and Initial Conditions
The steady-state solution for EOF of a single fluid is first obtained; it is used as the initial condition for solving the time-dependent solution of two-fluid electroosmotic displacement flow. The boundary conditions for the steady-state single fluid EOF numerical model are listed in Table 5. The simulated domain consisted of 24,000 quadrilateral elements (1000 and 24 edge elements, respectively, in the axial and radial directions). The fluid flow velocity and pressure were discretized with linear elements. While the ion concentrations and applied electric potential were discretized with quadratic elements. Through steady-state simulation, a convergence test was carried out with a higher element number, i.e., 30,000 elements (25% increment). The results were found to have a percentage change less than 1%, which is negligible. Between subsequent iterations, a relative tolerance of less than 0.01 (i.e., approximately 1%) was used as the convergence criterion.

Boundary and Initial Conditions
The steady-state solution for EOF of a single fluid is first obtained; it is used as the initial condition for solving the time-dependent solution of two-fluid electroosmotic displacement flow. The boundary conditions for the steady-state single fluid EOF numerical model are listed in Table 5.  The effects of varying the concentration of the NaHCO 3 -NaCl solution pair and applied electric field will be numerically studied. Ion concentrations of the simulation domain, inlet, and outlet (see Figure 3) were set according to the experimental conditions (see Table 1). Electroneutrality condition was enforced to satisfy charge neutrality. Voltages of 62.5 V and 93.75 V were prescribed at the inlet (0 V at outlet) to establish electric fields of 125 V·cm −1 and 187.5 V·cm −1 , respectively. Electrical insulation condition was set at the microchannel wall to limit current flow within the channel.
EOF of solutions with dissimilar anion species exhibits concentration evolutions [33]. The NaHCO 3 -NaCl solution pair consists of dissimilar anions Cl − and HCO 3 − , and common cation Na + . Hence, concentration evolutions occur for displacement flow of the NaHCO 3 -NaCl solution pair, which follow the dissimilar anion species system. The evolved concentrations were determined from simulations, which can be found in Table 6. Knowing the evolved concentrations ensures that precise wall boundary conditions can be prescribed (see Table 7) for accurate EOF simulations.    Table 3).
Varying zeta potential functions (see Table 7) were specified at the wall boundary to capture the changes of zeta potential for different displacement flow directions. A formulation of these functions were based on the experimentally measured zeta potentials of the displacing, residing, and evolved concentration solutions for the various experimental flow conditions (see Tables 1, 3 and 6). Based on these experimentally measured zeta potentials, these functions can accurately reflect the changes of zeta potential according to the content in the microchannel depending on the displacement flow conditions. The inlet and outlet pressures were prescribed to be zero.
The time-dependent electroosmotic displacement flow was solved by modifying the boundary conditions at the inlet. With the steady-state solution for EOF of a single fluid set as the initial condition, ion concentrations at the inlet were ramped from the residing solution to the displacing solution (see Table 1) in an arbitrarily short time of 0.01 s for commencing the displacement process.

Results and Discussion
EOF behavior of the 0.5 mM NaHCO 3 -NaCl solution pair at the electric field of 125 V·cm −1 is shown in Figure 4. Normalization of the currents and times enables the comparison between the experimental and numerical simulation results. The normalized currents and times are calculated with Equations (9) and (10), respectively: where the initial currents of solutions are represented by I 0 with NaHCO 3 and NaCl subscripts, respectively, and the time to reach steady-state current for NaCl → NaHCO 3 (arrow indicates EOF direction) is represented by T SS NaCl → NaHCO3 . and outlet pressures were prescribed to be zero. The time-dependent electroosmotic displacement flow was solved by modifying the boundary conditions at the inlet. With the steady-state solution for EOF of a single fluid set as the initial condition, ion concentrations at the inlet were ramped from the residing solution to the displacing solution (see Table 1) in an arbitrarily short time of 0.01 s for commencing the displacement process.

Results and Discussion
EOF behavior of the 0.5 mM NaHCO3-NaCl solution pair at the electric field of 125 V·cm −1 is shown in Figure 4. Normalization of the currents and times enables the comparison between the experimental and numerical simulation results. The normalized currents and times are calculated with Equations (9) and (10), respectively: T* = T/T SS NaCl → NaHCO3 (10) where the initial currents of solutions are represented by I 0 with NaHCO3 and NaCl subscripts, respectively, and the time to reach steady-state current for NaCl → NaHCO3 (arrow indicates EOF direction) is represented by T SS NaCl → NaHCO3.  ) and T* = T/T SS NaCl → NaHCO3 respectively, where initial currents of solutions are represented by I 0 with NaHCO 3 and NaCl subscripts respectively, and time to reach steady-state current for NaCl → NaHCO 3 (arrow indicates EOF direction) is represented by T SS NaCl → NaHCO3 .
Our slip-velocity EOF model with modified wall boundary condition for zeta potential variation (see Section 3) predicts the experimental observations excellently (see Figure 4). With the ability to facilitate accurate zeta potential variation to accommodate the displacement and evolution of the solution concentration, our numerical model can capture the essence of the experimental phenomena, and its validity is affirmed.
From the normalized current-time curves (see Figure 4), two phases were observed (namely Phases 1 and 2) to be separated by an abrupt gradient change, which typically occurs for the EOF of solutions with dissimilar ionic species [32,33]. For NaCl → NaHCO 3 (the arrow indicates the EOF direction), a decrease and an increase in the current were observed for Phases 1 and 2, respectively (see Figure 4). In contrast, for 0.5 mM NaHCO 3 → 0.5 mM NaCl, an increase and a decrease in current were observed for Phases 1 and 2, respectively (see Figure 4). For both cases, the current-time curves were unable to attain the steady-state values of the displacing solutions at the end of Phase 2, indicating that the displacing solutions were unable to displace completely the residing solutions. Figure 5 shows the simulated evolution of ion concentration distributions for 0.5 mM NaCl → 0.5 mM NaHCO 3 . At normalized time T* = 0.136 (Phase 1), depletion of Na + and HCO 3 − ions occurs, with NaHCO 3 concentration reducing to 0.375 mM. The change in concentration can be predicted by the Kohlrausch regulating function (KRF) [33,51], which satisfies the current continuity and electroneutrality: where z i is the ion charge number, c i the ion concentration as a function of microchannel axial coordinate X and time T, and u m(i) the ion mobility.
observed for Phases 1 and 2, respectively (see Figure 4). In contrast, for 0.5 mM NaHCO3 → 0.5 mM NaCl, an increase and a decrease in current were observed for Phases 1 and 2, respectively (see Figure 4). For both cases, the current-time curves were unable to attain the steady-state values of the displacing solutions at the end of Phase 2, indicating that the displacing solutions were unable to displace completely the residing solutions. Figure 5 shows the simulated evolution of ion concentration distributions for 0.5 mM NaCl → 0.5 mM NaHCO3. At normalized time T* = 0.136 (Phase 1), depletion of Na + and HCO3 − ions occurs, with NaHCO3 concentration reducing to 0.375 mM. The change in concentration can be predicted by the Kohlrausch regulating function (KRF) [33,51], which satisfies the current continuity and electroneutrality:

KRF(X) = ∑[zici(X,T)/um(i)]
where zi is the ion charge number, ci the ion concentration as a function of microchannel axial coordinate X and time T, and um(i) the ion mobility. Normalized time T* = T/T SS NaCl → NaHCO3 , where T SS NaCl → NaHCO3 is time to reach steady-state current for NaCl → NaHCO 3 . Normalized X* = X/L, where X is axial coordinate and L microchannel length.
An amount of 0.5 mM NaCl has a KRF value of 1.59 × 10 7 mol·V·s·m −5 , and 0.5 mM NaHCO 3 has a KRF value of 2.12 × 10 7 mol·V·s·m −5 . Anions migrate upstream in opposition to the EOF flow. When 0.5 mM NaCl → 0.5 mM NaHCO 3 , the migration of HCO 3 − causes a change in NaHCO 3 concentration (see Figure 5b). The NaHCO 3 concentration is lowered from 0.5 mM to 0.37 5mM to match the KRF of 0.5 mM NaCl (solution concentration derived from known KRF of 0.5 mM NaCl and ion mobilities of Na + and HCO 3 − from Table 4, see Equation (11)), with a corresponding drop in conductivity (see Table 2). Therefore, a decrease in current was observed during Phase 1 for 0.5 mM NaCl → 0.5 mM NaHCO 3 (see Figure 4). When the residing solution 0.5 mM NaHCO 3 was flushed out of the microchannel (at T* = 0.227 (Phase 1), see Figure 5c), it was indicated by the abrupt gradient change from the current-time curve (see Figure 4).
In Phase 2 of 0.5 mM NaCl → 0.5 mM NaHCO 3 , displacement of the reduced concentration 0.375 mM NaHCO 3 by 0.5 mM NaCl continues (see Figure 5d). Since 0.5 mM NaCl has higher conductivity than 0.375 mM NaHCO 3 (see Table 2), an increase in current was detected for Phase 2 (see Figure 4). However, 0.375 mM NaHCO 3 cannot be completely flushed out, and the normalized X* is unchanged at 0.75 when T* = 1 (see Figure 5d), where X* = X/L with X the axial coordinate and L the microchannel length. The gradient of EOF induced a convective flux (term on the right of Equation (6)) experienced by the Cl − anions, which gradually reduces as the flow propagates downstream. This is due to the displacement of NaHCO 3 out of the microchannel; due to its high pH, the NaHCO 3 solution has high absolute zeta potential value (see Tables 2 and 3) that are supposed to facilitate faster EOF. Thus, the gradient of electromigrative flux (third term left of Equation (6)) cancels out that of the convective flux for Cl − anions. This resulted in the interface staying unchanged at X* = 0.75 when T* = 1 (see Figure 5d), and the current stabilized thereafter (see Figure 4). According to the discussion in our previous study [32], for 0.5 mM NaCl → 0.5 mM NaHCO 3 , as a result of the low conductivity of NaHCO 3 (see Table 2), the residing ion-depletion region (0.375 mM NaHCO 3 ) has a higher electric field than 0.5 mM NaCl (displacing electrolyte). Thus, as HCO 3 − ions diffuse to the boundary of the displacing electrolyte, they decelerate due to the low electric field in the displacing electrolyte region. This results in a sharp and constant interfacial region (see Figure 5).
The simulated ion concentration distributions when 0.5 mM NaHCO 3 → 0.5 mM NaCl are shown in Figure 6. The 0.5 mM NaHCO 3 is unable to flow out from the inlet reservoir to the microchannel for the entire displacement process, i.e., the interface stays at X* = 0. The HCO 3 − anions experience much stronger gradient of electromigrative flux than the EOF induced convective flux. This is because the initial residing solution NaCl has a low absolute zeta potential value (due to its low pH value) (see Tables 2 and 3) that generates slower EOF. infinite source for the HCO3 − ions, as well as an infinite sink for the incoming Cl − (migrating upstream), a diffusive mixture region developed near the vicinity of the inlet with three ion types, i.e., Na + , HCO3 − , and Cl − (see Figure 6), which fulfill both current continuity and electroneutrality. This diffusive-interface-induced concentration evolution has been discussed in detail in our previous investigation on EOF with dissimilar anion species [33]. The equilibrium concentration after the second evolution in the microchannel can be determined with [33]: where I SS (Exp) is the steady-state experimental current, EL the applied electric field along the microchannel, F the Faraday constant, AL the channel cross-sectional area, and ni the molecular formula. Normalized time T* = T/T SS NaCl → NaHCO3, where T SS NaCl → NaHCO3 is time to reach steady-state current for NaCl → NaHCO3. Normalized X* = X/L, where X is axial coordinate and L microchannel length.
The numerical simulation of 0.5 mM NaHCO3 → 0.5 mM NaCl demonstrates that the NaCl concentration decreases from 0.665 mM to 0.619 mM at the end of Phase 2 when T* = 1 (see Figure 6d). Electroosmotic displacement flow of 0.619 mM NaCl and 0.665 mM Normalized time T* = T/T SS NaCl → NaHCO3 , where T SS NaCl → NaHCO3 is time to reach steady-state current for NaCl → NaHCO 3 . Normalized X* = X/L, where X is axial coordinate and L microchannel length.
At T* = 0.227 (Phase 1 of 0.5 mM NaHCO 3 → 0.5 mM NaCl), accumulation of Na + and Cl − ions occurs, and the increase of NaCl concentration takes place at 0.665 mM (see Figure 6b). The KRF value of 0.5 mM NaHCO 3 is 2.12 × 10 7 mol·V·s·m −5 , and the KRF value of 0.5 mM NaCl is 1.59 × 10 7 mol·V·s·m −5 . To match the KRF of 0.5 mM NaHCO 3 , the upstream migration of Cl − causes the NaCl concentration to increase from 0.5 mM to 0.665 mM (solution concentration derived from known KRF of 0.5 mM NaHCO 3 and ion mobilities of Na + and Cl − from Table 4, see Equation (11)), which leads to an increase in conductivity. Hence, an increase in current was observed for 0.5 mM NaHCO 3 → 0.5 mM NaCl in Phase 1 (see Figure 4). The residing solution 0.5 mM NaCl was flushed out of the microchannel (at T* = 0.455 (Phase 1), see Figure 6c); this was indicated by the abrupt gradient change from the current-time curve (see Figure 4).
For 0.5 mM NaHCO 3 → 0.5 mM NaCl in Phase 2, a second concentration evolution happens (see Figure 6d). Initially, the inlet reservoir and the microchannel were filled with 0.5 mM NaHCO 3 and 0.5 mM NaCl, respectively. With the inlet reservoir acting as an infinite source for the HCO 3 − ions, as well as an infinite sink for the incoming Cl − (migrating upstream), a diffusive mixture region developed near the vicinity of the inlet with three ion types, i.e., Na + , HCO 3 − , and Cl − (see Figure 6), which fulfill both current continuity and electroneutrality. This diffusive-interface-induced concentration evolution has been discussed in detail in our previous investigation on EOF with dissimilar anion species [33]. The equilibrium concentration after the second evolution in the microchannel can be determined with [33]: where I SS (Exp) is the steady-state experimental current, E L the applied electric field along the microchannel, F the Faraday constant, A L the channel cross-sectional area, and n i the molecular formula.
The numerical simulation of 0.5 mM NaHCO 3 → 0.5 mM NaCl demonstrates that the NaCl concentration decreases from 0.665 mM to 0.619 mM at the end of Phase 2 when T* = 1 (see Figure 6d). Electroosmotic displacement flow of 0.619 mM NaCl and 0.665 mM NaCl progresses during Phase 2. As such, a slight current decrease was observed, and the current stabilized after the displacement process (see Figure 4). Through calculation with Equation (12) based on the experimental data, the NaCl equilibrium concentration is 0.640 ± 0.002 mM, and is reasonably similar to the simulated NaCl equilibrium concentration of 0.619 mM.
To maintain flow continuity over the entire microchannel during the electroosmotic displacement flow process, pseudo-pressure is generated due to non-uniform zeta potential. This resulted velocity profile resembles a combination of electroosmotic flow and pressure driven flow distributions, and deviates from the usual plug-like EOF profile. The simulated normalized velocity V* = (V − V Avg )/V Avg and pressure P along the normalized radial r* = r/R for normalized X* = 0.25, 0.5, and 0.75 of 0.5 mM NaCl-0.5 mM NaHCO 3 in two different flow directions are presented in Figures 7 and 8, respectively. The average velocity V Avg = Q/A L = Q/(πR 2 ), where Q is the flow rate obtained by integrating the radial velocity over the channel cross sectional area A L , and R is the channel radius.
For 0.5 mM NaCl → 0.5 mM NaHCO 3 , when normalized time T* = 0, V* and P for the different X* values are zero across r*, see Figure 7a. This is because of the uniform zeta potential of the residing NaHCO 3 at the initial state (see Figure 5a); plug-like EOF profile is observed, with the velocity having zero deviation (i.e., V* = 0) from the average velocity. As the displacement process progresses, as shown in Figure 5b-d, concentration evolutions occur with different fluid segments along the microchannel having different and thus non-uniform zeta potential; this generates different wall driving force. To maintain fluid flow continuity, internal fluid pressure is generated, which changes along the axial length. Negative pressure (back pressure) is generated due to EOF slowing down (reduction of velocity) as the displacement flow progresses. This fluid pressure variation causes the flow velocity profiles to evolve and deviate from the plug-like EOF profile, as shown in Figure 7b-d. For 0.5 mM NaHCO 3 → 0.5 mM NaCl, when T* = 0, V*, and P for the different X* values are zero across r* (see Figure 8a) due to the uniform zeta potential of the residing NaCl at the initial state (see Figure 6a). However, NaHCO 3 is unable to flow out from the inlet, and only slight concentration evolutions are shown in Figure 6b-d. Therefore, the pressure does not vary significant, and only small deviations of the flow velocity profile from the plug-like EOF profile are observed in Figure 8b-d. pressure driven flow distributions, and deviates from the usual plug-like EOF profile. The simulated normalized velocity V* = (V − VAvg)/VAvg and pressure P along the normalized radial r * = r/R for normalized X * = 0.25, 0.5, and 0.75 of 0.5 mM NaCl-0.5 mM NaHCO3 in two different flow directions are presented in Figures 7 and 8, respectively. The average velocity VAvg = Q/AL = Q/(πR 2 ), where Q is the flow rate obtained by integrating the radial velocity over the channel cross sectional area AL, and R is the channel radius. , where Q is flow rate obtained by integrating radial velocity over microchannel cross sectional area AL, and R is channel radius. Normalized r * = r/R. Normalized time T* = T/T SS NaCl → NaHCO3, where T SS NaCl → NaHCO3 is time to reach steady-state current for NaCl → NaHCO3. Normalized X*= X/L, where X is axial coordinate and L channel length.
For 0.5 mM NaCl → 0.5 mM NaHCO3, when normalized time T* = 0, V* and P for the different X* values are zero across r*, see Figure 7a. This is because of the uniform zeta potential of the residing NaHCO3 at the initial state (see Figure 5a); plug-like EOF profile is observed, with the velocity having zero deviation (i.e., V* = 0) from the average velocity. As the displacement process progresses, as shown in Figure 5b-d, concentration evolutions occur with different fluid segments along the microchannel having different and thus non-uniform zeta potential; this generates different wall driving force. To maintain fluid flow continuity, internal fluid pressure is generated, which changes along the axial length. Negative pressure (back pressure) is generated due to EOF slowing down (reduction of velocity) as the displacement flow progresses. This fluid pressure variation causes the flow velocity profiles to evolve and deviate from the plug-like EOF profile, as shown is time to reach steady-state current for NaCl → NaHCO 3 . Normalized X*= X/L, where X is axial coordinate and L channel length. different X* values are zero across r* (see Figure 8a) due to the uniform zeta potential of the residing NaCl at the initial state (see Figure 6a). However, NaHCO3 is unable to flow out from the inlet, and only slight concentration evolutions are shown in Figure 6b-d. Therefore, the pressure does not vary significant, and only small deviations of the flow velocity profile from the plug-like EOF profile are observed in Figure 8b-d. , where Q is flow rate obtained by integrating radial velocity over microchannel cross sectional area AL, and R is channel radius. Normalized r * = r/R. Normalized time T* = T/T SS NaCl → NaHCO3, where T SS NaCl → NaHCO3 is time to reach steady-state current for NaCl → NaHCO3. Normalized X*= X/L, where X is axial coordinate and L channel length.
The effect of varying the concentration of the solution pair was examined with 1 mM NaHCO3-NaCl solution pair (see Figure 9a). The electric field was kept the same as the case for the 0.5 mM NaHCO3-NaCl solution pair, in which the electric field strength was 125 V·cm −1 . The numerically predicted results match well with the experimental observations (see Figure 9a). The overall trend is almost identical to 0.5 mM NaHCO3-NaCl (see Figure 4) with an abrupt gradient change observed (separating Phases 1 and 2), and the current-time curve was unable to reach the current of the displacing solution. For 1 mM NaCl → 1 mM NaHCO3, the NaHCO3 concentration is reduced to 0.750 mM (according to Kohlrausch regulating function (KRF), see Equation (11)), which is two times that of 0.5 mM NaCl displacing 0.5 mM NaHCO3. Supposedly, the difference in concentration should affect the zeta potential and hence EOF strength. However, interestingly, the competing effect between the solution concentration and pH results in insignificant variation of the absolute zeta potential value, which can be seen from Tables 2 and 3. Hence, the displacement interface of 0.750 mM NaHCO3 stays unchanged at normalized X* = 0.75. is time to reach steady-state current for NaCl → NaHCO 3 . Normalized X*= X/L, where X is axial coordinate and L channel length.
The effect of varying the concentration of the solution pair was examined with 1 mM NaHCO 3 -NaCl solution pair (see Figure 9a). The electric field was kept the same as the case for the 0.5 mM NaHCO 3 -NaCl solution pair, in which the electric field strength was 125 V·cm −1 . The numerically predicted results match well with the experimental observations (see Figure 9a). The overall trend is almost identical to 0.5 mM NaHCO 3 -NaCl (see Figure 4) with an abrupt gradient change observed (separating Phases 1 and 2), and the current-time curve was unable to reach the current of the displacing solution. For 1 mM NaCl → 1 mM NaHCO 3 , the NaHCO 3 concentration is reduced to 0.750 mM (according to Kohlrausch regulating function (KRF), see Equation (11)), which is two times that of 0.5 mM NaCl displacing 0.5 mM NaHCO 3 . Supposedly, the difference in concentration should affect the zeta potential and hence EOF strength. However, interestingly, the competing effect between the solution concentration and pH results in insignificant variation of the absolute zeta potential value, which can be seen from Tables 2 and 3. Hence, the displacement interface of 0.750 mM NaHCO 3 stays unchanged at normalized X* = 0.75. For 1 mM NaHCO 3 → 1 mM NaCl, 1 mM NaHCO 3 cannot be displaced into the microchannel. The diffusive-interface-induced evolution near the inlet vicinity results in NaCl equilibrium concentration of 1.21 ± 0.012 mM (calculated based on experimental data with Equation (12)), which is approximately two times that of 0.5 mM NaHCO 3 displacing 0.5 mM NaCl. For 1 mM NaHCO3 → 1 mM NaCl, 1 mM NaHCO3 cannot be displaced into the microchannel. The diffusive-interface-induced evolution near the inlet vicinity results in NaCl equilibrium concentration of 1.21 ± 0.012 mM (calculated based on experimental data with Equation (12)), which is approximately two times that of 0.5 mM NaHCO3 displacing 0.5 mM NaCl. with NaHCO3 and NaCl subscripts respectively, and the time to reach steady-state current for NaCl → NaHCO3 (arrow indicates EOF direction) is represented by T SS NaCl → NaHCO3.
To study its effect on EOF with the 1 mM NaHCO3-NaCl solution pair, the electric field was increased by 50% to 187.5 V·cm −1 , as shown in Figure 9b. Good agreement is obtained between simulations and experimental results. The overall trend follows closely, and is similar to that of 1 mM NaHCO3-NaCl with an electric field of 125 V·cm −1 (see Figure 9a) and 0.5 mM NaHCO3-NaCl with an electric field of 125 V·cm −1 (see Figure 4). For 1 mM NaCl → 1 mM NaHCO3, despite increasing the electric field by 50% to 187.5 V·cm −1 , the NaHCO3 concentration is similarly reduced to 0.750 mM (not affected based on KRF, see Equation (11)) and without complete displacement by 1 mM NaCl. For 1 mM NaHCO3 → 1 mM NaCl with an electric field of 187.5 V·cm −1 , the NaCl equilibrium concentration is calculated to be 1.16 ± 0.050 mM (based on experimental data with Equation (12)), which is approximately the same as the case with the electric field of 125 V·cm −1 . Changing the electric field strength has no influence on the concentration evolution for the NaHCO3-NaCl solution pair system, except producing a much faster electroosmotic displacement process.
The effect of varying NaHCO3 concentration was investigated experimentally by employing 0.5 mM, 3 mM or 5 mM NaHCO3, with the NaCl concentration fixed at 0.5 mM ) and T* = T/T SS NaCl → NaHCO3 respectively, where initial currents of solutions are represented by I 0 with NaHCO 3 and NaCl subscripts respectively, and the time to reach steady-state current for NaCl → NaHCO 3 (arrow indicates EOF direction) is represented by T SS NaCl → NaHCO3 .
To study its effect on EOF with the 1 mM NaHCO 3 -NaCl solution pair, the electric field was increased by 50% to 187.5 V·cm −1 , as shown in Figure 9b. Good agreement is obtained between simulations and experimental results. The overall trend follows closely, and is similar to that of 1 mM NaHCO 3 -NaCl with an electric field of 125 V·cm −1 (see Figure 9a) and 0.5 mM NaHCO 3 -NaCl with an electric field of 125 V·cm −1 (see Figure 4). For 1 mM NaCl → 1 mM NaHCO 3 , despite increasing the electric field by 50% to 187.5 V·cm −1 , the NaHCO 3 concentration is similarly reduced to 0.750 mM (not affected based on KRF, see Equation (11)) and without complete displacement by 1 mM NaCl. For 1 mM NaHCO 3 → 1 mM NaCl with an electric field of 187.5 V·cm −1 , the NaCl equilibrium concentration is calculated to be 1.16 ± 0.050 mM (based on experimental data with Equation (12)), which is approximately the same as the case with the electric field of 125 V·cm −1 . Changing the electric field strength has no influence on the concentration evolution for the NaHCO 3 -NaCl solution pair system, except producing a much faster electroosmotic displacement process.
The effect of varying NaHCO 3 concentration was investigated experimentally by employing 0.5 mM, 3 mM or 5 mM NaHCO 3 , with the NaCl concentration fixed at 0.5 mM and an electric field of 125 V·cm −1 , as shown in Figure 10. For NaCl → NaHCO 3 , the NaHCO 3 concentration is reduced to 0.375 mM for different NaHCO 3 concentrations with 0.5 mM NaCl (according to Kohlrausch regulating function (KRF), see Equation (11)). With the displacement of 0.375 mM NaHCO 3 stayed unchanged at normalized interface X* = 0.75, the current-time curves stabilized at approximately the same current value despite different NaHCO 3 concentrations, as shown in Figure 10a. For NaHCO 3 → NaCl, NaHCO 3 still cannot displace into the microchannel. The diffusive-interface-induced evolution near the inlet vicinity results in NaCl equilibrium concentrations for 0.5 mM, 3 mM and 5 mM NaHCO 3 to be 0.640 ± 0.002 mM, 2.59 ± 0.02 mM and 4.02 ± 0.04 mM, respectively (calculated based on experimental data with Equation (12)). The increase in the NaCl equilibrium concentration was captured by the increase in the steady-state current of the current-time curve, when NaHCO 3 concentration was increased, see Figure 10b.
NaHCO3 concentrations for 0.5 mM and 1 mM NaCl reduce to 0.375 mM and 0.750 mM, respectively (according to KRF, see Equation (11)). As such, the steady-state current of the current-time curve was approximately doubled for 1 mM NaCl, as compared to 0.5 mM NaCl, see Figure 11a. However, complete displacement was realized when 0.1 mM NaCl was employed to displace 3 mM NaHCO3, with the current stabilized at the steady-state current of 0.1 mM NaCl (see Figure 11a). For NaHCO3 → NaCl, the NaCl equilibrium concentrations for 0.5 mM and 1 mM NaCl are 2.59 ± 0.02 mM and 2.70 ± 0.01 mM, respectively (calculated based on experimental data with Equation (12)). Since the NaHCO3 concentration was the same (fixed at 3 mM), the current-time curves stabilized at approximately the same current value for 0.5 mM and 1 mM NaCl, as shown in Figure 11b. While NaHCO3 completely displaced NaCl (with a concentration of 0.1 mM), the current stabilized at the steady-state current of 3 mM NaHCO3, see Figure 11b. Through lowering the NaCl concentration, the absolute zeta potential is increased, see Table 3; this enables stronger EOF for the displacement process. Investigation on the influence of varying NaCl concentration was experimentally conducted with 0.1 mM, 0.5 mM and 1 mM NaCl, where the NaHCO 3 concentration was fixed at 3 mM and electric field at 125 V·cm −1 , see Figure 11. For NaCl → NaHCO 3 , the NaHCO 3 concentrations for 0.5 mM and 1 mM NaCl reduce to 0.375 mM and 0.750 mM, respectively (according to KRF, see Equation (11)). As such, the steady-state current of the current-time curve was approximately doubled for 1 mM NaCl, as compared to 0.5 mM NaCl, see Figure 11a. However, complete displacement was realized when 0.1 mM NaCl was employed to displace 3 mM NaHCO 3 , with the current stabilized at the steady-state current of 0.1 mM NaCl (see Figure 11a). For NaHCO 3 → NaCl, the NaCl equilibrium concentrations for 0.5 mM and 1 mM NaCl are 2.59 ± 0.02 mM and 2.70 ± 0.01 mM, respectively (calculated based on experimental data with Equation (12)). Since the NaHCO 3 concentration was the same (fixed at 3 mM), the current-time curves stabilized at approximately the same current value for 0.5 mM and 1 mM NaCl, as shown in Figure 11b. While NaHCO 3 completely displaced NaCl (with a concentration of 0.1 mM), the current stabilized at the steady-state current of 3 mM NaHCO 3 , see Figure 11b. Through lowering the NaCl concentration, the absolute zeta potential is increased, see Table 3; this enables stronger EOF for the displacement process. Figure 10. Experimental observations of varying NaHCO3 concentration, i.e., 0.5 mM, 3 mM and 5 mM, with the NaCl concentration fixed at 0.5 mM and an electric field of 125 V·cm −1 . (a) NaCl → NaHCO3 (arrow indicates EOF direction), and (b) NaHCO3 → NaCl. Figure 11. Experimental observations of varying NaCl concentration, i.e., 0.1 mM, 0.5 mM and 1 mM, with the NaHCO3 concentration fixed at 3 mM and an electric field of 125V·cm −1 . (a) NaCl → NaHCO3 (arrow indicates EOF direction), and (b) NaHCO3 → NaCl.

Conclusions
Electroosmotic displacement flow involving multiple fluids exhibits EOF hysteresis, i.e., flow direction-dependent behavior. Thus far, no study has been conducted on the hysteresis effect for a solution pair system comprising of dissimilar ionic species and substantial pH difference; their rather different characteristics tend to result in a complex hysteretic behavior. In this investigation, the NaHCO3-NaCl solution pair was chosen as a model system to examine the hysteresis phenomenon.
The EOF of the NaHCO3-NaCl solution pair was carried out through current monitoring experiments. Finite element numerical simulations based on slip velocity model with modified wall boundary condition were performed for the evaluation of the displacement flow behaviors and ion concentration distributions. Quantitative agreements were achieved between experimental and simulation results.
For NaCl → NaHCO3 (the arrow indicates EOF direction), a concentration evolution of NaHCO3 happens. The displacement of the original residing and evolved NaHCO3 concentrations with high absolute zeta potential values (due to high pH values) by NaCl with low absolute zeta potential value causes the EOF to be reduced. As a result, NaCl is not fully displaced within the microchannel due to the gradient of the electromigrative flux cancelling that of the convective flux.
For NaHCO3 → NaHCO3, NaHCO3 cannot displace into the microchannel. This rather surprising and counter-intuitive outcome is a result of the stronger gradient of electromigrative flux than convective flux, as NaCl has low absolute zeta potential value (due to a low pH value) that generates slow EOF. Hence, evolution of the diffusive-interfaceinduced concentration of NaCl occurs.

Conclusions
Electroosmotic displacement flow involving multiple fluids exhibits EOF hysteresis, i.e., flow direction-dependent behavior. Thus far, no study has been conducted on the hysteresis effect for a solution pair system comprising of dissimilar ionic species and substantial pH difference; their rather different characteristics tend to result in a complex hysteretic behavior. In this investigation, the NaHCO 3 -NaCl solution pair was chosen as a model system to examine the hysteresis phenomenon.
The EOF of the NaHCO 3 -NaCl solution pair was carried out through current monitoring experiments. Finite element numerical simulations based on slip velocity model with modified wall boundary condition were performed for the evaluation of the displacement flow behaviors and ion concentration distributions. Quantitative agreements were achieved between experimental and simulation results.
For NaCl → NaHCO 3 (the arrow indicates EOF direction), a concentration evolution of NaHCO 3 happens. The displacement of the original residing and evolved NaHCO 3 concentrations with high absolute zeta potential values (due to high pH values) by NaCl with low absolute zeta potential value causes the EOF to be reduced. As a result, NaCl is not fully displaced within the microchannel due to the gradient of the electromigrative flux cancelling that of the convective flux.
For NaHCO 3 → NaHCO 3 , NaHCO 3 cannot displace into the microchannel. This rather surprising and counter-intuitive outcome is a result of the stronger gradient of electromigrative flux than convective flux, as NaCl has low absolute zeta potential value (due to a low pH value) that generates slow EOF. Hence, evolution of the diffusive-interfaceinduced concentration of NaCl occurs.
The aforesaid flow characteristics are independent of the applied electric field. However, through lowering the NaCl concentration, the absolute zeta potential value is increased, which enables EOF to increase for achieving complete displacement. While varying the NaHCO 3 concentration has negligible impact on the displacement process.
The outcomes of this investigation could provide a proper understanding of the flow behavior of inhomogeneous solutions with dissimilar pH and ion species for micro-/nanofluidic applications, such as isotachophoresis (ITP) and field amplified sample stacking (FASS).
Author Contributions: Y.C.L. and A.E.L. conceived the project, analyzed the results, wrote and reviewed the paper; A.E.L. carried out the experimental and numerical investigations. All authors have read and agreed to the published version of the manuscript.

Funding:
The research work was supported by Nanyang Technological University (NTU) with Grant No. 001274-00001.

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