Sample dispersion in isotachophoresis with Poiseuille counterflow

Isotachophoresis (ITP) is a special mode of electrophoretic transport that has already been described by Kohlrausch 1 . It has found widespread applications in the analytical sciences (for an overview, see reference2 ) and is mainly used as a technique for the preconcentration and separation of samples. ITP relies on consecutively stacking a high mobility leading electrolyte (LE) and a low mobility trailing electrolyte (TE) in a capillary or channel as sketched in figure 1. Upon application of an electric field the ions arrange in the order of their electrophoretic mobility, forming a sharp transition zone between them that migrates along the capillary according to the electrophoretic velocity of the ions. Accordingly, there is a corresponding change in electric field across this transition zone, such that the high and low mobility species move at the same speed. A high mobility ion diffusing backwards across the transition zone into the low mobility region will experience this higher electric field and is thus transported back to its own region and vice versa for a low mobility ion diffusing into the region of high mobility ions. It is this balance between diffusion and electrophoretic migration that determines the width of the transition zone. Depending on whether the stacking occurs for the anions or the cations, anionic or cationic ITP is observed. Sample ions with an electrophoretic mobility in between the LE and TE ion mobilities are sandwiched between the two electrolytes and can form a zone of only a few micrometers width, depending on the total amount of sample.

Very early in the development of electrophoretic separation techniques researchers applied a counterflow opposing the electrophoretic sample migration to increase the time span available for the separation of analytes or to reduce the length of the separation passage 3 . While this principle was first applied to capillary electrophoresis, Everaerts et al. 4 introduced a counterflow to balance the sample migration in ITP experiments. Originally, pressure-driven counterflow was applied, but later electroosmotic counterflow was considered too 5 . In the past decades, counterflow-balanced ITP has found widespread applications, for example in the context of sample preconcentration for improving the sensitivity of capillary electrophoresis 6 ^-8 or to control the elution of analytes 9 .

Despite the importance of the method, it is rather poorly understood how the counterflow affects the analyte distribution in a sample focused by ITP. Urbánek et al.10 have noted in their experiments that pressure-driven counterflow increases sample dispersion. A few efforts have been made to compute the effects of counterflow on ITP. The one-dimensional models that have been developed 11 ^,12 may be suitable for flow profiles resembling a plug flow (such as electroosmotic flow in capillaries that are much wider than the Debye layer thickness), but do certainly not account for the complex convection-diffusion-electromigration phenomena occurring when a Poiseuille counterflow is applied. More successful one-dimensional models considering convective dispersion in ITP were developed with the aim to describe dispersion in a setting with a mismatch in electroosmotic flow Saville_1990 ; Ghosal_2003 ; 21 ; 22 .

In this article, we present the first computational fluid dynamics (CFD) study of sample dispersion occurring when ITP is balanced by a Poiseuille counterflow. For this purpose we have numerically solved the coupled Nernst-Planck and charge conservation equations using a finite-volume method. We analyze the most important dependencies of the width of the area-averaged sample distribution on the dimensionless parameters governing the problem. We show that in ITP with short sample zones and at large enough (but not unrealistic) Péclet numbers the flow considerably broadens the sample distribution compared to pure ITP without counterflow. Furthermore, we show that a one-dimensional Taylor dispersion model reproduces the data obtained with CFD reasonably well.

The migration of charged species in an electrolyte under an external electric field is governed by convection, diffusion and electromigration. The mass conservation of the ionic species leads to the Nernst-Planck equation,

Here $C_{i}$, $D_{i}$, $\mu_{i}$, $z_{i}$ denote, respectively, the concentration, diffusion coefficient, mobility and valency of the $i^{th}$ ionic species. Our model system comprises four ionic species: ions in the TE, sample and LE having a valency of the same sign and a common counter-ion with valency of opposite sign. The subscript $i$ denoting these species is chosen as $t$, $s$, $l$, and $0$, respectively. The mobility and diffusivity are related via the Nernst-Einstein relation $\mu_{i}=D_{i}F/RT$, where $F$ is the Faraday constant, $R$ the gas constant and $T$ the absolute temperature. V is the flow velocity field and $\textbf{E}=-\nabla\phi$ the electric field. The electric current density and charge density are defined as $\textbf{j}=F\sum_{i}{z_{i}N_{i}}$ and $\rho_{e}=F\sum_{i}{z_{i}C_{i}}$, where $\textbf{N}_{i}=-D_{i}\nabla C_{i}+C_{i}(\textbf{V}+z_{i}\mu_{i}\textbf{E})$ denotes net flux of the $i^{th}$ ionic species with $z_{i}=1$ for $i=t,s$ or $l$ and $z_{0}=-1$.

