## Abstract

Electrophoresis techniques are characterized by concentration disturbances (or waves) propagating under the effect of an electric field. These techniques are usually performed in microchannels where surface conduction through the electric double layer (EDL) at channel walls is negligible compared with bulk conduction. However, when electrophoresis techniques are integrated in nanochannels, shallow microchannels or charged porous media, surface conduction can alter bulk electrophoretic transport. The existing mathematical models for electrophoretic transport in multi-species electrolytes do not account for the competing effects of surface and bulk conduction. We present a mathematical model of multi-species electrophoretic transport incorporating the effects of surface conduction on bulk ion-transport and provide a methodology to derive analytical solutions using the method of characteristics. Based on the analytical solutions, we elucidate the propagation of nonlinear concentration waves, such as shock and rarefaction waves, and provide the necessary and sufficient conditions for their existence. Our results show that the presence of surface conduction alters the propagation speed of nonlinear concentration waves and the composition of various zones. Importantly, we highlight the role of surface conduction in formation of additional shock and rarefaction waves which are otherwise not present in conventional electrophoresis.

## 1. Introduction

Electrophoresis techniques such as capillary zone electrophoresis (CZE), field-amplified sample stacking (FASS) and isotachophoresis (ITP) are used widely for separation and in some cases, preconcentration of ionic species from sample mixtures [1,2]. Separation and preconcentration in electrophoresis techniques are based on differential electrophoretic mobilities of ionic species, that is, speed of ions in the electrolyte solution per unit applied electric field [3]. Electrophoresis techniques are typically characterized by linear and nonlinear concentration waves (disturbances in species concentration) propagating under the effect of electric field. Linear concentration waves are observed in CZE where analyte concentrations are small compared with those of background electrolyte. Therefore, propagation of analyte zones in CZE does not alter the local electrical conductivity and electric field, resulting in linear waves. For finite amplitude concentration disturbances such as those observed in ITP and during electromigration dispersion in CZE, nonlinearity in electromigration results in formation of shock and rarefaction waves [3–6]. The nonlinearity results due to the dependence of electromigration flux on local electric field, which in turn depends on the species concentrations through Ohm's law [5]. In CZE, nonlinear waves are undesirable as they result in unwanted dispersion of analyte zones [5,6] whereas in ITP, self-sharpening shock waves enable simultaneous preconcentration and separation of analytes between zones of high-mobility (LE) ions and low-mobility trailing electrolyte (TE) ions (figure 1*a*) [1,2]. The self-sharpening nature of shock waves separating analyte zones in ITP counter molecular diffusion of analytes and can enable over million-fold preconcentration of analytes [7].

Typically, electrophoresis experiments are performed in microchannels or microcapillaries with characteristic internal dimension *h* of ^{−1}) [1,2]. For electrolyte ionic strength of *ρ*_{w}/(*zcFh*) [9,10]. Here, *z* and *c* are, respectively, the characteristic valence and concentration of ions in bulk solution, *F* is the Faraday's constant and *ρ*_{w} is the surface charge density. With surface charge density of ^{−2}), in conventional electrophoresis experiments the contribution of surface conduction in the EDL to overall conduction through the electrolyte is negligible

Recent technological advancements in nanofabrication have led to integration of electrophoresis techniques, such as FASS and ITP, in shallow microchannels and nanochannels with characteristic internal dimensions of

Propagation of nonlinear concentration waves during electrophoresis in binary electrolytes due to the competing effects of surface and bulk conduction has been modelled by Mani *et al.* [9,10] and experimentally validated by Zangle *et al.* [11]. However, the model of Mani *et al.* [9,10] is not applicable for analysing multi-species electrophoretic processes as it does not account for nonlinearities induced by electromigration of multiple species. Note that, multi-species electromigration is nonlinear even in the absence of surface conduction. On the other hand, existing mathematical models for multi-species electrophoretic systems [4,15] are not applicable for electrophoresis systems where surface conduction competes with bulk conduction. This is because, the existing models assume non-zero electrophoretic mobilities and hence they cannot account for surface charge which results from immobile ionic species at the channel surface. To our knowledge, no mathematical model exists which accounts for the effects of surface conduction on bulk electrophoretic transport in multi-species electrolyte systems.

We here present a novel mathematical model for electrophoretic transport in multi-species electrolyte systems considering the effects of both surface and bulk conduction. Our model is applicable for electrophoretic systems consisting of arbitrarily large number of co-ions, having same charge polarity as the surface charge, and a single counter-ion. We also provide a methodology for obtaining analytical solutions to the governing equations of model. To this end, we reformulate the governing equations for electrophoretic transport in terms of Riemann invariants and solve them using the method of characteristics. Based on the analytical solutions to the governing equations, we highlight the effects of surface conduction on propagation of nonlinear waves in electrophoresis. In particular, we show that unlike in conventional electrophoresis, the zone concentrations do not obey the Kohlrausch regulating function [16] when surface conduction is prominent. Using the example of ITP, we show that surface conduction alters the zone compositions and propagation speed of concentration waves. Moreover, surface conduction results in formation of new concentration shock or rarefaction waves which are otherwise not present in conventional electrophoresis techniques.

## 2. Mathematical modelling

The species transport equations which govern the transport of ions due to advection, electromigration and diffusion in an electrolyte solution consisting of (*N*+1) ionic species are given by
*c*_{i}, *μ*_{i} and *D*_{i}, respectively, denote the concentration, electrophoretic mobility and molecular diffusivity of ionic species *i* in the bulk solution. In this equation, **u** is the bulk fluid velocity and **E** the applied electric field. Here we have assumed that all species in the electrolyte solution are fully ionized and we have neglected the effect of ionic strength on electrophoretic mobilities of ions [17]. Therefore, we treat electrophoretic mobilities and molecular diffusivities of all ionic species as constants.

