Numerical Simulation of the Influence of Non-Uniform ζ Potential on Interfacial Flow

Zeta potential (ζ potential) is a significant parameter to characterize the electric property of the electric double layer (EDL), which is important at the solid–liquid interface. Non-uniform ζ potential could be developed on a chemically uniform solid–liquid interface due to external flow. However, its influence on the flow has never been concerned. In this investigation, we numerically studied the influence of non-uniform 2D ζ potential on the flow at the solid–liquid interface. It is found, that even without any external electric field and only considering the influence of 2D ζ potential distribution, swirling flow can be generated near EDL, according to the rotational electric volume force. The streamwise vortices, which are important in the turbulent boundary layer, are theoretically predicted in this laminar flow model when considering the 2D distribution of ζ potential, implying the necessity of considering the origin of streamwise vortices of the turbulent boundary layer from the perspective of electrokinetic flow. In addition, the ζ potential distribution can promote the wall shear stress. Therefore, more attention must be paid to shear-sensitivity circumstances, like biomedical, medical devices, and in vivo. We hope that the current investigation can help us to better understand the effect of charge distribution on interfacial flow and provide theoretical guidance for the development of related applications in the future.


Introduction
Electric double layer (EDL) is a kind of physical structure that exists generally in the two-phase interface or even three-phase interface (including solid, liquid, and gas).Taking the solid-liquid interface as an example, when the electrolyte solution flows through the channel, charges will be redistributed on the channel wall.Under the influence of the surface charge, ions in the solution with the opposite charge of the surface charge, i.e., counter-ions, will gather on the surface of the solid, causing the corresponding positive or negative charges to be reordered.Two thin layers with different electric properties can be formed on the interface, including a Stern layer and diffuse layer, constituting the important and common interface electric structure-EDL.
The Gouy-Chapman-Stern-Grahame (GCSG) model is a widely accepted EDL model, which holds that EDL is composed of a Stern layer and a diffuse layer [20][21][22][23].The ions in the Stern layer are considered to be adsorbed and fixed on the solid surface, while the ions in the diffuse layer but near the Stern layer cannot move freely due to the strong electrostatic force locally.As the distance from the Stern layer is beyond a certain level, the electrostatic force is not enough to restrain the thermal motion of the ions, and the ions diffuse randomly in the solution, which also shows that the diffuse layer as a dynamic structure can slide on the charged surface.The position where the ions are ready to move is called the shear plane.The ζ potential is the potential at the shear plane [23].
The ζ potential is a key parameter to characterize the electric structure in EDL.It plays an important role in interface science (e.g., interfacial reaction [24,25], self-assembly [26,27], etc.) and electrokinetic flow dynamics [28][29][30].However, limited by electrochemical analysis techniques [31][32][33][34] on ζ potential, the hypothesis of chemically uniform interfaces with unique ζ potential was usually considered in past research.Our understanding of local ζ potential and its influence on diverse disciplines is nearly blank.
In recent years, with the in-depth study of nanoscale interface phenomena and the development of new materials with complex surfaces and superstructures, researchers have gradually begun to pay attention to the influence of non-uniform surface charge on the electrical properties of the interface and explore the factors affecting the surface charge distribution and the corresponding potential.
On one hand, relevant research shows that chemically homogeneous slab channel walls exhibit non-uniform ζ potential under the action of an outflow field.For example, Lis et al. [35] using vibrational sum frequency generation (v-SFG) spectroscopy preliminarily reveal that surface potentials vary with external flows temporally.Subsequently, Werkhoven et al. [36] showed theoretically that a pressure-driven flow can induce a strong heterogeneous surface charge and ζ potential on a chemically homogeneous channel wall.Ober et al. [37] show that liquid flow along the surface of CaF 2 creates a reversible space charge gradient by v-SFG spectroscopy as well.In a current investigation, relying on a novel fluorescence photobleaching electrochemistry analyzer (FLEA), Meng et al. [38] reported a two-dimensional ζ potential induced by external flow at the solid-liquid interface.
On the other hand, non-uniform ζ potential can lead to unpredicted outcomes in diverse fields.Paratore et al. [39] demonstrated that the flow that surrounds a surface with nonuniform surface charge density (ζ potential as well) can induce complex flow patterns without a physical constraint.Yang et al. [40] revealed that gold superlattices with nanoscale variations in ζ potential distribution can significantly augment electrochemical reactions, underscoring the profound impact of ζ potential heterogeneity on electrochemical processes.
In the investigation, according to the previous experimental and theoretical reports, we hope to show how the non-uniform ζ potential distribution influences the electric volume force and its influence on the flow, through the numerical simulation method.