The charge conservation equation is

Under electro-neutrality13 , i.e. $\rho_{e}=0$, the time-dependent and convective terms in (2) vanish. By multiplying equation (1) by $Fz_{i}$ and taking the sum over all species the reduced charge conservation equation can then be written as

where the ionic conductivity is $\nu=(F\sum_{i}{\mu_{i}z_{i}^{2}C_{i}})$. The right-hand side term in (3) is known as the diffusion current and its contribution is insignificant at all locations, except for the transition zones where large concentration gradients occur. For microchannels with charged walls, electroneutrality is violated within a region of the order of the Debye screening length, $\lambda_{D}$, much smaller than the channel dimensions. In such cases equation (3) is retained in the channel bulk, while the effect of the charge layer is incorporated via an effective wall-velocity boundary condition taking into account the electroosmotic flowSchoenfeld_2009 . Since here we focus on sample dispersion by Poiseuille flow we assume uncharged walls.

In this article we consider a two-dimensional situation, i.e. the isotachophoretic transport and the Poiseuille flow occur between parallel plates. Translated to realistic microfluidic setups this means that very shallow channels with a width much larger than their depth are assumed. We impose a Poiseuille counterflow with an average speed equal and opposite to ITP migration, as

where $H$ is the depth of the channel and $U^{ITP}$ is the velocity of ITP transport for the case that no counterflow is applied. The situation is depicted in figure 1. Thus, on average the electromigration of sample ions is just compensated by convection, meaning that the average speed of the ITP transition zone (the “interface” between LE and TE) will be zero.

We non-dimensionalize the $X$- and $Y$-coordinates by the depth $H$ of the channel, i.e. set $x=X/H$ and $y=Y/H$. Further, all concentrations are non-dimensionalized by the bulk LE concentration, $C^{\infty}_{l}$, and the electric potential by $\phi_{0}=RT/F$. We choose the isotachophoretic velocity $U^{ITP}$ as velocity scale. Time is scaled by $\tau=H/U^{ITP}$, i.e. $\hat{t}=t/\tau$. The Nernst-Planck equation (1) can then be written in non-dimensionalized form as

where $c_{i}=C_{i}/C^{\infty}_{l}$ are the non-dimensional concentrations and the Péclet number is defined as $\mathrm{Pe}=U^{ITP}H/D_{s}$. The charge conservation equation (3) becomes

Under no counterflow, all the zones move with a constant isotachophoretic velocity $U^{ITP}=\mu_{t}E_{t}=\mu_{s}E_{s}=\mu_{l}E_{l}$, where $E_{i}$ ($i=t,s$ or $l$) are the local electric field strengths in the TE, sample or LE zone, respectively, given as

where $l_{i}=L_{i}/L$. $L_{i}$ ($i=t,s$ or $l$) is the portion of the channel filled by the respective electrolyte and $L$ is the total length of the channel. $E_{0}=\Phi/L$ is the average electric field due to the voltage drop $\Phi$ along the channel. Here we have assumed that the sample forms a distinct zone where its concentration remains constant and the overlap region between the zones are negligible. Since in ITP applications the channel is long compared to the section occupied by sample we take the limit $l_{s}\to 0$ in equation (7), which then is also valid for the electric field in the TE and LE even in the case of a dispersed sample zone. Also, we assume that the sample zone is situated in the middle of the channel, i.e. we choose $l_{l}=l_{t}=0.5$.

From the conservation of electric current along with electro-neutrality and the negligible diffusive flux at the far-field, a relationship between the TE and LE concentrations can be obtained as

where the suffix $\infty$ stands for the concentrations far away from the interfaces between the electrolytes 1 ^,14 . We will refer to these relations as the Kohlrausch conditions. Note that the second of these relations is only valid in the case where the sample concentration forms a distinct zone with constant concentration, i.e. when it is only marginally affected by diffusive and convective dispersion.