The nature of concentration waves in electrophoresis and their propagation are primarily governed by electromigration. Molecular diffusion only acts to diffuse the gradients, while bulk fluid flow advects all the waves at same velocity. In a majority of electrophoresis experiments, particularly those performed in microfluidic chips, applied electric field is high (^{4} V m^{−1})) and experiment time is short (*A* and internal perimeter *P*. As the concentration gradients in electrophoresis are primarily in the axial direction, the governing equations (2.1), can be simplified using the above assumptions to
*x* denotes the axial coordinate and *E* the axial electric field. The spatio-temporal evolution of species concentrations depends on the local electric field in the electrolyte solution. The local electric field in bulk solution and EDL can be modelled using the Poisson's equation of electrostatics. However, the coupled set of equations consisting of the system of conservation laws (2.2) and the Poisson's equation is stiff. To avoid solving a stiff system, we here use the approximation proposed by Mani *et al.* [9], wherein EDL is modelled as a delta distribution of counter-ions at the channel walls having opposite charge polarity compared with the surface charge. The counter-ions in the EDL screen the surface charge and hence the bulk solution outside the EDL is assumed to be electroneutral. This approximation has been referred to as the ‘Leaky Membrane Model’ by Yaroshchuk [18] and Dydek & Bazant [19] in their analysis of concentration polarization in non-ideally perm-selective (‘leaky’) charged porous media due to electrophoretic transport in binary electrolytes. This assumption is illustrated in figure 1*a*,*b*, where the EDL is approximated as a zone of vanishingly small and uniform thickness *δ*. The EDL and bulk solution can be considered as two parallel conductive layers through which current can flow as shown in the equivalent circuit diagram in figure 1*c*. That is, the total current *I*_{T} through the electrolyte solution is the sum of currents through the EDL *I*_{EDL} and the bulk solution *I*_{B}, respectively. As the bulk solution is electroneutral, the electric field in bulk *E* in the absence of diffusive current is given by Ohm's law [3],
*σ*_{B} is the local electrical conductivity of the bulk solution, *F* the Faraday's constant and (*A*−*Pδ*) the cross-sectional area of channel excluding the EDL. Next, using the electroneutrality assumption *σ*_{B} can be simplified in terms of concentrations of only *N* species,

Equation (2.3) relates the local electric field (*E*) and the current in bulk solution (*I*_{B}). However, in electrophoresis experiments only total current, *I*_{T}=*I*_{B}+*I*_{EDL}, can be measured and often maintained constant. Therefore, to find a relationship between electric field and total current, we consider the contribution of surface conduction to the total current. We model the EDL as a narrow zone of thickness *δ* consisting of only counter-ions; we take (*N*+1)th species as the lone counter-ionic species and species 1 to *N* as co-ions. Noting that the counter-ions in EDL completely screen the wall surface charge density *ρ*_{w}, the counter-ion concentration in the EDL is given by *c*′_{N+1}=−*ρ*_{w}/(*z*_{N+1}*Fδ*). Therefore, the species transport equation for the counter-ions, equation (2.2), in the EDL can be written as
*E* and *I*_{EDL} which resembles the Ohm's law
*Pδ* is the cross-sectional area of EDL. Next, using equations (2.3) and (2.6), and noting that *I*_{T}=*I*_{B}+*I*_{EDL}, we obtain a relation between electric field (*E*) and total current (*I*_{T}),

Finally, using this expression for electric field in the equation (2.2) along with the relations for *σ*_{B} (equation (2.4)) and *σ*_{EDL} (equation (2.5)), we obtain a coupled set of *N* equations for electrophoretic transport of co-ions in the presence of surface conduction,
*A*−*Pδ*)≈*A* and introducing a new variable, *u*_{i}=*z*_{i}(*μ*_{i}−*μ*_{N+1})*Fc*_{i}, the governing equations can be written in a compact form given by
*α*_{i}=*μ*_{i}*I*_{T}/*A* and *k*=−*μ*_{N+1}*ρ*_{w}*P*/*A*. For our calculations, we assume that *α*_{i}>0, which can be ensured by reversing the *x*-direction in case *μ*_{i}*I*_{T}<0. Note that *u*_{i} are non-negative scalar quantities and can be converted back to species concentrations *c*_{i} using *c*_{i}=*u*_{i}/*β*_{i}, where *β*_{i}=*z*_{i}(*μ*_{i}−*μ*_{N+1})*F*. We note that *σ*_{B} and *k* always have non-negative values. Therefore, the presence of surface conduction in equation (2.9) reduces the propagation speed of concentration waves for same magnitude of total current. We will discuss the effect of surface conduction on wave propagation speed in §4. In the absence of wall surface charge (*ρ*_{w}=0) or large *A*/*P* ratio, the system of conservation laws (2.9) reduces to those for electrophoretic transport in large microchannels where surface conduction is negligible [4].