Theory
The momentum transport process of EK flow can be described by the Navier-Stokes equation with the electric volume force term, which is expressed as: where ρ f is the fluid density, P is the pressure on the fluid, µ is the dynamic viscosity.
is the velocity vector, with u 1 , u 2 and u 3 being the velocity components in x 1 , x 2 , and x 3 directions, respectively, and x1 , x2 and x3 the corresponding unit vectors in Cartesian coordinates (see Figure 1).
where ρ e is electric charge density, → E is the vector of electric field strength.In homogeneous dielectric media with incompressible approximation, the variation of ε is ignored, thus with ρ e and → E being: where φ is the electric potential in the EDL, ε = ε 0 ε r is electric permittivity with ε 0 being the vacuum permittivity, and ε r is the relative permittivity of fluid.When the potential changes slowly, the potential distribution in the EDL can be approximately described by Poisson-Boltzmann theory as: where k B is the Boltzmann constant, T is temperature, e is the elementary charge, and z i is the relative valence state of ions in electrolytes.λ = εk B T/2N A e 2 ∑ i c i z 2 i is the Debye length, which is a characteristic length of the EDL.N A is Avogadro constant, and c i is the concentration of ions in electrolytes.
where  is electric charge density,  ⃗ is the vector of electric field strength.In homogeneous dielectric media with incompressible approximation, the variation of  is ignored, thus  ⃗ =   ⃗ , with  and  ⃗ being: where  is the electric potential in the EDL,  =   is electric permittivity with  being the vacuum permittivity, and  is the relative permittivity of fluid.When the  = ( ,  ) potential changes slowly, the potential distribution in the EDL can be approximately described by Poisson-Boltzmann theory as: where  is the Boltzmann constant,  is temperature,  is the elementary charge, and  is the relative valence state of ions in electrolytes. =   2  ∑   ⁄ is the Debye length, which is a characteristic length of the EDL. is Avogadro constant, and  is the concentration of ions in electrolytes.Under steady state, divide both sides of Equation ( 1) by  , and then calculate the curl on both sides with incompressible fluid relation.Thus,, we have: where  ⃗ = ∇ ×  ⃗ is the vorticity,  =   ⁄ is the kinematic viscosity coefficient, and  ⃗ = ∇ ×  ⃗ =   +   +   is the driving term of vorticity equation.According to Maxwell's equations and assuming the magnetic field is steady, by substituting Equations ( 3) and ( 4) into  ⃗ , we find For convenience, let  =   ⁄ ,  =     ,  =      (with , ,  are either 1 or 2), from Equation (5), it is obtained Under steady state, divide both sides of Equation ( 1) by ρ f , and then calculate the curl on both sides with incompressible fluid relation.Thus, we have: where u is the vorticity, ν = µ/ρ f is the kinematic viscosity coefficient, and is the driving term of vorticity equation.According to Maxwell's equations and assuming the magnetic field is steady, by substituting Equations (3) and (4) into → T , we find For convenience, let where and ψ = 4k B T/ez i is thermal potential.Furthermore, after substituting Equations ( 8)- (10) into → F e and → T , we have and It can be seen when ζ is non-uniform in 2D, ζ i , ζ ij , and ζ ijk can be non-zero.All the three components in Equation (12) can be nonzero, and thus, the driving term of vorticity is: Therefore, based on Equation (6), the electric volume force → F e with curl can drive the fluid to produce a swirling flow.It is worth noting that, since there is no external electric field applied, we only consider the Coulomb force in Equation (1).The influence of other electric volume forces, e.g., according to electrothermal or dielectric effect, are believed to be absent in this investigation.

Numerical Simulation
Based on the theory above, we can make a reasonable prediction that in the steadystate flow field, such as pressure-driven flow, there may be a non-uniform two-dimensional ζ potential distribution at the solid-liquid interface, which will induce the generation of local electrostatic charge and EVF, resulting in a swirling flow near the interface.To validate this conjecture, we established a mathematical model of non-uniform ζ potential distribution and proceeded with a numerical simulation by Comsol Multiphysics on Equations ( 1)- (5).The influence of non-uniform ζ potential distribution at the solid-liquid interface on the flow is investigated under steady-state approximation with the value of relative tolerance is 1 × 10 −6 .