In solving the Nernst-Planck equation, the net ion fluxes through the channel walls are set to zero, i.e. $\textbf{N}_{i}\cdot\textbf{n}=0$ for $i=t,l$ or $s$, where n is unit outward normal. The electric potential is subjected to insulating boundary condition ($\nabla\phi\cdot\textbf{n}=0$) along the wall. Both the left and right boundaries of the computational domain are placed sufficiently far away from the sample zone such that the distribution of sample ions, $c_{s}$, is not affected by their presence. The concentrations far away from the transition zones are governed by the Kohlrausch condition (8).

Due to the presence of different electrolytes in different portions of the channel, a sharp change in conductivity will appear across the transition zones. If there is no bulk flow, the width of such a transition zone is of the order of 15

where $\Delta E$ is the change in electric field strength across the zone. Upon application of a voltage drop $\Phi$ across the channel, a locally uniform electric field appears sufficiently far away from each transition zone. This allows for calculation of the electric potential across the channel via charge conservation (3). In other words, owing to the electro-neutrality assumption we are spared the effort of solving the Poisson equation from which the electric potential is usually derived.

The electroneutrality assumption13 relies on the fact that its inevitable violation dictated by a changing electric field due to Gauss’s law, $\rho_{e}=\varepsilon_{e}\mathrm{div}\mathbf{E}$, is much smaller than the total amount of charges present, $\sim F|z_{i}|c_{i}$. For ITP applications the ratio between the two is typically of the order of $10^{-4}-10^{-5}$ (cf. supplementary material for Baier_2011 ; 22 ) which allows us to safely neglect the impact of $\rho_{e}/F$ on the concentrations. However, this is not the only place where the charge density enters, since in terms of Maxwell stresses it also contributes as a body force, $\sim\rho_{e}\mathbf{E}$, in the momentum equation for the fluid Lin_2004 ; 22 . This issue of Maxwell stresses in the context of ITP in a similar setting was discussed in Baier_2011 , where it is concluded that the body force term becomes important for large applied fields with very narrow transition zones and hence large field gradients, something that is not achieved with the strong convective dispersion considered here.

In order to make quantitative predictions using this model, the electrophoretic mobilities of all species, concentration of the leading electrolyte, voltage drop across the channel, and depth of the channel must be specified. Note that the TE concentration can be evaluated by using the Kohlrausch conditions (8). The electric field strengths at the ends of the channel can be adjusted via relation (7). We use $\varepsilon_{e}=78.5\times 8.85\cdot 10^{-12}\,\mathrm{C/(Vm)}$, $T=300\,\mathrm{K}$ and $\phi_{0}=25.9\,\mathrm{mV}$.

We developed a computer code to compute the governing equations based on the numerical algorithm as described below. The non-dimensional equations for ion transport and electric field are computed in a coupled manner through the finite volume method 17 . The computational domain is subdivided into a number of control volumes. When the electromigration in the Nernst-Planck equations dominates the electro-diffusional transport of ions, the transport equations show hyperbolic characteristics. Due to the hyperbolic nature of the ion transport equations, we adopt a higher-order upwind scheme QUICK18 (Quadratic Upwind Interpolation Convection Kinematics) to discretize the electromigration and convection terms in the ion transport equations. The diffusion flux at the control volume interfaces is estimated by a linear interpolation of variables between the two neighbors to either sides of the control volume interfaces. Details of the spatial discretisation scheme can be found in the supplementary material. An implicit first-order scheme is used to discretize the time derivatives present. At every time level we solve the ion transport equations iteratively, as the ion transport equations and the charge conservation equation are coupled. The iteration procedure starts with a guess for the electric potential at each cell center. At every iteration, the elliptic PDE for the charge conservation equation is integrated over each control volume through the finite volume method. The discretized equations are solved by a line-by-line iterative method along with the successive-over-relaxation (SOR) technique. The iterations are continued until the absolute difference between two successive iterations becomes smaller than the tolerance limit $10^{-6}$ for concentration as well as for the electric potential. The details on discretization of the governing equations are provided in the supplementary material.