To derive the governing equations (2.9), we have made several simplifying assumptions which warrant further discussion. Owing to the assumption of thin EDLs, our model is applicable when channel thickness is sufficiently small to ensure surface conduction competes with bulk conduction, but not so small that EDL thickness is comparable with channel thickness. For example, for a circular channel with diameter of 1 μm and surface charge density of 0.3 C m^{−2} [20], the contribution of surface conduction, *k*=−*μ*_{N+1}*ρ*_{w}*P*/*A*, to the electromigration flux in equation (2.9) is comparable with that of bulk conduction for species concentrations up to order 1–10 mM. As such species concentrations are typical of many electrophoresis techniques, effects of surface conduction must be taken into account. Moreover, the EDL thickness at these concentrations is order 1–10 nm, which justifies the thin EDL assumption. In our model, we have accounted for additional conductivity associated with EDL by considering EDL and bulk solution as two resistors in parallel (figure 1*c*). A rigorous derivation of such homogenized ion-transport equations for binary electrolyte systems can be found in [21]. We also note that our ‘leaky’ membrane-model based on the assumption of uniform background charge is similar to the classical Teorell–Meyer–Sievers model [22,23] of ion-transport through ion-exchange membranes. However, the Teorell–Meyer–Sievers model was not used to predict ion-concentration shock waves prior to the work of Mani *et al.* [9]. In our analysis, we have also neglected the dispersive effects of EOF. This is justified because, as predicted by Dydek *et al.* [24] and demonstrated experimentally by Nam *et al.* [25], the hydrodynamic dispersion and associated overlimiting current due to EOF is negligible for narrow microchannels with internal dimensions of order 1 μm. Moreover, the dispersive effects of EOF are important primarily in the ion-depleted region next to an electrode or at the junction of microchannel and nanochannel [24,25]. However, we note that the surface conduction due to electro-osmotic convection can be comparable with surface conduction due to electromigration. Recently, Nielsen & Bruus [26] showed that for a binary electrolyte in a circular microchannel with thin EDLs, electro-osmotic surface convection can enhance surface conduction by a factor of (1+2*Pe*), where *Pe* is a normalization Peclet number of order 0.2. As the effect of EOF is to only amplify surface conduction, our model yields qualitatively correct predictions. Moreover, for quantitative comparison of model predictions with experimental data the surface conductivity *k* in equation (2.9) can be used as a single fitting parameter, as done by Zangle *et al.* [11], which can account for surface conduction due to both electromigration and EOF.

### (a) Reduction to Riemann invariants

The governing system of equations (2.9) bear a striking similarity with those describing multi-component chromatography [3,27,28]. Analysis of such systems was pioneered by Rhee *et al.* [27] and Kuznetsov [28] using the method of characteristics, and we here follow their approach to construct analytical solutions to the system of conservation laws (2.9). We begin by reformulating the set of equations (2.9) in terms of Riemann invariants so as to construct analytical solutions using the the method of characteristics. We assume that the co-ionic species are numbered in order of the magnitude of their unique electrophoretic mobilities, i.e. |*μ*_{1}|<|*μ*_{2}|<⋯<|*μ*_{N}|. We first show that the system of conservation laws (2.9) is strictly hyperbolic. To this end, we reformulate equation (2.9) as
*δ*_{ij} is the Kronecker delta. We note that here repeated indices do not imply Einstein summation. The eigenvalues *ξ* of the Jacobian matrix *A*_{ij} can be obtained by expanding the determinant det(*A*_{ij}−*ξδ*_{ij})=0, which yields the following characteristic equation

As the electrophoretic mobilities of co-ionic species are unique and ordered (*α*_{1}<*α*_{2}⋯<*α*_{N}), the roots of equation (2.11) *R*_{1}≤*α*_{1} and *α*_{m−1}≤*R*_{m}≤*α*_{m}, *m*=2,…,*N*; there are no repeated roots. Therefore, the eigenvalues *ξ*_{m}=*R*_{m}/(*σ*_{B}+*k*) of the Jacobian matrix *A*_{ij} are distinct and ordered as *ξ*_{1}<*ξ*_{2}<⋯<*ξ*_{N}. This implies that the system of conservation laws (2.9) is strictly hyperbolic [29]. In other words, we can reformulate the system of conservation laws (2.9) in terms of Riemann invariants. We note that for a certain class of electrolytes consisting of weak acids and bases, called oscillating electrolytes [30], the eigenvalues of the Jacobian of flux vector are complex valued. Such electrolytes exhibit unstable electrophoretic transport of ions. However, for fully ionized acids and bases, as we have assumed here, the eigenvalues are always real valued and distinct as shown above.

Interestingly, for the current system of conservation laws ∂*R*_{m}/∂*u*_{i} is the left eigenvector of the Jacobian matrix *A*_{ij} corresponding to the eigenvalue *ξ*_{m} for *m*=1,…,*N*. That is,
*u*_{i} and noting that ∂*R*_{m}/∂*u*_{i} is collinear with the *m*th left eigenvector of *A*_{ij}, *R*_{m}/∂*u*_{i} and using (2.12), we obtain
*R*_{m} remains constant along the characteristic curves given by *dx*/*dt*=*ξ*_{m}.

To completely write the set of equations (2.13) in terms of the Riemann invariants, we need to express *ξ*_{m}=*R*_{m}/(*σ*_{B}+*k*) in terms of the Riemann invariants alone. This necessitates an explicit relation between *σ*_{B}+*k* and the Riemann invariants. To obtain this relation, we expand the characteristic equation (2.11) to get
*σ*_{B}+*k* and Riemann invariants,

We note that these set of equations differ from the case of conventional electrophoresis where surface conduction is negligible. In the absence of surface conduction, one eigenvalue is zero and the Riemann invariant corresponding to zero eigenvalue is the Kohlrausch regulating function [16], whereas in the presence of surface conduction none of the eigenvalues are zero. Therefore, unlike in conventional electrophoresis, the zone concentrations do not obey the Kohlrausch regulating function when surface conduction is prominent.