Pressure-Driven Laminar Flow
In this simulation model, the material is liquid water, with electric conductivity and pH values are σ = 100 µS/cm and pH = 9, respectively.The thickness of the EDL on the bottom of the microchannel, i.e., Debye length, is estimated to be 11 nm.The specific physical parameters are shown in Table 1.The basic flow is a pressure-driven steady flow with incompressible conditions.The bulk flow Reynolds number R e = Ud/ν < 1, indicating the basic flow is a typical laminar flow.Here, d = A/2(w + h) is the hydraulic diameter of the flow, U = Q/A is the bulk flow velocity, Q is flow rate, and A is the cross-sectional area of the computational regime.
In the laminar flow interface, the two sides of the rectangular model that are perpendicular to the x 1 direction are set to the inlet and outlet, respectively.At the entrance, the fully developed flow condition has been applied, with the average velocity U av = 6.78 × 10 −5 m/s.At the exit, the boundary condition is suppressing backflow (which adjusts the outlet pressure to reduce the amount of fluid entering the domain through the exit), with static pressure P 0 = 0 Pa.According to the symmetry of the geometric model and the rationality of the physical model, a no slip boundary condition has been imposed on the other four walls (including the upper wall, bottom wall, and two side walls), where the fluid velocity is zero.Electric volume force → F e has been applied to the entire domain of the model (see Equation (5), similar to potential φ, and → F e is dominant in the EDL at the bottom of microchannel), as shown in Figure 1.

Electrostatics Module
In this investigation, there is no external electric field applied, and the electric volume force → F e is entirely generated by the 2D distribution of ζ potential.To study the effect of the electrical volume force → F e due to the non-uniform ζ potential on the flow field, the electrostatic interface has been added into the simulation model.Furthermore, since the electric field is steady, the effect of the magnetic field is negligible.
The charge is conserved throughout the domain, and the electroosmotic potential φ calculated by Equation ( 5) was applied in the entire physical model.To further reduce the computational cost, we calculate the electric volume force by defining custom variables and analytic functions according to Equations ( 2)-( 5).
As mentioned above, in order to investigate the effect of non-uniform ζ potential distribution on interfacial flow, a mathematical model of ζ potential distribution (as plotted in Figure 2) was established to characterize the ζ potential distribution in the solid-liquid interface at the bottom of the microchannel only (as shown in Figure 1) based on the recent report of Meng et al. [38], as where of the model (see Equation ( 5), similar to potential , and  ⃗ is dominant in the EDL at the bottom of microchannel), as shown in Figure 1.

Electrostatics Module
In this investigation, there is no external electric field applied, and the electric volume force  ⃗ is entirely generated by the 2D distribution of ζ potential.To study the effect of the electrical volume force  ⃗ due to the non-uniform ζ potential on the flow field, the electrostatic interface has been added into the simulation model.Furthermore, since the electric field is steady, the effect of the magnetic field is negligible.
The charge is conserved throughout the domain, and the electroosmotic potential  calculated by Equation ( 5) was applied in the entire physical model.To further reduce the computational cost, we calculate the electric volume force by defining custom variables and analytic functions according to Equations ( 2)-( 5).
As mentioned above, in order to investigate the effect of non-uniform ζ potential distribution on interfacial flow, a mathematical model of ζ potential distribution (as plotted in Figure 2) was established to characterize the ζ potential distribution in the solid-liquid interface at the bottom of the microchannel only (as shown in Figure 1) based on the recent report of Meng et al. [38], as where  is the classical ζ potential value at the interface of water-fused silica [44].T could be calculated accordingly by Equations ( 3), ( 4), ( 11) and (12).
Finally, it is worth noting that due to the choice of steady-state study method, initial values of velocity, pressure, and potential are not considered when the boundary conditions are set.