A steady state solution is achieved by taking a sufficient number of time steps until the concentration distributions remain unchanged with time. The initial condition for ITP with counterflow is governed by the solution of the corresponding ideal ITP case (without convection). The solution for ideal ITP is also obtained based on the algorithm as described above and the parabolic velocity profile is imposed after a steady isotachophoretic zone is formed. We find that the steady-state state is archived after a short transition phase. The time evolution and formation of steady-state is illustrated in figure S4 of the supplementary material. However, in this paper all results presented correspond to the situation where a steady state has been reached.

As the variables within the transition zone between two electrolytes change more rapidly than elsewhere, a nonuniform grid spacing along the $x$-direction and uniform spacing along $y$-direction is chosen. The grid size $\delta x$ within the transition zones is relatively small and is increased with a constant increment as we move away from the transition zone. The smallest $\delta x$ is chosen such that the Courant-Friedrich-Levy (CFL) criterion is satisfied. For the present computations, $\delta x$ is taken to be $0.005$ for the finest grid and it is $0.01$ for coarse grids. Here $\delta y$ is taken to be $0.005$ and $\delta t=0.001$. A grid independency test and validation of our algorithm by comparing with the analytical solution for ideal ITP (Goet et al. 15 ) is presented in figure S3 in the supplementary material. The efficiency of the present numerical code in resolving the transient zones for peak and plateau-mode ideal ITP in steady-state is illustrated in figure S5 of the supplementary material.

To analyze the dispersion effect on the sample when a counterflow is applied, we consider a one-dimensional model based on Taylor-Aris dispersion 19 ^,20 . Taylor’s original work was concerned with dispersion solely due to the pressure driven flow through a channel. In that case the mechanism at work is the interplay of convective dispersion and lateral diffusion. Convective dispersion alone would result in a dispersion linear with time, but its effect is limited by the lateral diffusion over the cross section of the channel. Combined, these phenomena result in a dispersion proportional to the square root of time, such that the area-averaged sample distribution can be considered to spread under the influence of an effective axial diffusion. For two reasons it is a priori unclear whether such a Taylor-Aris dispersion model can be applied to the situation considered in this article. First, electromigration appears as an additional transport process that results in a sharpening of the sample distribution. Secondly, an assumption implicit to the Taylor-Aris model is the infinite time limit, i.e. the model is valid only on time scales large compared to the time a molecule takes to diffusively sample the channel cross section. In turn, this means that it is usually applicable only to comparatively broad (in axial direction) sample distributions with $\sigma/H\gg\mathrm{Pe}$, where $\sigma$ is the length scale of a transition zone. For this reason one would a priori not expect that a description based on an effective axial diffusivity is applicable. Nevertheless, such models have been successfully used previously in similar situations Saville_1990 ^,21 ^,22 . Therefore we will investigate the applicability of this model to the case considered here.

The core of the Taylor-Aris dispersion model is the effective axial diffusion due to convection. In particular, the diffusion coefficient is replaced by a term dependent on the Péclet number squared

where $i=t,s$ or $l$. For a parallel-plates channel the value of the constant $\beta$ is $1/210$. The dispersion model is based on the area-averaged ion concentrations

Likewise, we write area-averaged transport equations for the ion concentrations in which the diffusivity is replaced by the effective axial diffusivity according to Taylor and Aris

The axial electric field is determined by solving the corresponding one-dimensional charge conservation equation assuming electro-neutrality throughout

where $\overline{\nu}$ is the area-averaged conductivity. The equation for the electric field is coupled with the transport equations (12). The numerical solution of the species transport equations above is achieved by employing a higher-order accurate upwind differencing scheme. The equation for the electric field is solved along with the ion transport equations in a coupled manner as outlined before.

We consider the LE and common ion diffusivities as fixed values throughout this study, setting $D_{l}=D_{0}=7.0\times 10^{-10}\,\mathrm{m^{2}/s}$. The diffusivities of the TE and sample ions are prescribed by specifying the two non-dimensional parameters $k_{1}=D_{l}/D_{t}$ and $k_{2}=D_{l}/D_{s}$. The corresponding mobilities are related to the diffusivities via the Nernst-Einstein relation, $\mu_{i}=D_{i}F/RT$. Since the sample mobility lies between the mobilities of TE and LE, these parameters can be varied subject to the restriction $k_{1}>k_{2}>1$. The scale of the applied electric field is set to $E_{0}=10^{5}\,\mathrm{V/m}$ from which the field in the respective buffers can be inferred using equation (7). The bulk LE concentration is taken as $C^{\infty}_{l}=10^{-3}\,\mathrm{M}$, so that the Debye layer thickness ($\sim 10\,\mathrm{nm}$) is much smaller than the depth of the channel considered here ($>5\,\mathrm{\mu m}$). The bulk TE concentration can be inferred from the Kohlrausch condition (8), which in our example becomes $C^{\infty}_{t}/C^{\infty}_{l}=2/(1+k_{1})$. A corresponding upper bound can be obtained for the maximum of the sample concentration, which due to dispersion will usually be much smaller in the cases considered here, however.