### (b) Relation between species concentration and Riemann invariants

To solve the governing equations (2.16) we need to convert the initial conditions in terms of species concentrations to Riemann invariants and then convert the resulting solution in terms of Riemann invariants back to species concentrations. This necessitates an explicit relation between the species concentrations *c*_{i} and the Riemann invariants. As the Riemann invariants *N* roots of polynomial equation *L*(*u*,*R*)=0 (2.11) of order *N*, the polynomial *L*(*u*,*R*) can be written as
*σ*_{B}+*k*) has been chosen to match the coefficient of *R*^{N} in equation (2.11). Finally, substituting *R*=*α*_{i} in the above equation, we obtain an explicit expression for species concentration *c*_{i} in terms of the Riemann invariants,

## 3. Analytical solution using the method of characteristics

The set of conservation laws (2.9) (or (2.16)) are a system of quasi-linear hyperbolic equations. This suggests that multi-species electrophoretic transport in the presence of surface conduction can exhibit nonlinear concentration waves such as shock and rarefaction waves. Moreover, the system of conservation laws is strictly hyperbolic which precludes the existence of undercompressive and overcompressive shocks. For the current system of conservation laws (2.16), the conditions for genuine nonlinearity [29] hold for all characteristic fields. That is,
*A*_{ij} corresponding to the *m*th eigenvalue *ξ*_{m}. The conditions of genuinely nonlinearity are analogous to the convexity of flux in scalar equations, and they ensure that for a Riemann problem the characteristics are either compressing or expanding. We can therefore construct solutions for the Riemann problem using classical shock and expansion waves.

To illustrate the effect of surface conduction on the propagation of nonlinear concentration waves in electrophoresis, we first present a generalized solution to the Riemann problem with piecewise constant concentrations and a single discontinuity in the initial condition. Based on the generalized solution, we illustrate the nature of nonlinear concentration waves that form when surface conduction competes with bulk conduction in multi-species electrophoresis. In general, for a multi-species electrophoresis system consisting of *N* species, the Riemann solution consists of at most *N* nonlinear waves corresponding to *N* families of characteristics. These nonlinear waves can be either shock or expansion waves depending upon the initial conditions. We note that the conditions of genuine nonlinearity for all characteristic fields preclude formation of contact discontinuities.

### (a) Shock waves

A shock wave corresponding to *m*th characteristic family, termed *m*-shock, is formed due to intersection of *m*-type characteristics on *x*–*t* plane as shown schematically in figure 2*a*. As the eigenvalues are ordered 0<*ξ*_{1}<*ξ*_{2}<⋯<*ξ*_{N}, the *i*-type characteristics have higher slope in *x*–*t* plane than the *m*-shock if *i*<*m*. Therefore, the characteristics of families 1 to *m*−1 cross the *m*-shock from right as shown in figure 2*b*, whereas the characteristics of families *m*+1 to *N* have a smaller slope than that of the *m*-shock and hence they cross the *m*-shock from left (figure 2*c*). As all the characteristics, except the *m*-type, do not intersect in the vicinity of *m*-shock the values of *R*_{i} for *i*≠*m* can be found by tracing back the values of *R*_{i} to the initial conditions. The value of *R*_{m}, however, changes across *m*-shock as *m*-type characteristics on either side of *m*-shock can be traced back to different initial conditions. Therefore, the change in species concentrations across *m*-shock results only from the change in *m*th Riemann invariant.

We now obtain the conditions on Riemann invariants for existence of shock waves. The generalized solution for a shock wave must satisfy the integral form of conservation laws (2.9) called the Hugoniot conditions,
*f*]=*f*^{+}−*f*^{−} is the difference between downstream (denoted by superscript +) and upstream values (denoted by superscript −) of any physical quantity *f* across the shock and *v* is the shock speed. Therefore, knowing the values of Riemann invariants and concentrations across the shock, the shock speed *v* can be obtained using equation (3.2). Moreover, the generalized solution for *m*-shock must satisfy the following stability conditions
*m*-type characteristics on either side of *m*-shock intersect.

To translate the stability conditions (3.3) for the existence of shock waves in terms of Riemann invariants, we follow an approach which is analogous to that employed by Zhukov [4] for shocks in ITP, albeit in the absence of surface conduction. For our discussion, we would require the following identities
*u*) and Riemann invariants *R* in the polynomial *L*(*u*,*R*) given by equation (2.11). Equations (3.4) and (3.5) imply that *v*(*σ*^{−}_{B}+*k*) and *v*(*σ*^{+}_{B}+*k*) are the roots of polynomials *L*(*u*^{+},*R*) and *L*(*u*^{−},*R*), respectively. Furthermore, using equations (2.17) and (3.6) we observe that all the remaining roots of the polynomials *L*(*u*^{+},*R*) and *L*(*u*^{−},*R*) are identical. As the roots of polynomial *L*(*u*,*R*) are the Riemann invariants, all but one Riemann invariants remain the same across the *m*-shock. Note that if *L*(*u*^{+},*R*) and *L*(*u*^{−},*R*), respectively. As only one Riemann invariant can change across the shock, we take *j*=*i*=*m*. For the *m*-shock to exist, the stability conditions (3.3) must hold. Using *R*_{m}]<0. As all other Riemann invariants are identical across the *m*-shock [*R*_{i}]=0 ∀ *i*≠*m*.