Modeling and Meshing
The computational regime is a cuboid with dimensions in the order of micrometers.In the x 1 -x 2 plane, the geometric model is built around the origin O(0, 0).In the investigation, the potential φ (see Equation ( 5)) decreases exponentially with x 3 and is dominant only in the EDL at the solid-liquid interface.To capture the influence on flow dynamics, grid cells with nanoscale resolution are required.However, to be able to capture the evolution of flow, a large computational domain is required.To compromise these two requirements, the computation regime is L = 5 µm long, w = 4 µm wide, and h = 1 µm high.Thus, we can capture the influence of micrometer variation of ζ potential on the flow.In the modeling process, the geometric model is divided into three regions in the x 3 direction, including a large bulk region (h b = 0.9 µm) in the center of the channel and two symmetrical thin-layer regions (h t = 50 nm) near the bottom and upper wall.
In the meshing process, mapping operation is used in the x 1 -x 2 plane, and sweeping operation is carried out upward from the bottom to the upper wall in the x 3 direction and downward from the ceiling to the under wall in the opposite of the x 3 direction, respectively, in the two thin-layer regions.The grid can be denser in the target area by the arrangement of geometric sequence.Detailed grid distribution of each region is shown in Figure 3.
1 µm high.Thus, we can capture the influence of micrometer variation of ζ potential on the flow.In the modeling process, the geometric model is divided into three regions in the  direction, including a large bulk region (ℎ = 0.9 µm) in the center of the channel and two symmetrical thin-layer regions (ℎ = 50 nm) near the bottom and upper wall.
In the meshing process, mapping operation is used in the  - plane, and sweeping operation is carried out upward from the bottom to the upper wall in the  direction and downward from the ceiling to the under wall in the opposite of the  direction, respectively, in the two thin-layer regions.The grid can be denser in the target area by the arrangement of geometric sequence.Detailed grid distribution of each region is shown in Figure 3.During the sweeping operation in the  direction, the grid in the thin-layer region is divided down to the nanometer scale to ensure resolution.The minimum and maximum grid sizes are 1.33 nm and 2.66 nm, respectively.While in the bulk region, the grid is much coarser with a fixed size of 180 nm.In the  - plane, both bulk region and thin-layer During the sweeping operation in the x 3 direction, the grid in the thin-layer region is divided down to the nanometer scale to ensure resolution.The minimum and maximum grid sizes are 1.33 nm and 2.66 nm, respectively.While in the bulk region, the grid is much coarser with a fixed size of 180 nm.In the x 1 -x 2 plane, both bulk region and thinlayer region have the same minimum grid size and the same maximum grid size.In the x 1 direction, the minimum and maximum grid sizes are 6.68 nm and 20 nm, respectively.In the x 2 direction, the minimum and maximum grid sizes are 4.42 nm and 13.26 nm, respectively.Hexahedral grid elements are used for all regions in the model.Detailed parameters for complete grid statistics of the model are shown in Table 2.

Grid Independence Analysis
To ensure computational validity and stability, we conducted a grid independence analysis by systematically varying the grid size in all dimensions while keeping other parameters such as boundary conditions and solvers unchanged.The relative errors (δ) of the simulated velocity field between the models of different grid quantities (N) are plotted in Figure 4.It can be seen that δ decreases asymptotically towards zero, with the increase in grid quantity.The minimum value is 0.24%.The results show that when the grid quantity is equal to or beyond N = 7, 106, 000, δ is within 1%, indicating a negligible dependency of computational results on grid density [45].Consequently, to strike a balance between computational efficiency and precision, we adopted a grid resolution of N = 7, 106, 000 for subsequent analyses in this study.

Grid Independence Analysis
To ensure computational validity and stability, we conducted a analysis by systematically varying the grid size in all dimensions while rameters such as boundary conditions and solvers unchanged.The re the simulated velocity field between the models of different grid quanti in Figure 4.It can be seen that  decreases asymptotically towards zero in grid quantity.The minimum value is 0.24%.The results show that w tity is equal to or beyond  = 7,106,000,  is within 1%, indicating a ency of computational results on grid density [45].Consequently, to s tween computational efficiency and precision, we adopted a grid 7,106,000 for subsequent analyses in this study.

Pressure-Driven Basic Flow
In the previous investigations on interfacial flow, the influence o distribution and  potential distribution was not taken into account.T common Poiseuille flow) is driven only by pressure.