In ITP experiments one considers two modes of operation, peak and plateau mode, depending on the amount of sample present in the system. ”Plateau mode” refers to the case where the sample-zone is wide compared to the transition zones, i.e. the sample concentration distribution forms a constant plateau with blurred boundaries towards LE and TE. ”Peak mode” refers to the other extreme of a very short sample zone, where the two transition zones at both sides of the sample overlap or more precisely, when the sample is entirely within the diffused interface between LE and TE. In this situation the influence of the analytes on the overall conductivity is typically negligible. Everything else being equal, the difference between peak and plateau mode lies in the amount of sample present. In the present case, owing to the strong additional dispersion due to convection, a larger amount of sample will still lead to a pronounced peak instead of a plateau, as will be elaborated below. In order to distinguish this regime from peak mode as observed without convective dispersion, we will call this situation “dispersed plateau mode”. This is the regime that we will mainly, but not exclusively, consider presently. It is different from plateau mode in that the area averaged sample distribution shows a pronounced peak and differs from peak mode in that the sample strongly affects the overall conductivity.

In order to obtain comparable results for the distortion of the sample region due to the applied Poiseuille flow, we fix the area-averaged amount of sample

present in the channel. In that way the area-averaged distribution without convective dispersion is independent of the channel height. Note that the exact form of the distribution will depend on the choice of other parameters such as $k_{1}$ and $k_{2}$ due to the different concentrations dictated by the Kohlrausch conditions (8) or the width of transition zones (9). Unless stated otherwise we will mainly use the values of ${\cal C}_{s,0}=40\,\mathrm{\mu mol/m^{2}}=40\,\mathrm{\mu m}\cdot C^{\infty}_{l}$ and ${\cal C}_{s,0}=4\,\mathrm{\mu m}\cdot C^{\infty}_{l}$ throughout; the former value leads to a distinct plateau without convective dispersion while it represents dispersed plateau mode for wide enough channels, while the latter corresponds to a sample distribution in the form of a Gaussian peak without convective dispersion, cf. figure S5 in the supplementary material.

We remark that without convective dispersion the width of the sample zone in plateau mode can easily be inferred from the Kohlrausch condition (8) and the amount of sample present, neglecting the width (9) of the diffusive zones, to $L_{s}\approx{\cal C}_{s}/C_{s}^{\infty}={\cal C}_{s}/C_{l}^{\infty}\cdot(1+k_{2})/2$, which gives a baseline for the width of the sample distribution.

In figure 2 the corresponding sample concentration profiles are shown for two different channel heights. As mentioned, all results correspond to a situation where a steady state has been reached. Figure 3 shows the corresponding area-averaged sample concentrations, both derived from the 2D model and from the 1D Taylor-Aris dispersion model (12) with effective diffusivity (see supporting information for a combined graph of the TE, sample and LE distributions in this situation). From these figures it becomes apparent that whether a plateau or peak is obtained now becomes dependent on the channel height. For the shallow channel a plateau is observable in the area-averaged concentrations, contrasting the distinct peak observed for the wider channel, characteristic for the “dispersed plateau mode”. Since in both cases the ITP velocity is the same, the larger dispersion is a direct consequence of the larger Péclet number for the wider channel. We also note that the 1D model based on the effective Taylor-Aris diffusion captures the area-averaged distribution quite well. This is particularly true for the narrower channel. In the wider channel the larger discrepancy is due to the slight focusing of the sample in the center region of the channel, which effectively results in a narrower peak than determined from the 1D model. We will further investigate this below by varying the total sample concentration present in the channel.