Therefore, the necessary conditions for existence of *m*-shock are [*R*_{m}]<0 and [*R*_{i}]=0 ∀ *i*≠*m*. To show that the converse is also true, we assume that [*R*_{m}]<0 and [*R*_{i}]=0 ∀ *i*≠*m*. Substituting the expression for *u*_{i}(*R*) given by equation (2.18) in the Hugoniot conditions (3.2), we obtain
*σ*_{B}+*k*) from equation (2.15) in equation (3.7), we obtain the dependence of speed *v* of discontinuity on the Riemann invariants,
*R*_{m}]<0. To summarize, *m*-shock exists if and only if [*R*_{m}]<0 and [*R*_{i}]=0 ∀ *i*≠*m*. Moreover, [*R*_{m}]<0 ensures that the stability conditions (3.3) are satisfied.

### (b) Rarefaction waves

If the conditions on Riemann invariants for existence of a shock wave do not hold, it can result in formation of a rarefaction wave. Figure 3 illustrates the formation of *m*-rarefaction corresponding to *m*th characteristic family. The *m*-rarefaction wave is characterized by *m*-type characteristics spreading out in the *x*–*t* plane. The edges of *m*-rarefaction wave are bounded by *m*-type characteristics corresponding to the neighbouring regions with uniform concentrations. Analogous to the case of shock waves, the characteristics of families 1 to *m*−1 cross the *m*-rarefaction from the right, as shown in figure 3*b*, whereas the characteristics of families *m*+1 to *N* have a smaller slope than that of *m*-rarefaction and hence they cross the *m*-rarefaction from left (figure 3*c*). As all characteristics, except the *m*-type, simply pass through the *m*-rarefaction wave without interacting with each other the values of corresponding Riemann invariants can be obtained by tracing back their values to the initial conditions. Therefore, the variation in species concentrations through the *m*-rarefaction wave results only from the change in *m*th Riemann invariant. The *m*th Riemann invariant varies continuously through the *m*-rarefaction wave.

For a Riemann problem, the governing equations (2.16) admit a similarity solution for rarefaction waves with *z*=(*x*−*x*_{0})/(*t*−*t*_{0}) as the similarity variable. Here (*x*_{0},*t*_{0}) are the coordinates of the centre of similarity solution. Substituting *R*_{m}=*R*_{m}(*z*) in the governing equations (2.16), we obtain
*N* non-trivial solutions for Riemann invariants of the type *ξ*_{m}(*R*)=*z* and *R*_{i}= constant for *i*≠*m*. Therefore, the variation of Riemann invariants through *m*-rarefaction wave is given by
*ψ*_{m} can be obtained using (2.16) and noting that *ξ*_{m}=*z* inside the *m*-rarerfaction wave,

The *m*-rarefaction wave of the form given by equation (3.10) connects two constant solutions of the equation (2.16) *R*^{−} and *R*^{+} cannot be connected by *m*-shock when *m*th characteristic field ensures that *ξ*_{m} varies monotonically with *z* within the *m*-rarefaction wave [29].

### (c) Construction of solution

Having obtained the conditions on Riemann invariants for existence of shock and rarefaction waves we can construct the generalized solution for the Riemann problem for system (2.16), that is, with piecewise constant initial values of Riemann invariants with a single discontinuity,
*N* nonlinear waves. Based on the earlier discussion we know that only one Riemann invariant *R*_{m} corresponding to the *m*th characteristic field can change across the *m*-shock or *m*-rarefaction wave. The remaining Riemann invariants are the same across the shock or rarefaction wave and can be found by tracing back their values along the characteristics to the initial data. Therefore, to construct the solution we assume *N* nonlinear waves emanating from the initial point of discontinuity. On both sides of the *m*th nonlinear wave (either shock or rarefaction wave) *i*<*m* and *i*>*m*, whereas *m*th nonlinear wave. After obtaining the values of Riemann invariants across all the nonlinear waves, each wave can be classified as shock or rarefaction based on the conditions derived above. That is, the nonlinear wave is *m*-shock if *m*-rarefaction wave if *m*th wave is a shock wave, the shock speed can be obtained using equation (3.8), whereas for rarefaction wave the solution can be completed using equation (3.10). Of course, if *k* then there will be no nonlinear wave corresponding to the *k*-th characteristic field.

When the initial conditions have only one discontinuity, the nonlinear waves do not interact with each other. For practical electrophoresis systems the initial conditions often have piecewise-constant distribution of concentrations (or Riemann invariants) with several discontinuities. In such cases, nonlinear waves emanating from different points of initial discontinuity can interact. The interaction of nonlinear waves can be analysed by treating the instant of wave interaction as a new initial condition. In the current work, we consider electrophoretic systems involving interaction of only shock waves. Therefore, the instant of shock interaction can be regarded as a new piecewise-constant initial condition of the form (3.12) for which the solution procedure has been provided above.

## 4. Results and discussion

In this section, we illustrate the effect of surface conduction on the propagation of nonlinear concentration waves through examples of various electrophoretic systems. We present analytical solutions to the set of conservation laws (2.9) for different initial conditions using the method of characteristics presented in §3. We begin by reviewing the formation of shock and rarefaction waves in binary electrolyte systems in the presence of surface conduction. Thereafter, we present analytical solutions for multi-species electrophoretic systems such as ITP and compare our results with similar electrophoretic systems but with negligible surface conduction.

### (a) Electromigration in binary electrolytes