Pressure-Driven Basic Flow
In the previous investigations on interfacial flow, the influence of interfacial charge distribution and ζ potential distribution was not taken into account.The basic flow (e.g., common Poiseuille flow) is driven only by pressure.
In this case, the velocity field distribution of the basic flow is shown in Figure 5a, where x 1 * = x 1 /l, x 2 * = x 2 /l, x 3 * = x 3 /λ with characteristic length l = h/2 = 0.5 µm and Debye length λ = 11 nm.We applied a small volume rate (Q = 0.001 µL/h), where the basic flow velocity is 1.2 × 10 −4 m/s at 0.5 µm from the wall, to simulate a relatively large wall strain rate (240 s −1 ).Slice and quiver diagrams are used to show the distribution of the intensity (where the magnitude of intensity is represented by the color) and direction of the fluid velocity → u in the x 1 -x 3 planes at different streamwise positions.The direction of the fluid velocity after normalization at each position is represented by the vector arrow.Due to the absence of the external force, the velocity of the basic flow shows an approximately parabolic distribution in the spanwise direction near the wall, as can be found in Figure 5b with x 1 * = 4.In contrast, as shown in Figure 5c, since the EDL is extremely thin, the streamwise component of velocity, i.e., u 1 , shows a linear variation with the dimensionless vertical position x 3 * .This is an approximation of the parabolic velocity distribution of Poiseuille flow in the thin EDL.
of the intensity (where the magnitude of intensity is represented by the color) and direction of the fluid velocity  ⃗ in the  - planes at different streamwise positions.The direction of the fluid velocity after normalization at each position is represented by the vector arrow.Due to the absence of the external force, the velocity of the basic flow shows an approximately parabolic distribution in the spanwise direction near the wall, as can be found in Figure 5b with  * = 4.In contrast, as shown in Figure 5c, since the EDL is extremely thin, the streamwise component of velocity, i.e.,  , shows a linear variation with the dimensionless vertical

Uniform ζ Distribution
When considering the presence of uniformly distributed ζ potential (i.e., ζ = ζ 0 ) at the solid-liquid interface, based on Equation ( 5), it can be seen that there is an electroosmotic potential distribution in Figure 6a.It shows a uniform distribution in the x 1 -x 2 plane and exponential decay in the vertical direction.Accordingly, → F e are all in the negative x 3 direction with decreasing magnitudes from bottom to top, as shown in Figure 6b.tion with decreasing magnitudes from bottom to top, as shown in Figure 6b. ⃗ is exactly zero everywhere.In other words, the uniform distribution of the ζ potential induces irrotational  ⃗ .Under the influence of  ⃗ , the velocity distribution of the flow field is shown in Figure 6c.Compared with Figure 5a, it is obvious that the direction of the velocity  ⃗ has not changed, which still dominated by the  -directional component  .Furthermore, the parabolic velocity distribution in the  and  directions are still present.From the perspective of vortex dynamics, the source term → T in Equation ( 6) does not exist in either the basic flow or the flow field considering uniform ζ potential.Therefore, it does not disturb the basic flow, nor drive the fluid to produce vorticity, and cannot cause the vorticity to deform (stretch, tilt, twist, etc.).

One-Dimensional ζ Distribution
In this subsection, we consider if an 1D distribution of ζ potential with a certain slope in the x direction is formed due to nonuniform interface charge redistribution by the basic flow.Let A * = 0 in Equation ( 14), we have: This subsequently means that the x 2 component of → T, i.e., T 2 ̸ = 0, as shown in Figure 7c.
cal directions through the electrical volume force, except for  and  .As a result, the velocity distribution of the flow field in Figure 7e shows a clear difference from Figures 5a and 6c.The dominancy of  is no longer always true; in contrast,  becomes dominant when 0 <  * < 1.Compared with the basic flow, both  (Figure 5c) and | ⃗| have exhibited significant increment.The parabolic distribution in the  direction of basic flow (see Figure 5b) is no longer present when considering 1D ζ potential distribution.Although the flow with 1D ζ potential distribution can induce the electric volume force with F e1 ̸ = 0 and F e3 ̸ = 0, in this case, we still have the streamwise and vertical components of → T , i.e., T 1 = T 3 = 0, which could be seen more intuitively in Figure 7d.According to Equation (6), ω 1 = ω 3 = 0.There is no vorticity generated in streamwise and vertical directions through the electrical volume force, except for T 2 and ω 2 .As a result, the velocity distribution of the flow field in Figure 7e shows a clear difference from Figures 5a and 6c.The dominancy of u 1 is no longer always true; in contrast, u 3 becomes dominant when 0 < x 3 * < 1.Compared with the basic flow, both u 1 (Figure 5c) and → u have exhibited significant increment.The parabolic distribution in the x 2 direction of basic flow (see Figure 5b) is no longer present when considering 1D ζ potential distribution.