Before doing so, we investigate the area-averaged distributions obtained with the two models for different ratios of the diffusivities, cf. figure 4 (a) when the sample amount is ${\cal C}_{s}/C^{\infty}_{l}=40\,\mathrm{\mu m}$. Corresponding results for the reduced sample amount i.e., ${\cal C}_{s}/C^{\infty}_{l}=4\,\mathrm{\mu m}$ is shown in figure 4 (b). In the cases where the sample diffusivity is very close to the diffusivity for either LE or TE one expects a very wide transition zone between the sample and the respective electrolyte, since hardly any electromigrative sharpening is present. On the other hand, for the transition towards the electrolyte with markedly different diffusivity a much narrower zone is expected. This is indeed seen in figure 4 (a) and (b). We remark that the three situations depicted correspond to different Péclet numbers, since the parameter varied is the sample diffusivity. Nevertheless, the total width of the sample zone is affected more strongly by the long tails towards the electrolyte with similar diffusivity than by the Péclet number. We will further investigate this fact below. Also note that again the Taylor-Aris model quite successfully captures the area-averaged concentration profiles, in particular towards the electrolyte with similar diffusivity. The sharper profile towards the electrolyte with dissimilar diffusivity in the 2D compared with the 1D case is again due to the more detailed 2D structure with a slight focusing in the channel center emerging in this model. Comparing figure 4 (a) and (b) we find that in the dispersed plateau mode the profiles for the sample distribution are largely independent of its amount.

To further investigate the validity of the 1D model, the total amount of sample present in the channel is varied. Figure 5 shows the 2D profiles obtained for amounts ranging from $1/2$ to $5$ times ${\cal C}_{s,0}$, the amount considered so far. For the smallest amount of sample, one again observes the focusing of the sample in the central region of the channel. Contrasting that, for the largest amount of sample present, the situation may be considered as plateau mode. The corresponding area-averaged concentration profiles are shown in figure 6, both for the 2D and 1D cases. Again, the 1D model is surprisingly effective in predicting the profiles, with the quality of the predictions deteriorating for small sample amounts due to the strong 2D structure of the profile, leading to a sharper focusing than predicted by the 1D model.

In order to quantify the dispersion, we define the width of the sample zone as the second moment of the sample distribution

where the center of the distribution is defined as

Additionally, we also measure the skewness of the distribution through the third normalised moment as

to further quantify the similarity between results for the sample distribution obtained within the 2D model and the 1D Taylor-Aris dispersion model. The value of skewness may be positive, negative or zero depending on whether the sample electrolyte is skewed to the right (positive skew), to the left (negative skew) or symmetric (zero skew).

The stage is now set for a more systematic study of the influence of the individual parameters on the form of the sample distribution. We use the width of the sample zone, $\sigma$, and its skewness, $\gamma_{1}$, in figures 8, 7 and 9 to investigate the dependence of the dispersion on diffusivity ratios of the individual electrolytes, the amount of sample present and the channel depth, respectively.

In particular, figure 7 shows the distribution’s width and skewness for several values of $k_{1}$ when $k_{2}$ is varied in the range $1<k_{2}<k_{1}$. In order to compare the results we have plotted them using the rescaled variable $(k_{2}-1)/(k_{1}-1)$ so that the curves span the same range in parameter space. In this figure the remark already made in the discussion of figure 4 is quantified: the sample dispersion is largest in the cases where the sample mobility is close to the mobility of one of the other electrolytes. This is reflected both in the increase of width, $\sigma$, as well as its skewness, $\gamma_{1}$, for values of $k_{2}$ close to either of the extreme cases 1 or $k_{1}$. The sign of the skewness reflects the fact that for $k_{2}$ close to 1, i.e. the diffusivity is very close to the LE diffusivity, the sample distribution has a long tail reaching into the LE, and similarly when the diffusivity is close to the TE diffusivity, i.e. $k_{2}$ close to $k_{1}$, the tail reaches to the left into the LE. We find the skewness to be almost zero, i.e. a symmetric sample distribution, when its mobility, $\mu_{s}$, is close to the harmonic average of $\mu_{t}$ and $\mu_{l}$. We note that the individual curves corresponding to different values of $k_{1}$ almost coincide showing that the dependence on $k_{1}$ is relatively weak when considering the rescaled $k_{2}$-values.