We begin by illustrating the effect of surface conduction on the evolution of conductivity gradients in a binary electrolyte due to applied electric field. This is of practical significance as gradients in binary electrolyte are observed during FASS and concentration polarization at the interface of a microchannel and nanochannel [13]. The latter phenomenon has found applications in shock-electrodialysis for deionization of water [10]. Nonlinear waves in binary electrolytes due to the competing effects of bulk and surface conduction have been analysed theoretically by Mani *et al.* [9] and observed experimentally by Zangle *et al.* [11]. Here we review the problem of nonlinear waves in binary electrolytes considered by Mani *et al.* [9] to demonstrate the application of our generalized mathematical model and solution methodology for multi-species electrophoretic systems. We will use this analysis to explain the results for multi-species electrophoretic systems in later sections.

We consider a binary electrolyte system consisting of a co-ionic (anion) and a counter-ionic (cation) species denoted by subscripts 1 and 2, respectively. Without loss of generality, the surface charge is assumed to have negative polarity. The governing equations (2.9) in terms of Riemann invariants for a binary electrolyte system simplify to
*R*_{1}. Therefore, the electrophoretic system can exhibit shock and rarefaction waves depending upon the initial conditions. This is in contrast with the case when surface conduction is negligible (in large microchannels) where concentration gradients in binary electrolytes remain stationary due to the zero eigenvalue of the electrophoretic system [16].

To demonstrate the effect of surface conduction on concentration wave propagation in binary electrolytes, we consider an idealized case of a binary electrolyte system with an initial co-ion concentration distribution *c*_{1}(*x*,0) shown in figure 4*a*. For our calculations, we assume that the channel has a uniform surface charge density with *ρ*_{w}*P*/*A*=−1.0×10^{6} C m^{−3}. This value can correspond to, for example, surface charge density of 0.25 C m^{−2} and a circular channel cross section with diameter of 1 μm. The binary electrolyte is sodium chloride and hence chloride (*z*=−1, *μ*=−79.1×10^{−9} m^{2} V^{−1} s^{−1}) is the co-ion and sodium (*z*=1, *μ*=51.9×10^{−9} m^{2} V^{−1} s^{−1}) is the counter-ion. The total current density, *I*_{T}/*A*, through the solution is taken as 2450 *A* m^{−2} with electric field vector pointing towards the left. This value is typical of current densities used in electrophoresis experiments [32]. The initial conditions shown in figure 4*a* can be written in terms of co-ion concentration *c*_{1} and Riemann invariant *R*_{1} (using equation (2.11)) as
*c*_{H} and *c*_{L} denote high and low co-ion concentrations, respectively.

Following the solution procedure outlined in §3, this initial condition results in formation of two nonlinear waves emanating from the initial points of discontinuity, (*x*_{1},0) and (*x*_{2},0), as shown in figure 4*b*. These waves separate the *x*–*t* plane into three zones I, II and III depicted in figure 4*b*. Across the trailing wave originating from (*x*_{1},0), [*R*_{1}]=[*α*_{1}*k*/(*β*_{1}*c*_{1}+*k*)]<0 as *c*_{H}>*c*_{L}. Therefore, the trailing wave is a 1-shock wave. Knowing the value of Riemann invariant *R*_{1} across the 1-shock in zones I and II, the shock speed *v*_{1S} can be calculated using equation (3.8) to get
*x*_{2},0), *ξ*_{1}. As the values of *R*_{1} in zones II and III remain unchanged from the initial conditions, propagation speeds of the edges of 1-rarefaction wave are given by
*R*_{1} inside the rarefaction wave is given by equation (3.10). Therefore, the distributions of Riemann invariant *R*_{1} and co-ion concentration for *t*>0 are given by

The concentration distribution of co-ion at a later time (*t*=20 s) is shown in figure 4*c*. The left interface, which is a shock wave, remains sharp and propagates towards the right, whereas the right interface, which is a rarefaction wave, disperses over time while propagating towards the right. These results can be explained qualitatively by noting that the wave speed *R*_{1} varies inversely with *c*_{1}. That is, the point with higher concentration has lower wave speed and vice versa. Therefore, the positive gradient in the initial co-ion concentration remains sharp and propagates as a shock wave whereas the negative gradient in co-ion concentration disperses to form a rarefaction wave as the leading edge of a wave with lower concentration has a higher wave speed than the trailing edge with higher concentration and lower wave speed.

### (b) Shock and rarefaction waves in discontinuous electrolyte system

Next, we analyse the propagation of concentration waves due to the competing effects of bulk and surface conduction in an electrolyte system consisting of two co-ionic species and one counter-ionic species. We assume that the initial distribution of co-ions has a discontinuity as shown schematically in figure 5*a*. Such an electrolyte system is used in ITP, where co-ion (1) is the TE and co-ion (2) is the LE. Although in ITP, LE ion has higher mobility than the TE ions, here we consider a generalized case where co-ion 1 may have higher or lower mobility than co-ion 2.

Upon application of an electric field, the discontinuous initial conditions shown in figure 5*a* give way to two nonlinear waves originating from the point of discontinuity (*x*_{0},0) (figure 5*b*). These waves divide the *x*–*t* plane into three zones (I–III). Following our discussion on the construction of a solution in §3c, the two nonlinear waves can be either shock or rarefaction waves depending on the values of Riemann invariants (*R*_{1} and *R*_{2}) in zones I–III. To obtain the concentrations and the Riemann invariants in various zones, we begin by writing the initial conditions in terms of Riemann invariants. Knowing the initial species concentrations and using equation (2.11), we obtain the initial values of Riemann invariants (table 1). Here, we have denoted the initial values of Riemann invariants in front and behind the initial discontinuity with subscripts + and −, respectively. We note that during the assignment of initial values of Riemann invariants to first and second characteristic fields we have assumed that *R*_{1} and *R*_{2} in zones I–III by tracing back their values to the initial conditions as shown in figure 5*b*. Thereafter, we obtain the species concentrations in all the zones using equation (2.18) (table 1).