Two-Dimensional ζ Distribution
In this subsection, the influence of the 2D distribution of ζ potential on flow is investigated.In Equation ( 14), we take k = 222 1/m, A * = 0.1, ζ 0 = −36 mV, and k 2 = 1/π as an example.The ζ potential exhibits a 2D distribution on a solid-fluid surface, as plotted in Figure 2.
Based on Equations ( 2)-( 5), a 3D electroosmotic potential φ distribution is generated in the EDL at the bottom of the microchannel (as shown in Figure 8a), which in turn generates local electric charge and induces a non-uniform electric volume force.As shown in Figure 8b,  T shows an approximate symmetry with respect to the center of the channel, but in opposite directions (see Figure 8d).Two maxima can be observed at x 2 * = −1.5 and 1.5.

The fluid velocity distribution under the driven of rotational
→ F e is shown in Figure 9a.The vertical velocity component u 3 was observed near the bottom within one Debye length (i.e., 0 < x 3 * < 1).While u 1 component begins to dominate at x 3 * ≥ 1.This phenomenon suggests that even though most of the electric volume force is offset by pressure gradient, such a large volume force will still have a significant impact on the velocity field.At the height of x 3 * ≈ 2.5 and higher, → u is concentrated in the center of the microchannel in the spanwise to reach the maximum, which shows significant difference from those in Figures 5a, 6c and 7e.
In order to further explore the influence of 2D ζ potential on the interfacial flow, the velocity distribution in the x 1 -x 2 plane that one Debye length from the bottom (x 3 * = 1) is shown in Figure 9b by a surface diagram.It can be seen that the velocity → u of the flow, considering the effect of → F e caused by 2D ζ potential distribution, is about 3 folds higher than the basic flow.It exhibits an obvious increasing trend from upstream to downstream in the streamwise direction.In this plane, u 3 , rather than u 1 , is dominant among the three components of the flow velocity, indicating that the electric volume force (after being mostly balanced by the pressure term of Equation ( 1)) can cause the development of flow at the bottom wall.This implies, when considering the EK flow driven by 2D ζ potential, the momentum transport near the interface can be significantly promoted.
symmetry with respect to the center of the channel, but in opposite directions (see Figure 8d).Two maxima can be observed at  * = −1.5 and 1.5.
The fluid velocity distribution under the driven of rotational  ⃗ is shown in Figure 9a.The vertical velocity component  was observed near the bottom within one Debye length (i.e., 0 <  * < 1).While  component begins to dominate at  * ≥ 1.This phenomenon suggests that even though most of the electric volume force is offset by pressure gradient, such a large volume force will still have a significant impact on the velocity field.At the height of  * ≈ 2.5 and higher, | ⃗| is concentrated in the center of the microchannel in the spanwise to reach the maximum, which shows significant difference from those in Figures 5a, 6c and 7e.Streamwise vortices commonly exist in flow boundary layers, e.g., the well-known coherent vortex structures in the turbulent boundary layer.However, the origin of such vortices remains unclear.We hope to give an insight into this problem through the vision of electrokinetic mechanism, and provide a new approach to control the generation and evolution of vortices in boundary layers.

Discussion
In this section, the influence of the parameters (such as slope k, dimensionless amplitude A * , initial zeta potential ζ 0 , and wavenumber k 2 ) on the flow are systematically investigated by numerical simulations.

Influence of k
In this subsection, the influence of slope k on the flow is investigated at A * = 0.1, ζ 0 = −36 mV, and k 2 = 1/π.In Figure 10a, it is obvious that u 1 shows a linear increment with x 3 * , which is an approximation of the parabolic distribution caused by fluid viscosity in a thin layer.When the slope k gradually increases, a stronger variation of ζ potential along the streamwise direction is induced.The streamwise component of electric volume force gradually increases as well (see Figure 10b), accompanied by the increasing u 1 and velocity gradient ∂u 1 /∂x 3 .Thus, the wall shear stress is: which considering 2D ζ potential distribution can be significantly stronger than that of basic flow, as can be seen more clearly in Figure 10c.
= −36 mV, and  = 1  ⁄ .In Figure 10a, it is obvious that  shows with  * , which is an approximation of the parabolic distribution caused in a thin layer.When the slope  gradually increases, a stronger varia along the streamwise direction is induced.The streamwise component force gradually increases as well (see Figure 10b), accompanied by the i velocity gradient   ⁄ .Thus, the wall shear stress is: which considering 2D ζ potential distribution can be significantly stro basic flow, as can be seen more clearly in Figure 10c.All in all, it can be seen, when taking the influence of 2D ζ po wall shear stress can be significantly enhanced.For example, at  = −36 mV, and  = 1  ⁄ , the maximum  can be up to 1.07  9a.The increasing k 2 also decreases the streamwise velocity u 1 and its gradient ∂u 1 /∂x 3 (Figure 12c).Therefore, a smaller τ of flow can be concluded.
All in all, it can be seen, when taking the influence of 2D ζ potential into account, the wall shear stress can be significantly enhanced.For example, at k = 222 1/m, A * = 0.1, ζ 0 = −36 mV, and k 2 = 1/π, the maximum τ can be up to 1.074 Pa, which is about 2.5 folds larger than the basic flow.A stronger wall shear stress indicates a stronger interaction between wall and flow.This could be important in many shear-sensitive circumstances, especially on biological materials.One example is that the flow shear stress on cells and tissues may be larger than previously expected.Studying the different control parameters of 2D potential distribution can help us better understand the effect of charge distribution on interfacial flow, and according to Equation ( 16), we can even realize of adjustment of wall shear stress by changing these parameters.
tion between wall and flow.This could be important in many shear-sensitive circumstances, especially on biological materials.One example is that the flow shear stress on cells and tissues may be larger than previously expected.Studying the different control parameters of 2D potential distribution can help us better understand the effect of charge distribution on interfacial flow, and according to Equation ( 16), we can even realize of adjustment of wall shear stress by changing these parameters.

