Temperature-gradient induced massive augmentation of solute dispersion in viscoelastic micro-flows

Enhancing solute dispersion in electrically actuated flows has always been a challenging proposition, as attributed to the inherent uniformity of the flow field in absence of surface patterns. Over the years, researchers have focused their attention towards circumventing this limitation, by employing several fluidic and geometric modulations. However, the corresponding improvements in solute dispersion often turn out to be inconsequential. Here we unveil that by exploiting the interplay between an externally imposed temperature gradient, subsequent electrical charge redistribution and ionic motion, coupled with the rheological complexities of the fluid, one can achieve up to one order of magnitude enhancement of solute dispersion in a pressure-driven flow of an electrolyte solution. Our results demonstrate that the complex coupling between thermal, electrical, hydro-dynamic and rheological parameters over small scales, responsible for such exclusive phenomenon, can be utilitarian in designing novel thermally-actuated micro and bio-microfluidic devices with favorable solute separation and dispersion characteristics.


Introduction
Integrating multiple fluidic processes into a single platform has become progressively important in modern lab-on-a-chip devices where separation and mixing often turn out to be two most critical processes (Ghosal 2004;Hunter 1981;J H Masliyah & S Bhattacherjee 2006;R. F. Probstein 1994;Stone et al. 2004;Stroock et al. 2002;Whitesides 2006). With rapid advancement in micro-fabrication technologies, a large number of research efforts have been dedicated towards developing strategies for improved fluidic mixing or separation (Anderson et al. 2000;Chang & Yang 2008;Ghosh & Chakraborty 2012;Glasgow et al. 2004;Karniadakis G. 2005;Sugioka 2010;Zhang et al. 2006). Towards achieving enhanced mixing in micro-devices, diffusion and dispersion are undoubtedly the two most common phenomena. Accordingly, significant research interest in this domain has evolved over the past years, with a vision of employing different flow actuating mechanisms as well as geometric alterations in the fluidic pathways, so as to achieve the desired functionalities.
Although flow actuation using electric fields has wide spectrum of applications in both engineering and medical domain (Bandopadhyay & Chakraborty 2012;Becker & Gärtner 2000;Berli 2010;Das et al. , 2018Das S., and Chakraborty S. 2007;Garcia et al. 2005;Haeberle & Zengerle 2007;Mandal et al. 2012;Mark et al. 2010;Nguyen et al. 2013;Ohno et al. 2008;Van Der Heyden et al. 2005;Zhao 2011), one of the major aspects of the classical electroosmotic flow (EOF) in presence of homogenous interfacial conditions is the existence of the uniform velocity profile which arises when the electrical double layer (EDL) becomes very thin compared to the channel dimension (Ghosal 2004;J H Masliyah & S Bhattacherjee 2006). This results in a plug-type velocity distribution, thus reducing the extent of mixing significantly.
Hydrodynamic dispersion is the band broadening of the solute which mainly arises from the non-uniformity in the flow field (Ajdari et al. 2006;Arcos et al. 2018;Aris 1956Aris , 1959Barton 1983;Chatwin 1970Chatwin , 1975Chatwin & Sullivan 1982;Datta & Ghosal 2008;Dutta 2008;Ghosal 2006;Ghosal & Chen 2012;Jansons 2006;Mazumder & Das 1992;Ng & Yip 2001;Rana & Murthy 2016;Smith 1982;Sounart & Baygents 2007;Taylor 1953;Watson 1983;Zholkovskij & Masliyah 2004). Under ideal circumstances, velocity profile of EOF does not contribute to the shear-induced axial dispersion because of the flatness of the velocity profile as opposed to the case of Poiseuille flow (which is parabolic in nature) (Gaš et al. 1997; The distortion of the local equilibrium of ions in the EDL creates a departure from the classical Boltzmann distribution of ions which has direct consequences on the potential distribution. This alteration in charge distribution influences the flow dynamics via electrokinetic forcing, where the intricate coupling between thermal and electrical effects are already prevalent through the aforesaid property variation and an additional dielectrophoretic force. Further, in the absence of any external electric field, the net ionic current turns out to be zero. This condition is itself another source of non-linearity in the analysis where both conduction and streaming current undergo drastic alteration under the influence of finite temperature difference. In addition, it is noteworthy to mention that, when an electrolyte solution is subjected to an imposed temperature gradient, a thermo-electric field is induced by virtue of the movement of ions in response to the thermal driving force, commonly known as Soret effect (Dietzel & Hardt 2016Ghonge et al. 2013;Maheedhara et al. 2018;Zhang et al. 2019;Zhou et al. 2015). Moreover, additional form of thermo-electric field can be induced within the system because of the diffusivity difference between the ions even if their Soret coefficients remain the same. Besides, the mode of application of thermal gradient can cause significant rearrangement of ions within the EDL and hence, the subsequent potential distribution for a transverse temperature gradient may not necessarily be same with that of longitudinal temperature, altering the hydrodynamics in a rather profound manner. Considering the aforementioned intricacies in coupling thermal, electrical and hydro-dynamical effects in micro-confinements, research efforts towards addressing various aspects of thermo-solutal convection of electrolyte solutions have turned out to be relatively inadequate, despite having widespread applications in processes like water treatment, charge separation, zeta-potential determination, waste heat recovery, and energy conversion (Barragán & Kjelstrup 2017;Dietzel & Hardt 2016Jokinen et al. 2016;Li & Wang 2018;Sandbakk et al. 2013;Würger 2008Würger , 2010Xie et al. 2018). As such, the research focus in this domain has been directed primarily towards incorporating non-isothermal effects as a secondary force in the alteration of hydrodynamics of simple fluids (Chakraborty 2006;Garai & Chakraborty 2009;Huang & Yang 2006;Keramati et al. 2016;Maynes & Webb 2003;Sadeghi et al. 2011;Sánchez et al. 2018;Tang et al. 2003;Xuan et al. 2004a;Xuan 2008;Xuan et al. 2004b;Yavari et al. 2012). Therefore, except for some limited physical scenarios, however, such an exclusive effect has not been utilized to a significant practical benefit (Dietzel & Hardt 2016Zhang et al. 2019).
Recently, incorporation of thermal gradient has emerged as an alternative tool in augmenting dispersion where interplay between thermal and electrical effects over small length scales, almost exclusively, dictates the flow physics (Chen et al. 2005;Mukherjee et al. 2019;Sánchez et al. 2018). In addition, it may be noted that the emergence of new generation medical devices, complex bio-fluids have more prominently come into the paradigm of microfluidics (Berli 2010;Berli & Olivares 2008;Olivares et al. 2009;Zhao & Yang 2011. Such fluids exhibit strikingly distinct behavior compared to the fluids obeying Newton's law of viscosity (Brust et al. 2013;De Loubens et al. 2011;Fam et al. 2007;Moyers-Gonzalez et al. 2008;Owens 2006;Silva et al. 2017;Vissink et al. 1984). Some recent studies have demonstrated that the constitutive behavior of these biological fluids has close resemblance with the rheology of viscoelastic fluids and therefore, the inclusion of fluid rheology and viscoelasticity in dispersion characteristics has gained significant attention lately (Arcos et al. 2018;Brust et al. 2013;Hoshyargar et al. 2018;Mukherjee et al. 2019). While considering the thermally induced electrokinetic flow of viscoelastic fluids, an additional source of non-linearity crops in as mediated by the constitutive behavior of the fluid (Afonso et al. 2009(Afonso et al. , 2013Coelho et al. 2012;Ferrás et al. 2016;Ghosh et al. 2016;. Moreover, the degree of viscoelasticity, which is determined by using physical properties like fluid viscosity and relaxation time, is a strong function of the prevalent thermal gradient, augmenting the complexity of the problem to a large extent (Bautista et al. 2013;Mukherjee et al. 2019).
To the best of our knowledge, dispersion characteristics of thermally induced electrokinetic transport of complex fluids in microfluidic environment, where the temperature gradient is solely used for flow manipulation, has not been addressed in the literature. Here, we report the effect of an external temperature gradient on the dispersion characteristics of an electrolyte solution in a parallel plate microchannel. We subsequently discuss the charge redistribution upon applied thermal gradient, subsequent perturbation on the fluid motion and its implications on hydrodynamic dispersion, considering both Newtonian and viscoelastic fluids.
Our results reveal that by combining some of the electrokinetic, thermal and fluidic parameters coupled with rheological aspects, it is possible to achieve massively augmented solute dispersion, while for some combinations significant enhancement in streaming potential (compared to the solely pressure-driven flow) can also be obtained. We believe that the present analysis can be used as a fundamental basis in the design of thermally actuated micro and biofluidic devices demanding improved solute dispersion where the interplay between electromechanics, thermal effect, hydrodynamic and rheological aspects in narrow confinement can be coupled together to a beneficial effect.

Problem formulation
We consider non-isothermal electrokinetic flow of a binary 1:1 symmetric electrolyte solution through a parallel plate microchannel. We choose rectangular Cartesian co-ordinate system where x and y co-ordinates represent longitudinal and transverse directions respectively while the origin is placed at the centreline of the channel. Length scales in the two directions are l and h respectively where the half-channel height (h) is very small compared to the channel length (l), i.e. h << l or β = h / l << 1. We have employed two different types of thermal gradients, (a) Case 1: axially applied temperature gradient and (b) Case 2: temperature gradient applied in the transverse direction. First, we will briefly discuss about the formulation of the main focus of the present analysis, i.e. obtaining the hydrodynamic dispersion coefficient which directly depends on the flow velocity which in turn is influenced strongly by the charge distribution modulated by external temperature gradient. Hence, formulation regarding temperature field is presented first and the resulting potential distribution and velocity fields are presented subsequently. Final part of the present analysis is the inclusion of fluid rheology to examine its role in altering the dispersion characteristics.

Dispersion coefficient
We consider the hydrodynamic dispersion characteristics of a combined pressure-driven and thermal-gradient driven electrokinetic flow of an electrolyte solution. According to the definition of band-broadening phenomenon, hydrodynamic dispersion coefficient (D ̅ eff ) depends on the rate of change of variance of solute displacement band as its ability to incorporate any change in the variance of the mean concentration. For a band of non-adsorbing solute flowing in a rectangular channel, the velocity of the centre mass becomes equal to the mean velocity of flow, i.e. avg u dx dt ′ = . Using the descriptions of h' and x', one can rewrite the expression of dispersion coefficient as Here, u avg is the mean velocity averaged cross-sectionally. As reported in the literature, (van Deemter et al. 1956), the plate height is related to the mean velocity a ( ) In equation (2), h min is the minimum plate height for a given flow condition and D is diffusivity.
On the right side, the contribution of the molecular diffusion is represented by the first term while the second term indicates the contribution due to the non-uniformity in the flow field. By following some recent studies (Arcos et al. 2018;Hoshyargar et al. 2018 where ρ, C p and k are density, specific heat capacity and thermal conductivity of the fluid respectively. The left hand side of equation (5) (5) can be written as is the viscosity of the fluid. Here, we have neglected the variation of ρ and C p with temperature in obtaining the temperature distribution, whereas other thermo-physical properties like viscosity (µ), electrical permittivity (ε) and thermal conductivity (k) are considered to be temperature-dependent. The following forms of temperature dependences are assumed: In equation (6)  ~ O(10 -10 ) Fm -1 , D ~ 10 -9 m 2 s -1 , κ ~ O(10 7 ) m -1 becomes O(10 -9 ) which is further multiplied by the quantity β 2 where β << 1 and hence, the heat generation due to streaming field can be neglected. Also, most of the terms in viscous dissipation component involve β 2 and diminish for β << 1 . The remaining term is , which in our analysis comes out to be O(10 -4 ) and remains insignificant in determining the temperature field.

Case (a): Temperature gradient applied in the axial direction
For axially applied temperature gradient, the simplified energy equation subjected to the aforesaid assumptions is given below where the dimensionless form of thermal conductivity is written as and solution: ( The comparison between the exact solution and asymptotic solution for temperature distribution is presented in figure 2 where the lines show the results of exact solution with symbols representing predictions of asymptotic approach. With increasing C k , enhanced temperature sensitivity of the fluid thermal conductivity results in a departure from the linear variation of temperature in the axial direction. At higher C k (i.e. at C k = 10), this distribution becomes parabolic in nature where a deviation between the asymptotic and exact solution can be noticed (as shown in the inset). However, as reported in the literature, the relative change of thermal conductivity with temperature in typical electrolyte solution is (Dietzel & Hardt 2017) for which the value of C k turns out to be of the order of unity and hence, a linear temperature distribution can be assumed safely as an approximation. Throughout our analysis, we maintain C k = 1 while presenting the results.

Potential distribution
In contrast to the conventional electrokinetics problem, in presence of thermal gradient, ions no longer remain in equilibrium and one cannot consider Boltzmann distribution assumption while obtaining the potential distribution. Here, we need to find the ionic number concentration first In equation (12), Ti D and * i µ are the thermophoretic and electrophoretic mobilities respectively with ( ) Here, Ti S is Soret coefficient of ions defined as ( ) x y ψ being the potential induced within EDL. Considering this, the distribution of i n subjected to symmetry condition at the channel centreline (i.e., at ȳ = 0, 0 i n y ∂ ∂ = ) and number density being equal to bulk number concentration in electroneutral region ( i i n n ∞ = at 0 Since this shows an exponential dependence, it may seem like similar to the well-known Boltzmann distribution. However, the bulk ionic concentration ( ) i n ∞ is not constant, instead it is varying axially in presence of the axial temperature gradient. To obtain this dependence, one need to equate the first term of right side of equation (13) with zero, so that Since electro-neutrality is prevalent in the bulk (i.e. 0 x φ ∂ ∂ = ), equation (15) Contrary to the traditional electrokinetic studies, here EDL thickness (i.e. inverse of eff κ ) no longer remains constant and becomes a function of the axial co-ordinate ( ) . For small values of surface potential and γ (i.e. small imposed temperature difference), we can use Debye- ; along with the exponential term being linearized as In this context, it is worth mentioning that in presence of a thermal gradient, the surface potential (or zeta-potential) no longer remains constant. Instead, it becomes a strong function of the imposed temperature difference by following a linear relationship. Some previous experimental studies have demonstrated this temperature-dependence of zeta potential where it not only depends on temperature but also other factors like surface reactions (Ghonge et al. 2013;Ishido & Mizutani 1981;Reppert 2003;Revil et al. 2003). However, this surface reaction has relatively weaker influence on altering the non-isothermal flow dynamics, as highlighted in a recent study by Dietzel et al. (Dietzel & Hardt 2016). Consideration of this surface equilibrium and other associated reaction kinetics will unnecessarily complicate the present theoretical framework and hence, only linear dependence of zeta-potential with temperature is considered, i.e. at 1 y = ± , being the temperature sensitivity of zeta potential. Equation (18) subjected to the above boundary condition results in the following two equations for the potential The coefficients of equation (20) are shown in Section A1 of the supplementary material where the modified potential distribution for temperature-dependent of zeta potential is also presented.
Reported experimental results of temperature dependence of zeta potential (Reppert 2003) is depicted in figure 3a (i) which clearly shows that zeta potential (ζ) is indeed strongly dependent on temperature by following a linear relationship. These data are fitted in the form of (this is the same form of zeta potential variation with temperature which we present as and this approximates quite well with the experimental data. As can be seen from this figure, depending on electrolyte concentration, this sensitivity with temperature is increased almost twice as the slope (m) of fitted line increases from 0.0112 to 0.0214. In the dimensionless form, the value of C ζ turns out to be varying from ~ 3 to 6 (approximately).
Accordingly, in our analysis, we have chosen C ζ to lie between 0 and 4 while presenting the results, with C ζ = 0 representing the case of constant zeta potential. Figure 3a (ii) shows the potential distribution in the transverse direction where the inset clearly shows the effect of C ζ on the potential distribution. Effect of C ζ is noticeable only in close vicinity of channel walls (where EDL is located) while remains ineffective in the electro-neutral region. As C ζ is increased, more is the ionic redistribution within EDL and the magnitude of zeta potential may increase up to twice (at γ = 0.1 as shown in the inset) as observed in figure 3a (iii) for C ζ = 4.
Another key factor in altering the potential distribution is C ε representing the change in electrical permittivity with temperature. Increasing C ε from 1 to 10 indicates significant reduction in electrical permittivity whose effect is reflected in the charge distribution via the permittivity-induced component ( ) As observed in figure 3a (iv), increasing C ε FIGURE 3a. (i) Reported experimental results of zeta potential variation with temperature (Reppert 2003), (ii) potential distribution in the y-direction for both constant and temperaturedependent zeta potential, variation of the same for (iii) different C ζ and for (iv) different C ε .
creates a deviation in potential profile close to the channel walls while remains unaffected in the bulk.

Velocity distribution
Because of the application of the externally imposed temperature gradient (∆T), a thermo-electric field is induced within the microchannel by virtue of the difference in Soret coefficients between cations and anions. Besides, one can expect another form of thermoelectricity which can be induced due to the difference in ionic diffusivities of the ions even if their Soret coefficients remain same. Also, since the thermo-physical properties of the fluid are strongly dependent on temperature, there will be drastic alteration in the hydrodynamics in presence of ∆T. For example, the variation of electrical permittivity with temperature is manifested through the contribution of the dielectrophoretic body force with concomitant induction of axial pressure gradient thereby strongly influencing the advection motion of fluid. For steady, incompressible flow, the velocity distribution for a Newtonian fluid is described by continuity equation 0 ∇ ⋅ = v along the following momentum equation The left hand side of equation (21) is zero because of negligible inertial effect consideration which occurs only when the Reynolds number (Re) associated with flow is very less compared to unity, i.e. Re << 1. Here, p is the hydrodynamic pressure, In equation (22) (23) is now solved by following the aforesaid asymptotic approach where the governing equations at different orders of perturbations are as follows To obtain the flow field, equation (24) is subjected to no-slip condition at the surface ( ) = 0 at 1 u y = and symmetry condition ( ) 0 at 0 u y y ∂ ∂ = = at the channel centreline. Using this, the velocity distribution are given in the following two equations Equation (25) represents the perturbation-free flow field which is the well-known expression for combined pressure-driven and streaming field induced electrokinetic flow while equation (26) indicates the contribution of the thermal perturbation to the flow field. The coefficients of equation (26) and 0  1 48 The coefficients of equation (28) can be found in Section A4 of the supplementary material.
Here, it is necessary to mention about the two parameters involved in the final expression of the streaming potential described by equations (27)-(28), χ and C D . χ represents the molecular diffusivity difference between the co-ions and counter-ions represents the sensitivity of diffusivity of ions with temperature (here diffusivity D is assumed to be obeying the following relationship ( ) [ ] Similarly, the dispersion coefficient can also be expressed as where the estimation of 1 D (i.e. dispersion coefficient at O(γ)) involves the knowledge of O(γ) flow velocity 1 u .

Case (b): Temperature gradient applied in the transverse direction
For temperature difference applied in the transverse direction, considering all previouslymentioned assumptions and simplifications (for the case of axially applied thermal gradient), the governing equation of temperature field reads as and the corresponding temperature profile (asymptotic solution) is given by Now, the governing equation for the transport of ionic species reads as which results the following distribution of the ionic number concentration within the EDL This distribution is now used in the Poisson equation to obtain the potential distribution. The governing equations for the potential distribution are as follows For potential distribution, unlike the case of axial temperature gradient, no closed form solution is possible. So, one needs to either employ the asymptotic approach (approximate solution) or can solve it numerically.
The coefficients of equation (36) are reported in Section B1 of the supplementary material.
Comparing the potential distribution profiles for longitudinal and transverse thermal gradients one can understand that the basic difference is the introduction of an asymmetry through the temperature distribution which has its immediate effect on the potential distribution within the EDL. With increasing C ζ , the magnitude of surface potential gets amplified significantly in the hot region while remaining unaffected in the cold region thus creating more asymmetry in the y-direction. This further influences the fluid advective motion through electrokinetic forcing thereby altering strongly the velocity distribution. Knowing the potential distribution, the flow field can now be obtained from the following governing equations ( ) The solution of these two equations subjected to no-slip boundary condition at the surfaces yields FIGURE 3b. Potential distribution in the y-direction (at γ = 0.1) for transverse temperature gradient. Inset shows the zoomed view towards the top wall (i.e. close to ȳ = 1).
The coefficients of equation (40) and ( ) with the coefficients shown in Section B3 of the supplementary material.

Fluid rheology
Now, we focus on the alteration of hydrodynamics and associated streaming field caused by the One important thing to note here that, in presence of temperature gradient, not only fluid viscosity but also fluid relaxation time start to become temperature-dependent and interestingly, the extent of viscoelasticity of a fluid is governed by these two parameters -viscosity and relaxation time. In equation (43)

Results and discussions
Since this analysis involves large set of parameters, numerous results can be obtained by combining all pertinent parameters. However, for improved readability, we have highlighted some key results involving velocity distribution, induced streaming potential, volumetric flow rate, and finally the dispersion coefficient (which in turn depends on previous parameters). For representing the results, we have employed the following ranges of the involved parameters: 0 ≤ γ (thermal perturbation parameter) ≤ 0.1, 10 -1 ≤ C µ (sensitivity coefficient of viscosity) ≤ 10, 10 -1 ≤ C ε (sensitivity coefficient of electrical permittivity) ≤ 10, C k (sensitivity coefficient of thermal First we have shown the validity of range of some parameters involved in our analysis. As already mentioned, both approximate and exact solutions are obtained for axially applied thermal gradient. Figure 4(i) depicts the comparison of volumetric flow rates between the asymptotic (approximate) and exact solutions as a function of C µ which denotes the sensitivity of fluid viscosity with temperature. Here, C µ is varied over two decades ranging from 10 -1 to 10 1 while keeping other parameters constant. As evident, asymptotic solution closely approximates the exact solution up to C µ = 3. Beyond this critical value of C µ , large deviation between these two solutions is observed with asymptotic approach under-predicting the results significantly at higher C µ . In figure 4(ii), comparison of the same has been shown for varying γ, which denotes the ratio of the imposed temperature difference with respect to the reference temperature, i.e. a quantification of degree of thermal perturbation imposed. In this context, it is necessary to note that our asymptotic analysis is corrected up to first order in γ. Therefore, our findings are able to capture the linear thermal effect only. For lower C µ , the asymptotic solution approximates the exact solution reasonably well while at higher C µ , deviation takes place between the two solutions with asymptotic solution under-predicting volumetric flow rate beyond γ = 0.042. In the similar way, comparisons between asymptotic and numerical solutions have been made for the case of transverse temperature gradient which for the sake of conciseness are presented in Section D of the supplementary material.

Effect of axial temperature gradient
In figure 5(a), the variation of streaming potential ratio (E r ) is plotted against some key parameters. Here, streaming potential ratio (E r ) is defined as the ratio of the streaming potential induced due to the combined action of externally imposed thermal and pressure gradient to that counter-ions) is higher compared to the anions which leads to more migration of counter-ions towards upstream. As a result, there will be an accumulation of counter-ions in the upstream while in the downstream there will be more co-ions. This segregation of ions in different axial locations creates an axial separation between them thus creating stronger induced streaming FIGURE 5a. Dependence of streaming potential ratio (E r ) with (i) ∆S̅ T , (ii) χ, (iii) C µ , (iv) C ε , (v) γ and (vi) 0 κ respectively. field, while negative value of χ means more mobility of anions thereby leads to more accumulation of co-ions in the upstream section resulting weaker streaming field. Here we have shown the variation of streaming potential ratio (E r ) for χ ranging from -0.2 to 0 (figure 5a (ii)).
Here, the value of χ is ion-specific. For typical electrolyte solutions like aqueous KCl, NaCl, LiCl, the value of χ ranges between ~ -0.3 to 0 (approximately) (Zhang et al. 2019) and we have chosen the range of χ accordingly while presenting the results.
On closely observing the governing equations involving velcoity profile and induced streaming field one can understand that the contribution of C µ comes through the viscous resistance term in the fluid advective motion and the subsequent advective current calculation involved in the electroneutrality condition. Since C µ indicates the temperature sensitivity of viscosity with temperature, increasing C µ means viscosity becomes more susceptible to any change in temperature thus resulting strong reduction in fluid viscosity. Therefore, the viscous resistance to the flow reduces to a great extent. Thus, the induced streaming field due to migration of ions upon applied pressure gradient and the streaming field due to migration of ions upon induced concentration gradient assist each other (both are induced in the direction from the hot region to the cold region) resulting significant augmentation in the streaming potential ratio (E r ). Here, increasing C µ from 1 to 10 results in ~ 3 times (shown by the green colored solid line) augmentation of the streaming potential ratio (E r ), as visible from figure 5a (iii). Now, examining the equation describing the potential distribution (equation (17)), it can be observed that unlike the conventional electrokinetic problem, here the EDL thickness (λ D ) no longer remains constant. Instead, the effective EDL thickness (λ Deff ) becomes a strong function of temperature and departs significantly from its reference value (of isothermal condition) in the following way: Keeping other parameters fixed, with increasing C ε (which is the temperature sensitivity of electrical permittivity) streaming potential ratio (E r ) decreases sharply by following an exponential thinning behavior. Since electrical permittivity decreases with temperature, the thickness of EDL also decreases which results lesser penetration of diffuse layer of EDL to the bulk, so more is the region of electroneutrality (which means the region where there is equal number of counter-ions and co-ions). Now, the strength of streaming current involved in the streaming potential estimation depends on this degree of penetration. In the region 10 -1 ≤ C ε ≤ 1, the rate of reduction of E r with C ε is not significant, it decreases slightly from 1.2 to 1.1. However, this reduction gets amplified beyond C ε ~ 1, E r decreases sharply and falls to ~ 0.2 at C ε = 10 (figure 5a (iv)).
For fixed value of the aforesaid factors like ∆S̅ T , χ, C µ , C ε , since this analysis is restricted up to order γ, it is expected to get linear dependence of perturbation parameter on the flow field and streaming potential. For constant zeta potential (C ζ = 0), E r is increased up to ~ 1.12 times as γ is varying from 0 to 0.1. However, for temperature-dependent zeta potential, its enhanced sensitivity with temperature results in significant increment of streaming potential where E r increases up to ~ 1.3 and ~ 1.85 times for C ζ = 1 and C ζ = 4 respectively (figure 5a (v)). Now, the degree of confinement is incorporated within the parameter 0 κ , i.e. inverse of the thickness of EDL. The variation of streaming potential ratio (E r ) with 0 κ is shown in figure   5a (vi). Increasing 0 κ at one side increases the region of electroneutrality (ensuring a reduction of streaming current), while on the other hand, the strength of acting electrokinetic forces gets modulated and net E r is decided by these two factors. For constant zeta-potential case (i.e. C ζ = 0), E r first decreases with 0 κ up to 0 κ = 1.7, then increases slowly up to 0 κ = 3 beyond which it falls sharply resulting only ~ 1.12 times increment of E r visible at higher 0 κ ( 0 κ = 10).
Now we recall the boundary condition involved in evaluating the potential distribution, However, as described earlier in the variation of streaming potential ratio (E r ), increasing C µ also ensures faster migration of the ions caused by the imposed pressure gradient or the induced concentration gradient which also generates more electrokinetic flow in the reverse direction.
Therefore, the net throughput depends on the rate of lessening of viscous resistance as well as the enhanced induced streaming field. For lower strength of pressure gradient ( 0 0.01 p x − ∂ ∂ = ), an enhancement in the flow rate ratio (Q r ) up to ~ 3.35 times is observed as C µ is increased from 1 to 10. Increasing the strength of pressure gradient lowers the magnitude of the net flow rate with the dependence on C µ remaining qualitatively similar.
The alteration in the flow rate ratio (Q r ) with increasing C ε is highlighted in figure 5c (iv) where the reduction in the streaming potential owing to the temperature-sensitivity of electrical permittivity causes significant enhancement (~ 2.87 times) in the net throughput for purely thermally driven flow. However, the impact of C ε on Q r becomes less in presence of pressure gradient and interestingly, a cross-over between plots of

Effect of transverse temperature gradient
Here the primary driving force for flow actuation is the axially applied pressure gradient while the contribution from the imposed thermal gradient (∆T) is secondary through physical property alteration and introducing a permittivity-variation induced body force (this term is C y y Keeping χ constant, the velocity distribution for varying C µ is shown in figure 6a (ii).
Since C µ signifies the sensitivity of fluid viscosity with temperature, higher the value of C µ , higher is the reduction of the viscosity, lower is the resistance to drive the flow. Because of the strengthened advective current, streaming potential ratio (E r ) is increased slightly as C µ is increased from 1 to 5. For lower C µ , the velocity profile is parabolic in nature while at higher C µ , the enhanced sensitivity of fluid viscosity creates strong departure from the parabolic distribution with the maxima being shifted towards the hot region. Increasing 0 κ ensures lesser penetration of EDL in the bulk resulting in lowering the net streaming potential thereby enhancing the magnitude of the flow velocity to some extent as depicted by the dotted lines in figure 6a (ii). Figure 6a (iii) shows the transverse variation of the velocity field with increasing C ε .
Here, C ε indicates the sensitivity of permittivity with temperature. The contribution of C ε primarily comes through the alteration in the potential distribution upon increasing C ε and through the permittivity-induced forcing term in the fluid momentum transport. Here, streaming potential ratio (E r ) decreases (inset I of figure 6a (iii)) slightly with increasing C ε resulting small increment in the velocity magnitude. Similar to the effect of χ, here also increasing 0 κ from 1 to 10 makes the effect of C ε on the flow field inconsequential as evident from inset II (where zoomed view of maximum velocity is presented) of figure 6a (iii).
The variation of the velocity field in the y-direction for different ∆S̅ T is shown in figure   6a (iv) where ∆S̅ T means the difference in thermophoretic mobilities between cations and anions.
Increasing ∆S̅ T implies higher thermophoretic mobility of counter-ions than co-ions which leads to preferential migration of the counter-ions towards the cold region. This leads to an ionic redistribution resulting an asymmetry in the potential distribution within EDL. On observing the momentum equation, one can understand that the role of ∆S̅ T comes through the charge distribution alteration (via the modulated EDL thickness and the term which further perturbs the fluid advection motion. This in turn influences the advection current and the induced streaming field. Accordingly, the streaming potential ratio (E r ) increases from ~ 1.8 times to ~ 2.05 times as one increases ∆S̅ T from 0 to 1 following a linear dependence. The enhanced streaming field induces more backward flow thus lowering the magnitude of the flow velocity as observed in figure 6a (iv) for ∆S̅ T = 1.
As previously discussed in the case of axial thermal gradient, any perturbation to the flow field is strongly reflected in the associated dispersion characteristics because of its strong depedence on the non-uniformity of the flow velocity. By inspecting the velocity distributions demonstrated by figure 6a, one general observation can be made that the effect of most of the parameters (except C µ ) on flow field is noticeable only under higher confinement (i.e. at lower 0 κ ) and becomes negligible at higher 0 κ . Apart from 0 κ , C µ is another important parameter whose effect on flow field is far more significant compared to other parameters (C ε , ∆S̅ T ) not only by altering the magnitude of the flow velocity but also creating strong departure from the parabolic distribution. Accordingly, in figure 6b we have incorporated the variation of dispersion coefficient ratio (D ̅ eff ) with two parameters C µ and 0 κ . Similar to axial ∆T, drastic reduction in flow resistance with increasing C µ is reflected in figure 6b (i) where D ̅ eff is increased up to 3 times as C µ is varying from 1 to 10.
The reduced strength of electrokinetic forces with increasing 0 κ lowers the thermoelectric perturbation to the flow field and hence, dispersion coefficient ratio (D ̅ eff ) (at γ = 0.1) reduces from 1.19 to 1.18 as 0 κ is changed from 1 to 10 with alteration being suppressed beyond 0 κ = 4.4. Corresponding results for lower γ are similar with saturation occuring at higher 0 κ .  κ for different values of (i) C ε and (ii) ∆S̅ T respectively. Now, the effect of two other parameters C ε and ∆S̅ T are shown in figure 6c where the reduced streaming potential with increasing C ε results in increasing the dispersion coefficient.
However, C ε turns out to be less effective at lower thermal perturbation (γ = 0.05) and noticeable effect can only be observed at higher γ (figure 6c (i)). Similarly, higher volumetric suppression due to enhanced streaming potential with increasing ∆S̅ T lowers the dispersion coeffficient with ∆S̅ T becoming influential only at higher strength of pressure gradient (figure 6c (ii)).

Effect of fluid rheology
Now, the inclusion of rheological aspect of fluid on streaming potential is highlighted in figure   7a where the variation of the streaming potential ratio (E r ) is shown with Deborah number ( relaxation time (C λ ) of fluid with temperature. For a fixed value of C µ and C λ , introducing thermal gradient (γ = 0.1) results a reduction in the net streaming potential with De where an increment of ~ 1.27 times in E r is observed as opposed to ~ 2.7 times increment for γ = 0 (figure 7a (i)). Interestingly, a cross-over at De = 0.3 takes place between the graphs of γ = 0 and γ = 0.1.
Below this critical De, the magnitude of the streaming potential is higher for combined temperature gradient and pressure-driven flow and beyond De = 0.3, it falls below the streaming potential in isothermal condition in which case the strongly pronounced shear thinning effect with increasing De dictates the flow physics creating faster rise in streaming potential. Now, at lower C µ (denoting lower temperature-sensitivity of fluid viscosity), this two aforesaid counteracting factors are of comparable magnitude resulting in a slight increase in E r (from 1.12 times to 1.28 times) as De is varying from 0 to 1, visible in figure 7a (ii). With increasing C µ , its effect starts to become dominant over C λ and E r undergoes significant enhancement up to ~ 3.86 times (at C µ = 3) compared to the Newtonian fluid. Similarly, increasing C λ from 0.25 leads to faster reduction in fluid relaxation time denoting elevated elasticity-mediated disturbance to the axial separation between ions thus lowering the net streaming potential. As evident from figure 7a (iii), E r decreases from 1.46 to 0.79 as C λ is increased from 0.25 to 3. Beyond C λ = 1, E r starts to fall with De (from its reference value ~ 1.12) and becomes less than unity at higher De. This tells us that if C λ is high, (i.e. rapid reduction of fluid relaxation time with temperature) employing a Newtonian fluid is more advantageous instead of a viscoelastic fluid as far as streaming potential generation is concerned. Now, as shown in the inset of figure 7a (ii), the influence of C λ on streaming potential gets dampened at higher C µ (C µ = 2) where E r reduces at a lower rate from 2.75 times to 2.3 times (with respect to a Newtonian fluid) with increasing C λ .
The flow rate ratio (Q r ) variation with Deborah number (De) is highlighted in figure 7b for two varying factors C µ and C λ respectively. At one side, increasing C µ induces more streaming potential leading to net volumetric suppression due to reverse electrokinetic flow, on the other hand, viscous resistance in the flow decreases significantly. Hence, net throughput through the microchannel is determined by their relative strengths. For lower value of C µ , Q r remains almost unaffected with the variation of De (at C µ = 1). However, with increasing C µ , Q r increases sharply with an enhancement up to ~ 2.03 times can be noticed at C µ = 3 (figure 7b (i)).
(i) (ii) FIGURE 7b. Volume flow rate ratio (Q r ) with De for different (i) C µ and (ii) C λ respectively.
Similarly, increasing C λ has inverse effect on the flow rate as Q r decreases with De from 1.78 times to 1.18 times as C λ is changed from 0.25 to 3. Since the net streaming potential reduces beyond C λ = 1, reduced volumetric suppression makes Q r still higher than unity at C λ = 3.
This section has reached its culmination where the dependence of the dispersion coefficient ratio (D ̅ eff ) is shown with Deborah number (De) for two crucial parameters C µ and C λ .
Similar to the variation of Q r , D ̅ eff is also highly sensitive to the variation in fluid viscosity and relaxation time. Looking into the constitutive form of a viscoelastic fluid one can realize that the inherent non-linearity in stress-tensor terms gets amplified with increasing C µ . Physically, the impact of physical property alteration is reflected more in viscoelastic fluids compared to Newtonian fluids and accordingly, D ̅ eff increases up to 1.9 times at higher C µ (at C µ = 3). Also, (ii) (i) (iii) FIGURE 7c. Variation of dispersion coefficient ratio (D ̅ eff ) with Deborah number (De) for (i) different C µ and (ii) different C λ . (iii) Variation of D ̅ eff for specific combination of C µ and C λ. . changing C λ can result significant alteration in the dispersion coefficient where D ̅ eff exhibits an inverse dependence on Deborah number (De) at higher C λ , as observed in figure 7c (ii). Overall, the dispersion coefficient ratio (D ̅ eff ) undergoes a reduction from 1.53 to 1.1 as C λ is increased from 0.25 to 3.
From the previous figure, it is clearly evident that increasing C µ and decreasing C λ turns out to be favorable as far as the enhancement of dispersion is concerned. When we combine this

Conclusions
We have considered thermally-modulated elctrokinetic transport to realize significant enhancement in solute dispersion of a complex fluid through a microfluidic channel. Although several techniques in the past have been deployed towards modulating the uniform velocity profile of electrically actuated flows, improved hydrodynamic dispersion still remains unexplored. In this context, the present study shows that combining the interplay between thermal and electrical effects coupled with fluid rheology, one can achieve up to one order of magnitude enhancement of dispersion coefficient in a pressure-driven flow of an electrolyte. This is mediated by breaking the equilibrium of the charge distribution within the electrical double layer upon imposed thermal gradient, subsequent modulation in thermo-physical properties, and eventual alterations in the fluid motion. We believe that such complex coupling between thermal, electrical, hydro-dynamic and rheological parameters in small scales can be exploited to a benefit in the design of thermally-actuated micro and bio-microfluidic devices demanding improved hydrodynamic dispersion.