Having obtained the values of Riemann invariants in all the zones, we can check for the conditions on Riemann invariants provided in §3a,b to classify the nonlinear waves as shock and rarefaction waves. That is, if *α*_{1}<*α*_{2} (or |*μ*_{1}|<|*μ*_{2}|). The latter condition is the same as that for existence of shock wave in conventional ITP. That is, the LE ion should have higher mobility than the TE ion for a shock to occur between the LE and TE zones. If the conditions for shock waves do not hold, the corresponding wave evolves as a rarefaction wave. Finally, to obtain the complete co-ion distribution we use equation (3.8) to calculate the shock speeds and similarity solution (3.10) to find the structure of the rarefaction wave.

To illustrate the application of the above analytical solution, we consider a practical anionic ITP system in a channel with negatively charged walls. Here, the TE ion (1) is methyl sulphonate (*z*=−1, *μ*=−50.6×10^{−9} m^{2} V^{−1} s^{−1}), the LE ion (2) is chloride (*z*=−1, *μ*=−79.1×10^{−9} m^{2} V^{−1} s^{−1}), and the background ion (3) is sodium (*z*=+1, *μ*=51.9×10^{−9} m^{2} V^{−1} s^{−1}). We assume that the channel walls have uniform charge density with *ρ*_{w}*P*/*A*=−1.0×10^{6} *C* m^{−3} and the total current density is 2450 *A* m^{−2}. We consider two cases, shown in figures 6*a* and 7*a*, with the same LE ion concentration but different TE ion concentrations. For both cases, application of an electric field results in the formation of a 2-shock wave as the mobility of the LE ion is greater than that of a TE ion. However, as discussed above, the first nonlinear wave can be a shock or rarefaction wave depending upon the initial concentration of the TE ion relative to that of the LE ion.

When the initial TE ion concentration is low (figure 6*a*), the first nonlinear wave is a 1-shock wave. This is characterized by a decrease in value of *R*_{1} and a corresponding increase in TE concentration from zone I to II as shown in figure 6*b*,*d*. The formation of 1-shock can be explained as follows. Upon application of an electric field, the TE ions displace the LE ions and simultaneously the concentration of TE ions adjusts to a new value so as to maintain continuity of current. This results in a positive gradient in TE concentration. As discussed in §4*a*, a positive concentration gradient in TE, which is locally a binary electrolyte, sharpens to form a shock wave.

In figure 6, we also present results for the case when surface conduction is negligible using the analytical solution provided by Zhukov [4]. When surface conduction is absent (*ρ*_{w}=0), ion transport in ITP results in formation of a propagating 2-shock wave separating LE and TE zones. Behind the 2-shock, the TE concentration adjusts to a new value so as to maintain the continuity of the current. The adjusted TE concentration is governed by the Kohlrausch regulating function [16] and is equal to *ρ*_{w}=0. Another difference between ITP with and without surface conduction is that the speed of 2-shock separating LE and TE zones is lower in the presence of surface conduction compared with the case when surface conduction is absent. This difference in shock speeds can be attributed to the fact that only a part of the total current flows through the bulk solution when surface conduction is present. Consequently, the magnitude of the electric field in the bulk solution is lower in the presence of a surface condition which reduces the speed of 2-shock.

In figure 7, we present results for the case when the initial TE concentration is higher compared with the adjusted TE concentration. The increase in TE concentration has no effect on the speed of 2-shock and the concentration of TE in the adjusted TE zone for cases with and without surface conduction. However, the negative concentration gradient that results from adjustment of TE causes the interface to disperse in the presence of surface conduction. Consequently, the interface separating the initial and the adjusted TE zones results in the formation of a rarefaction wave as shown in figure 7*b*,*d* whereas in the absence of surface conduction, the concentration gradient in TE is a stationary contact discontinuity as shown in figure 7*c*.

### (c) Preconcentration and separation in ITP

Lastly, we consider a practical case of simultaneous preconcentration and separation of analytes in plateau-mode ITP in the presence of surface conduction. In plateau-mode ITP, a mixture of analytes is initially injected between zones of LE and TE. Upon application of an electric field, the analytes segregate into purified zones and simultaneously their concentrations adapt to the concentration of LE, resulting in preconcentration [2]. The dynamics of this electrophoretic process involves formation and interaction of several nonlinear waves (shock and rarefaction waves). After the interaction of nonlinear waves subsides, the process is characterized by purified analytes zones between LE and TE zones. These zones are separated by sharp zone boundaries which are shock waves in ion concentrations.

We here consider an anionic ITP system wherein LE ion is chloride (*z*=+1, *μ*=−79.1×10^{−9} m^{2} V^{−1} s^{−1}), TE ion is methyl sulfonate (*z*=−1, *μ*=−50.6×10^{−9} m^{2} V^{−1} s^{−1}) and background ion is sodium (*z*=+1, *μ*=51.9×10^{−9} m^{2}V^{−1} s^{−1}). We consider two analytes, fluoride (*z*=−1, *μ*=−57.4×10^{−9} m^{2} V^{−1} s^{−1}) and nitrate (*z*=−1, *μ*=−74.1×10^{−9} m^{2} V^{−1} s^{−1}). The current density is taken to be 2450 Am^{−2}. We assume that the channel wall has a uniform charge density with *ρ*_{w}*P*/*A*=−1.0×10^{6} *C* m^{−3}. In our analysis, we denote TE, fluoride, nitrate, LE and background ions with subscripts 1 to 5, respectively. Figure 8*a* shows the initial conditions where the mixture of analytes is initially sandwiched between LE and TE zones. The generalized solution for this type of initial condition in the absence of surface conduction has been provided by Babskii*et al.* [3]. Figure 8*b*,*c* show the spatial variation of species concentrations at *t*=25 s and *t*=50 s, respectively, when surface conduction is absent. The analytes segregate into purified zones over time. The local TE concentration adapts to the Kohlrausch regulating function [16] set by the initial condition and the concentration gradients in the TE zone remain stationary in the absence of surface conduction.