Influence of Other Parameters
In addition to the parameters of the zeta potential distribution discussed above, there are also several other parameters that affect the EK flow at the solid-liquid interface, for example, the electric conductivity and viscosity of the fluid.On the one hand, the electric conductivity is closely related to Debye length, which affects the distribution of electric potential  and charge density  in the EDL.This in turn affects electric volume force  ⃗ and the flow as delineated in Equations ( 2)-( 5).On the other hand, the influence of viscosity is also important.The EK flow induced by non-uniform ζ potential is intrinsically a balance between EVF and viscosity [46,47], primarily due to the balance of electric

Influence of Other Parameters
In addition to the parameters of the zeta potential distribution discussed above, there are also several other parameters that affect the EK flow at the solid-liquid interface, for example, the electric conductivity and viscosity of the fluid.On the one hand, the electric conductivity is closely related to Debye length, which affects the distribution of electric potential φ and charge density ρ e in the EDL.This in turn affects electric volume force → F e and the flow as delineated in Equations ( 2)-( 5).On the other hand, the influence of viscosity is also important.The EK flow induced by non-uniform ζ potential is intrinsically a balance between EVF and viscosity [46,47], primarily due to the balance of electric volume force and the viscous effect.The changing of viscosity could affect both the velocity profile and wall shear stress τ near the wall.However, a comprehensive exploration on these parameters accompanied with zeta potential variation are beyond the scope of this manuscript and will be studied in future works.

Conclusions
In this investigation, the influence of the non-uniform 2D distribution of ζ potential on the flow at the solid-liquid interface has been investigated numerically.Since the electric volume force → F e can be rotational, it further drives the fluid to produce a swirling flow.The

Figure 1 .
Figure 1.The computational domain, coordinate system, and the boundary conditions of the physical model in the numerical simulation.

Figure 1 .
Figure 1.The computational domain, coordinate system, and the boundary conditions of the physical model in the numerical simulation.
0 is the classical ζ potential value at the interface of water-fused silica [44].The slope k represents how quickly the ζ potential increases in the streamwise direction (with ζ < 0).The dimensionless amplitude A * represents the variation of ζ potential along the spanwise direction.The wavenumber k 2 determines the possible periodicity of ζ potential in the spanwise direction.
The slope  represents how quickly the ζ potential increases in the streamwise direction (with  < 0).The dimensionless amplitude  * represents the variation of ζ potential along the spanwise direction.The wavenumber  determines the possible periodicity of ζ potential in the spanwise direction.

Figure 3 .
Figure 3.The meshing of solution domain: (a) the thin-layer region (the schematic is scaled up to 10 times in the  direction) and (b) the large bulk region.

Figure 3 .
Figure 3.The meshing of solution domain: (a) the thin-layer region (the schematic is scaled up to 10 times in the x 3 direction) and (b) the large bulk region.

Figure 5 .
Figure 5. Velocity distributions.(a) 3D distribution of fluid velocity → u of basic flow.The magni-

→T
is exactly zero everywhere.In other words, the uniform distribution of the ζ potential induces irrotational → F e .Under the influence of → F e , the velocity distribution of the flow field is shown in Figure6c.Compared with Figure5a, it is obvious that the direction of the velocity → u has not changed, which still dominated by the x 1 -directional component u 1 .Furthermore, the parabolic velocity distribution in the x 2 and x 3 directions are still present.