To further investigate the universality of the distribution width and skewness curves, we show in figure 8 the corresponding curves for different amounts of sample present in the system. In particular, we show results for a sample amount of ${\cal C}_{s}/C^{\infty}_{l}=40\,\mathrm{\mu m}$ and ${\cal C}_{s}/C^{\infty}_{l}=4\,\mathrm{\mu m}$ and find that both the width as well as the skewness of the distributions do not significantly change. As noted earlier, these two situations corresponds to plateau and peak mode in a situation of ITP without additional convective dispersion. This points towards the fact that the form of the distribution in dispersed plateau mode is not strongly dependent on the amount of sample present, the difference being mainly a scale factor proportional to the total amount, as already hinted at in the distributions shown in figure 4.

Both in figure 7 and 8 we show data obtained within the 2D calculation as well as for the effective 1D Taylor-Aris dispersion model. As can be seen, there is good agreement between the two calculations corroborating the use of the simplified model for predicting area-averaged sample distributions in cases with convective dispersion. Nevertheless, it is important to keep in mind that the actual sample distribution within the channel will usually be richer in features in the 2D case, in particular we mention the focusing of the sample in the channel center observed for example in figures 2 and 5.

Another interesting aspect is how the width of the sample zone changes when varying the channel depth $H$, and hence the Péclet number, as displayed in figure 9. From the effective diffusivity, eq. (10), one expects a sample zone width that scales as $\sim\mathrm{Pe}^{2}$ in dispersed plateau mode. For large enough $H$, such that the dispersed plateau mode is indeed observed, this is reproduced well by the 1D model, as it should be. However, in the region shown this limit is not fully reached yet and only for values of $H$ larger than $30\;\mathrm{\mu m}$ the points follow $\sigma\sim H^{2}$ for the 1D results. As can be seen from the figure, the 2D model predicts less sample dispersion, a fact that can again be attributed to the focusing of sample in the channel center. Figure 9 (b) shows that the skewness of the distribution is correspondingly amplified by the increasing channel height.

We have investigated sample dispersion in situations in which electromigration in isotachophoresis is balanced by a Poiseuille counterflow. The focus was put on small sample amounts, i.e. situations when the area averaged sample distribution shows a peak. It is found that at large Péclet numbers, the sample dispersion is significantly increased compared to a situation without counterflow. Since the amount of sample is such that without the counterflow present one would typically obtain a short sample zone in plateau mode, we have termed this regime “dispersed plateau mode”. One of the findings of this study is summarized in figures 7 and 8 which show the dispersion of the sample zone as a function of the sample ion mobility. In particular, the diagram shows in which window the mobility of the sample can be chosen relative to the mobility of the surrounding electrolytes before significant broadening of the sample zone is observed. In particular, we find that the area averaged sample distributions are only weakly dependent on $k_{1}=\mu_{l}/\mu_{t}$ but vary strongly with $k_{2}=\mu_{l}/\mu_{s}$ (or more precisely with the rescaled value $(k_{2}-1)/(k_{1}-1)$). These results are also observed to be relatively insensitive to the amount of sample present in the system, cf. figure 8, with the general form of the distribution remaining the same but the peak–height scaling with the total amount. These curves can thus be used to asses the degree of dispersion and skewness expected in an experiment. Similar results hold for the case in which it is desired to focus two sample zones stacked behind each other without the risk of too strong intermixing.

Another key result of this study is that a 1D area-averaged model for the sample distribution based on Taylor-Aris dispersion agrees well with the more detailed 2D model in many situations. In view of the original derivation of the effective diffusion coefficient, the quality of the agreement is surprising, since the complex interplay of convection, electromigration and diffusion that shapes a transition zone in ITP is not accounted for. Since the computational cost of the 1D model is much less than that of the 2D model, this opens the door for fast simulations of ITP processes that are superimposed by convection. This could be of considerable relevance for a number of processes involving sample preconcentration and separation. However, it is important to stress that the detailed spatial distribution of the sample ions generally is more complex than the simplified 1D picture suggests. In particular we observe in the 2D calculations that the sample has a tendency to become concentrated in the channel center. This will have an impact in situations where different ionic species are made to react within an ITP experiment, as for example described in Goet et al.Goet_2009 and Bercovici et al.Bercovici_2012 , when long reaction times are attained via a Poiseulle counterflow.