Next, we use the solution methodology discussed in §3*c* to construct the analytical solution when surface conduction is present. Figure 8*d* shows the spatiotemporal evolution of various zones and concentration waves on the *x*–*t* diagram. The *x*–*t* diagram shows interaction of shock waves which is accompanied by extinction of old zones and emergence of new zones. In figure 8*e*,*f*, we show the spatial variation of species concentrations at *t*=25 s and *t*=50 s, respectively, when surface conduction is taken into account. At *t*=25 s, the analytes are partially separated and there is a zone around *x*=50 mm where the analytes are still mixed. Moreover, there is an additional rarefaction wave at *x*=25 mm and a shock wave at *x*=42 mm in the TE zone. These nonlinear waves in TE zone are not present in conventional ITP where surface conduction is absent. Instead, in the absence of surface conduction the TE zone has contact discontinuities corresponding to the initial discontinuities in the Kohlrausch regulating function [16] as shown in figure 8*b*,*c*. At *t*=50 s, the analytes are completely segregated into purified zones, separated from other zones with propagating shock waves. Note that, the rarefaction and the shock wave in the TE zone shown in figure 8*e*,*f* also propagate towards the right. Such propagating concentration gradients will be detected during ITP experiments along with the gradients associated with analyte zones as they pass over a point-detector, such as conductivity and absorbance detectors [33]. Therefore, additional nonlinear waves resulting from the competition of bulk and surface conduction must be taken into account while interpreting electrophoresis experiments, wherein surface conduction is appreciable.

## 5. Conclusion

We have developed a theory for multi-species electrophoretic transport of ions in shallow microchannels and nanochannels, where surface conduction through EDL competes with bulk conduction. Our model is applicable for arbitrary large number of co-ions and a single counter-ionic species having opposite charge polarity compared with the surface charge. Based on our mathematical model, we show that the Kohlrausch regulating function [16], which is widely used to calculate zone concentrations in conventional electrophoresis, is not valid in the presence of surface conduction. This is because electrophoretic systems with appreciable surface conduction do not have zero valued eigenmobility (eigenvalue of the Jacobian of flux). For the same reason, such systems do not exhibit contact discontinuities in species concentrations. Instead, nonlinearities resulting from surface conduction and multi-species electromigration cause concentration gradients to propagate only as shock or expansion waves. To analyse the propagation of such nonlinear waves and construct analytical solutions to the governing equations, we have presented a generalized solution methodology based on the method of characteristics. In particular, we have derived the necessary and sufficient conditions for the existence of shock and expansion waves in electrophoretic processes where bulk and surface conduction compete with each other.

We have demonstrated the application of our model and solution methodology for analysing practical electrophoretic systems employing binary and multi-species electrolytes. In particular, we used the example of ITP to highlight the effects of surface conduction on electrophoretic processes. By comparing the concentration wave propagation in ITP with and without surface conduction, we showed that the propagation speed of ITP shocks decreases in the presence of surface conduction. This is because, for a fixed total current, surface conduction leads to reduction in current through the bulk solution which in turn reduces the shock speeds. Moreover, the competing effects of bulk and surface conduction result in lower levels of sample preconcentration in ITP compared with conventional ITP. Owing to the absence of zero eigenmobility of the electrophoretic system in the presence of surface conduction, the interface formed between the initial TE zone and the adjusted TE zone in ITP does not remain stationary. Instead, this interface propagates as a shock or an expansion fan, depending upon the initial TE concentration.

Although, here we have demonstrated the application of our model for analysis of electrophoresis techniques, our theory is equally applicable for understanding other electrokinetic processes where surface conduction is important, including shock-electrodialysis based water deionization [10] and concentration polarization in DC electro-osmotic pumps [34]. Current analyses of such processes are limited to those with binary electrolytes even though in practical applications one would expect the presence of more than two ionic species. Our theory which takes into account the combined effect of nonlinearities induced by surface conduction and multi-species electromigration can help in better understanding and design of these electrokinetic processes.

## Data accessibility

This article does not report any primary data. All computations were performed in Matlab (Mathworks).

## Authors' contributions

S.S.B. developed the mathematical model and the solution methodology. R.M. and M.K. obtained the analytical solutions and prepared the figures. S.S.B. drafted the manuscript. All authors reviewed the manuscript and gave final approval for publication. All authors contributed equally to this work.

## Competing interests

The authors declare no competing interests.

## Funding

The work was financially supported by the Industrial Research and Development (IRD) unit of Indian Institute of Technology Delhi under IITD/IRD/MI01068.

- Received September 16, 2015.
- Accepted January 25, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.