Figure 6 .Figure 6 .
Figure 6.Potential, electric volume force, and velocity fields generated by uniform ζ potential with  = −36 mV.3D distribution of (a) electroosmotic potential , (b) electric volume force  ⃗ , and (c) fluid velocity  ⃗.From the perspective of vortex dynamics, the source term  ⃗ in Equation (6) does not exist in either the basic flow or the flow field considering uniform ζ potential.Therefore,

)
Taking k = 222 1/m and ζ 0 = −36 mV as an example, the 1D distribution of ζ potential see Figure 7a, and the distributions of φ and → F e induced by this ζ potential are shown in Figures 7b and 7c, respectively.Based on the increase in |ζ| in the x 1 direction, the strength of → F e also shows an increasing trend from upstream to downstream.Compared with Figure 6b, although the x 3 component of → F e is still dominant, the x 1 component of the → F e is not zero anymore.

Figure 7 .
Figure 7. Electric field and flow field generated by 1D ζ potential with  = 222 1/m and  = −36 mV.(a) 1D distribution of ζ potential.3D distribution of (b) electroosmotic potential , (c) electric volume force  ⃗ , (d) the curl of electric volume force  ⃗ , and (e) fluid velocity  ⃗ under this 1D ζ potential distribution.4.4.Two-Dimensional ζ Distribution In this subsection, the influence of the 2D distribution of ζ potential on flow is investigated.In Equation (14), we take  = 222 1/m,  * = 0.1,  = −36 mV, and  = 1  ⁄ as an example.The ζ potential exhibits a 2D distribution on a solid-fluid surface, as plot-

F
e exhibits strong x 3 -directional component.The directions of → F e are mainly pointing to the bottom wall of the microchannel.The strength of → F e generally shows a decaying trend along x 3 , reaching a maximum at one Debye length above the bottom wall.The strength of → F e also exhibits variations along the streamwise direction and the spanwise direction.It increases gradually downstream and decreases from the center of the flow field, forming an approximate symmetry distribution along the spanwise direction.Different from the → F e induced by 0D (Figure 6b) and 1D (Figure 7c) distributions of ζ potential, in the model of 2D ζ potential, → F e can be rotational in 3D.All three components of → T can be non-zero, as shown in Figure 8c,d.→ T increases in the streamwise direction, which is consistent with that of → F e .In the spanwise, →

Figure 9 .
Figure 9. Velocity and vorticity fields generated by 2D ζ potential with  = 222 1/m,  * = 0.1,  = −36 mV and  = 1  ⁄ .(a) 3D distribution of fluid velocity  ⃗.(b) 2D velocity distribution of flow considering and without considering ζ potential distribution at  * = 1.(c) 3D distribution of the curl of fluid velocity  ⃗, denoted by | ⃗ * |, with  ⃗ * =  ⃗/  ⃗ .Here,  ⃗ is the maximum | ⃗| of basic flow.In order to further explore the influence of 2D ζ potential on the interfacial flow, the velocity distribution in the  - plane that one Debye length from the bottom ( * = 1) is shown in Figure9bby a surface diagram.It can be seen that the velocity | ⃗| of the flow, considering the effect of  ⃗ caused by 2D ζ potential distribution, is about 3 folds higher than the basic flow.It exhibits an obvious increasing trend from upstream to downstream in the streamwise direction.In this plane,  , rather than  , is dominant among the three components of the flow velocity, indicating that the electric volume force (after being mostly balanced by the pressure term of Equation (1)) can cause the development of flow at the bottom wall.This implies, when considering the EK flow driven by 2D ζ potential, the momentum transport near the interface can be significantly promoted.

Figure 11 .
Figure 11.Influence of ζ 0 on the flow, where x 1 * = 4 and x 2 * = 0, k = 222 1/m, A * = 0.1, and k 2 = 1/π.Plot of streamwise velocity component, u 1 , versus dimensional height x 3 * at different ζ 0 .5.4.Influence of k 2 At last, we discuss the influence of spanwise distributions of ζ potential on the flow, under k = 222 1/m, A * = 0.1, and ζ 0 = −36 mV.Larger k 2 indicates a wavier distribution of ζ potential in the x 2 direction, which induces a stronger gradient of ζ potential and more peak value positions of

Table 1 .
Properties of the material.

Table 2 .
Detailed information of the grid in the model.