Dynamic Stern layers in charge-regulating electrokinetic systems: three regimes from an analytical approach

We present analytical solutions for the electrokinetics at a charged surface with both non-zero Stern-layer conductance and finite chemical reaction rates. We have recently studied the same system numerically [Werkhoven {\em et al.}, Phys. Rev. Lett. {\bf 120}, 264502 (2018)], and have shown that an applied pressure drop across the surface leads to a non-trivial, laterally heterogeneous surface charge distribution at steady state. In this work, we linearise the governing electrokinetic equations to find closed expressions for the surface charge profile and the generated streaming electric field. The main results of our calculations are the identification of three important length and time scales that govern the charge distribution, and consequently the classification of electrokinetic systems into three distinct regimes. The three governing time scales can be associated to (i) the chemical reaction, (ii) diffusion in the Stern layer, and (iii) conduction in the Stern layer, where the dominating (smallest) time scale characterises the regime. In the reaction-dominated regime we find a constant surface charge with an edge effect, and recover the Helmholtz-Smoluchowski equation. In the other two regimes, we find that the surface charge heterogeneity extends over the entire surface, either linearly (diffusion-dominated regime) or nonlinearly (conduction-dominated regime).


I. INTRODUCTION
While the eld of electrokinetics is over a century old, interest in it has only grown since its foundations were laid by Helmholtz and Smoluchowski [1,2]. In addition to applications in well-established elds such as geology [3] and catalysis [4], the advent of micro-and nano-uidics renewed interest in electrokinetics [5][6][7][8][9] due to applications in, for example, blueenergy harvesting [10,11]. At the basis of all eletrokinetic systems is the interaction between uid ow and a charge current. In a closed-circuit setup, an imposed voltage drop over a channel with a charged wall induces uid ow via the electric body force in the Navier-Stokes equations, while in an open-circuit setup, an imposed pressure drop induces an electric eld. is streaming electric eld implies a voltage drop across the channel, the so-called streaming potential ∆Φ. In the linear response regime, the generated streaming potential is linearly related to the applied pressure drop ∆p via the Helmholtz-Smoluchowksi equation [1,2], where G is the channel conductivity, and η are the permi ivity and shear viscosity of the liquid, respectively, and ζ is the zeta potential -the electrostatic potential at the slipping plane of the charged surface of the channel. Since G, and η are material properties of the liquid, measurements of ∆Φ at a known ∆p allows one to use Eq. (1) to measure ζ, an important surface property. e presence of the charged walls is essential for electrokinetic phenomena, since charged surfaces induce an Electric Double Layer (EDL) in the uid adjacent to the surface. is di use layer of ions screens the charge of the surface. A uid ow (electric eld) through the charged EDL induces a charge current (body force), which in turn induces an electric eld ( uid ow). A detailed understanding of the surface charge is therefore vital to describe electrokinetic systems accurately. e total channel conductivity is typically decomposed as G = G b + G s /H, with G b the bulk conductivity of the uid, G s the surface conductivity, and H the channel height [12]. e signi cance of the surface conductivity is expressed by the Duhkin number Du=G s /(G b H) [12]. e surface conductivity G s can be further decomposed as G s = G EDL s + G S s . Here, G EDL s originates from the increased conductivity of the EDL with respect to bulk uid, as rst described by Bikerman [13,14], while G S s is the conductivity due to the mobile charges in the Stern layer [15], the quasi-2D layer in which the surface charge resides in. It is well known that for a wide variety of materials, including insulating materials such as glass or clay [16][17][18][19][20][21], the charges in the Stern layer are not static. It has even been previously shown that the mobility of the Stern-layer charges is comparable to the mobility of ions in bulk [16,22,23].
In the most common case, the surface charge in the Stern layer originates from a chemical reaction with dissolved ions in the uid, either via an adsorption or a desorption reaction. e surface charge is therefore not a xed quantity, but is determined by a charge regulation process [24]. In this work we will consider the desorption reaction SC S − +C + , where a neutral surface group SC releases a cationic counter ion C + leaving behind a charged, covalently bound surface group S − . is reaction represents, for example, a deprotonation reaction if we identify C + as a proton. In equilibrium, the balance of this reaction is given by the Langmuir desorption isotherm f = (1 + ρ C,s /K) −1 [24][25][26], with f the fraction of charged surface groups, K the chemical equilibrium constant, and ρ C,s the counter ion density at the surface. However, the Langmuir isotherm assumes chemical equilibrium at all times, as in previous theories of charge regulation in electrokinetic systems [27,28], and therefore does not take nite chemical reaction rates into account. We have recently shown, however, that the combination of nite chemical rates with a non-zero Stern-layer mobility leads to novel electrokinetic properties [29]. In particular, our numerical solutions showed that a lateral uid ow induces a heterogeneous surface charge on a chemically homogeneous, nite surface, which furthermore provided a rst-principles explanation for the experimentally observed in uence of uid ow on the surface chemistry [30]. In section II of this paper, we linearize the governing equations used in Ref. [29], and solve them analytically with the aid of several simplifying approximations. e analytical solutions exhibit qualitative agreement with the numerical results, but due to the underlying approximations we do not obtain full quantitative agreement. Nevertheless, our analytical approach allows us to identify three important length and associated time scales, summarised in section III A, as well as three distinct and qualitatively di erent regimes. For fast reaction rates, the reaction-dominated regime discussed in section III B, we reproduce a constant surface charge and the Helmholtz-Smoluchowski equation (1), while for nite reaction rates we obtain either a linear surface charge pro le (di usion-dominated regime, section III C), or a nonlinear pro le (conductiondominated regime, section III D). Our analysis allows one to categorise all electrokinetic setups in three regimes, such that the characteristics of the surface charge distribution can be predicted or tuned.

II. THE (LINEARISED) POISSON-NERNST-PLANCK EQUATIONS
We consider an electrokinetic system consisting of two water reservoirs connected by a rectangular channel with height H and length 2L. We denote the normal and lateral Cartesian coordinate by z ∈ [0, H] and x ∈ [−L, L], respectively, and assume translational invariance in the y direction. e top (z = H) and bo om (z = 0) surface of the channel carry chargeable surface groups, which for simplicity we assume to be equal such that the plane at z = 1 2 H is a symmetry plane of the system. Without loss of generality, we take the uid ow in the positive x direction. e reservoirs contain three di erent ionic species labelled by i = +, −, C, with valency z + = −z − = 1 and z C . Charge neutrality in the bulk demands i z i ρ b,i = 0, with ρ b,i the bulk concentration of ion i in the two reservoirs, i.e. we do not consider di usio-osmotic processes here. is xes the Debye screening length, as λ D = k B T /(e 2 i z 2 i ρ b,i ), the typical thickness of the EDL, with k B the Boltzmann constant, T the temperature, and e the elementary charge.
We denote the position and time dependent concentration and ux of the three ion species by ρ i (r, t) and J i (r, t), respectively, the electric potential by ψ(r, t) and the uid velocity and pressure by u(r, t) and p(r, t), respectively. ese quantities are governed by the Poisson-Nernst-Planck-Navier-Stokes (PNP-NS) equations [31], where D is the di usion constant (assumed for simplicity to be equal for all ionic species), and m the mass density of the uid. e PNP-NS equations combine the Poisson equation for the electrostatic potential, the incompressible Navier-Stokes equation for the uid ow, the Nernst-Planck equation for ionic transport and the continuity equation for the ion densities. ese equations are then to be coupled to a dynamic Stern layer. In our theoretical framework, we treat the density of surface charges σ(x, t), a 2D analogue of ρ i (r, t), as a dynamic variable. e surface charges are produced by a chemical reaction SC S − +C + , and therefore σ(x, t) is not necessarily a locally conserved quantity. However, since the total number of counter ions must be conserved, the production rate R of surface charges is equal to the counter-ion ux leaving the surface n s · J C,s , with n s the normal vector of the solid-liquid interface pointing into the liquid. Denoting the ux of surface charges in the x-direction by j σ (x, t), the 2D analogue of J i (r, t), we can write the continuity equation for σ(x, t) as e production rate of surface charges is governed by the chemical rate equations, where the surface charge production (annihilation) rate is proportional to the density of uncharged (charged) sites. Furthermore, we assume a Nernst-Planck-like equation for j σ , derived via a dynamical density functional theory in Appendix A, with the noticeable absence of convection in the Stern layer. e resulting equations that govern the surface dynamics, and therefore can be seen as the boundary conditions to the PNP-NS equations (2), are given by with D s is the di usion constant of the surface charges, z σ = −z C the valency of the surface charges, Γ the density of surface sites, and k ads and k des the adsorption and desorption rate constants, respectively. Here, and throughout this work, we use the subscript "s" to denote surface quantities, such that u s = 0 enforces a no-slip boundary condition and ρ C,s (x) ≡ ρ C (x, z = 0) the counter ion concentration at the surface. e di usion contribution (∝ ∂ x σ) to j σ must in general be adjusted by a factor (1 − σ/Γ) −1 as the surface groups cannot be multiply occupied [32]. Furthermore, we assume that D s does not depend the surface concentration σ. For high concentrations, however, the di usion constant depends non-trivially on the concentrtaion [33]. However, since at most a few percent of the total number of sites is charged, σ is su ciently small to safely assume D s to be constant. Similarly, we have assumed that the bulk di usion constant D in Eq. (2) is homogeneous throughout the system, while in general this depends on for example the distance to the surface [34]. For simplicity, we leave out these higher order e ects. In equilibrium all uxes vanish, and for J C = 0 we recover from Eq. (6) the standard Langmuir desorption equation σ = Γ/(1 + ρ C,s /K), with K = k des /k ads equal to the chemical equilibrium constant of the reaction SC S − +C + . e governing equations (2)-(6) cannot be solved analytically in general. In this article, however, we will show how to obtain approximate solutions to these equations. As is common in pressure-driven electrokinetic systems, we neglect not only the inertial terms in the Navier-Stokes equation (low Reynolds number, e ectively m = 0), but also the electric body forces on the uid [31], such that we can ignore the nal term on the right hand side of the Navier-Stokes equation (2). e la er approximation can be justi ed for setups with H λ D , and by realising that the body force is localised in the EDL while the pressure gradient extends over the entire channel height. is approximation is con rmed by our numerical calculations (see Appendix B), which show that including the electric body force has no signi cant e ect on the steady state surface charge pro le.
e Navier-Stokes equation then reduces to the Stokes equation, and is now decoupled from other quantities. is allows us to solve for u and p, resulting for an applied pressure drop ∆p to the standard Poisseuille ow, withx the unit vector in the x direction and ∂p ∂x = − ∆p 2L . e typical ow velocity can be estimated from Eq. (7). For pressure drops ∆p no larger than 1 bar, η ∼ 1 mPa s, and channel dimensions H ∼ 1 µm and L ∼ 10 µm, one arrives at u x (λ D ) < 10 −2 m/s for physically relevant salt concentrations in water.
In the following analysis we will focus solely on the net charge density ρ e (r, t) = i z i ρ i (r, t) and electric current J e (r, t) = i z i J i (r, t). e la er can be wri en as with ρ(r, t) = i z 2 i ρ i (r, t) the total local ionic strength. As is common in a linearised theory of electrokinetic systems, we assume here that ρ is constant throughout the system and equal to its bulk value, ρ = i z 2 i ρ i,b . At steady state, ∇ · J e = 0, together with the use of the Poisson equation to eliminate the electric potential in favour of the charge density, we obtain the governing equation for ρ e (x, z), We can simplify Eq. (9) using scaling arguments. From our previous work, we know that the uid ow will induce heterogeneities in the x direction, and therefore we estimate that ∂/∂x ∼ 1/L. We also know that ρ e reduces quickly to 0 within a few λ D in the z direction, and hence ∂/∂z ∼ 1/λ D . is allows us to de ne a few characteristic time scales for the EDL, .
Here, τ L is the characteristic time for an ion in bulk to di use the lateral length L, and is of the same order of magnitude as the rst term of Eq. (9). Additionally, τ EDL is the characteristic equilibration time of an EDL, and is of the same order of magnitude as the second and third term of Eq. (9). Lastly, τ adv is the characteristic time it takes for an ion in the EDL to be advectivelly transported from one end of the channel to the other. Since ρ e is only non-zero in the EDL, we evaluated u x at z = λ D . For our geometry, L λ D , and thus we can conclude that τ EDL τ L , meaning that the di usion in the lateral direction ( rst term Eq. (9)) is negligible compared to di usion in the normal direction (second term Eq. (9)). Moreover, since D/λ D ∼ 0.1m/s> u x (λ D ) for our parameter choice of interest, we have τ EDL τ adv implying that convection in the lateral direction (last term Eq. (9)) is negligible with respect to di usion in the normal direction. A er plugging in typical quantities, one indeed nds that τ EDL is of the order of (tens of) nanoseconds, while τ adv is of the order of miliseconds or larger. Eq. (9) now reduces to a simple di erential equation, where ζ(x) is an integration constant that remains to be found. As we will show below, we can identify ζ(x) = βe(ψ(x, 0) − ψ(x, 1 2 H)) as the (heterogeneous) dimensionless zeta potential at steady state. e second integration constant has been set to zero since ρ e (x, z → 1 2 H) = 0 from H λ D . Eq. (11) is analogous to the equilibrium linear Poisson-Boltzmann equation for the charge density, and is a direct consequence of τ EDL τ adv : the EDL is equilibrated in the z direction as convection is typically not strong enough to deform the EDL signi cantly.
To determine the integration constant ζ(x), we apply the boundary conditions Eqs. (3) and (4)- (6). To facilitate the calculations, we make use of the fact that for the majority of surfaces only a small fraction of the sites are charged, σ Γ. For instance, only a few percent of the surface groups of silica are charged under typical conditions (3 <pH< 11, 1 mM< ρ s <100 mM) [35]. With this assumption, the equation for σ(x) simpli es to Eq. (12) constitutes a di usion-conduction-reaction problem coupled to the 3D-channel via the in-plane electric eld ∂ψs ∂x and the counter ion density ρ C,s . As a consequence, three regimes will arise depending on which process is dominant.
is is reminiscent of a convection-di usion problem [36] with a linear/exponential density pro le for di usion/convection dominated systems. In our case, the role of convection in the Stern layer is played by conduction, but we will analogously nd a linear/exponential in the di usion/conduction limited regime.
Since the uid ow is in the positive x direction, the streaming electric eld −∂ψ s /∂x must have the same sign as the surface charge such that no net charge is transported between the two reservoirs. It is convenient to separate the sign and magnitude of the streaming electric eld. us, we de ne −βe∂ψ s /∂x = z σ E, where E is a positive quantity with dimensions of inverse length. To solve for σ(x), we further assume that both ρ C,s and E are spatially constant. While this is a valid assumption in simple electrokinetic systems, we have shown recently [29] that such approximations are no longer valid when both Sternlayer conduction and nite chemical rates are taken into account; we found that a heterogeneous surface charge leads to a heterogeneous streaming electric eld and counter ion density along the surface. Nevertheless, approximating E and ρ C,s to be spatially constant allows us to solve for σ(x) and determine ρ e , J e and E. In principle one could then reinsert the solutions in Eq. (12) and obtain an improved σ(x), ρ C,s and E. However, in this work we will refrain from applying an iterative scheme and aim for an analytical and qualitative understanding of the electrokinetic phenomena.
For a spatially constant E and ρ C,s , Eq. (12) is straightforward to solve. Since z 2 σ = 1, we nd with a ± integration constants and σ eq ≡ Γ(1+ ρC,s K ) −1 the equilibrium surface charge density. e wavenumbers k ± = 1 2 E ± 1 2 E 2 + 4λ −2 reac set the relevant lateral length scales, with λ reac ≡ √ D s τ reac (discussed in more detail in section III A), and τ reac = k des + k ads ρ C,s −1 the characteristic time scale of the chemical reaction. e amplitudes a ± can be determined by imposing that the surface current vanishes at the end-points of the charged surface, j σ (±L) = (−D s ∂ x σ +D s Eσ)| x=±L = 0. e integration constants a ± are then found to be Note that σ(x) does not depend on z σ , and that in equilibrium, E → 0, and hence σ(x) → σ eq . For non-zero E, we nd a double exponential pro le. is reduces to either a linear pro le for k ± L 1, or a single exponential pro le if k ± L 1. e nal unknown, the integration constant ζ(x), can be determined using Eq. (6), coupling σ(x) to ρ e (x, 0), see Appendix C. To facilitate this calculation, we rewrite the counter ion ux asẑ · J C (x, 0) = z Cẑ · J e (x, 0), where we have used n s =ẑ, z · J + (x, 0) =ẑ · J − (x, 0) = 0. Using the solutions to ρ e , Eq. (11), and σ, Eq. (13), we nd Here, we identi ed the dimensionless equilibrium zeta potential ζ eq = z σ 4πλ B λ D σ eq from linear Poisson-Boltzmann theory. Comparing Eq (13) with Eq. (15), we see that ζ(x) is proportional to the steady state surface charge, and indeed can be interpreted ζ(x) as the steady state zeta potential.
To determine E, we impose that at any position x no net current passes through any channel slice with normalx. is condition is a direct consequence of the vanishing divergence of J e and the open circuit geometry, hence where z σ j σ (x) is the net charge current through the Stern layer, determined using Eqs. (13) and (5), and due to symmetry we only integrate over half the channel height. Eq. (16) is a local condition, and we will in general nd E to depend on x (see below). is, however, contradicts our previous assumption that E is spatially constant, and it is at this point that our analytic approach is inconsistent. Nevertheless, the results of this analytic approach allow us obtain approximate solution, and agree qualitatively well with the full numerical solutions. Despite the inconsistency, this approach gives us physical insight in the system, and allows us to identify three separate regimes.

A. Governing time and length scales
We identify two physically important length scales, that appear in the de nition of the wavenumbers k ± , We can interpret the rst length scale λ reac as the typical distance a surface charge traverses di usively during a time τ reac , that is, the typical distance travelled between ad-and desorption. e conductive length scale λ cond can be interpreted as the distance a monovalent ion needs to travel in order to gain an energy equal to k B T due to the streaming electric eld. e dynamics of the system is fully determined by λ reac , λ cond , and the channel length L.
Alternatively, we can identify an equivalent time scale for each length scale using the surface di usion constant D s . Eq. (17) shows that τ reac is the equivalent time scale of λ reac . By introducing a conductive velocity v cond = D s E and a di usive velocity v dif = D s /L, we can transform L and λ cond into equivalent time scales, We can interpret τ dif as the characteristic time for a Stern-layer charge to di use across the channel length, and τ cond as the characteristic time a er which a Stern-layer charge has gained one thermal energy unit due to the streaming electric eld. Together with τ reac , the three time scales can be used equivalently to the three length scales to characterise the electrokinetic system, as we can express the ratio between every pair of length scales as the ratio between the two equivalent time scales: e three distinct characteristic times allow us to identify three electrokinetic regimes, de ned by the smallest time (or associated length). While the three length scales appear naturally in the analytical description, we found it more intuitive to consider the three time scales when considering the di erent dynamical regimes. In the reaction-dominated regime, characterised by the near-equilibrium of the adsorption/desorption process, τ reac τ cond , τ dif (and hence λ reac L, λ cond ), to be discussed in section III B, we obtain the standard Helmholtz-Smoluchowski picture with a constant surface charge and electric eld, except for a region of size λ reac around the edges at x = ±L. In the di usion-dominated regime, τ dif τ cond , τ reac (L λ cond , λ reac ), discussed in section III C, the surface charge is heterogeneous over the entire surface and linear in the lateral position x. Consequently, the streaming electric eld E is also heterogeneous, but we nd that the streaming potential approximately conforms to the Helmholtz-Smoluchowski equation Eq. (1). However, in the conduction-dominated regime τ cond τ reac , τ dif (and hence λ cond λ reac , L), to be discussed in section III D, Eq. (1) no longer holds, and both σ and E are heterogeneous and nonlinear function of x.

B. Reaction-dominated regime
In the rst regime, we consider systems where the chemical reaction rates are the fastest process in the system. In this regime, therefore, we expect to nd a constant surface charge and consequently the standard Helmholtz-Smoluchowski equation (1). In terms of time scales we have τ reac τ dif , τ cond , which implies that λ reac is the smallest length scale, i.e. λ reac L, λ cond . It should be noted that, since τ cond , τ dif ∝ D −1 s , a system without Stern-layer conduction (D s = 0) cannot be di usionor conduction-dominated and is in fact always in the reaction-dominated regime. e resulting equations in the reactiondominated regime do not depend on the ratio between λ cond and L, which we can therefore leave unspeci ed. In this limit, the wavenumbers can be approximated by k + = −k − = λ −1 reac . is simpli es the solution for σ(x), Eq. (12), the Stern-layer ux j σ (x), Eq. (5) (see Appendix C for the general expression), and the surface charge production rate R(x) = −n s · J C . Eq. (6), as R(x) = σ eq Eλ reac τ reac sinh x/λ reac cosh L/λ reac .
In Figs. 1(a) and 1(b) we plot σ(x) and j σ (x) respectively, for several values of L/λ reac smaller and larger than unity. From λ reac L, we can deduce from Eqs. (20)-(22) that σ(x) ≈ σ eq , j σ (x) ≈ D s σ eq E and R ≈ 0 for all x except within a few λ reac away from x = ±L. In Fig. 1, we show this deviation near x = ±L for both σ(x) and j σ (x), which is claerly visible for several values of L/λ reac . is edge e ect is a direct result of the boundary condition j σ (±L) = 0, that stems from the fact that our surface has nite length. In order for a non-zero j σ to develop, counter ions must adsorb at the inlet and desorb at the outlet, which is only possible if σ deviates from σ eq . is explains the heterogeneities of σ shown in Fig. 1 that persist even for large L/λ reac . Even in the classical Helmholtz-Smoluchowski se ing, a nite Stern-layer conduction implies that at the edges of the surface σ deviates from its equilibrium value. e range of this inhomogeneity is given by λ reac , as can be seen in Fig. 1. Consequently, the Stern-layer current and surface charge pro le are constant up to a few λ reac from the edges of the surface. e amplitude of the relative deviation is interestingly given by λ reac /λ cond . is edge e ect exists purely due to the surface charge discontinuity at x = ±L, but λ reac is nevertheless not to be confused with the healing length = HDu introduced by Khair and Squires [37], which also arises in the absence of Stern-layer conduction.
(a) (b) Figure 1: Surface charge density σ(x) (a) and Stern-layer current j σ (x) (b) between channel inlet (x = −L) and outlet (x = L), in the limit λ reac λ cond according to Eqs. (20) and (21) for varying values of L/λ reac . For (a) we used λ reac /λ cond = 0.2 (j σ does not depend on this quantity). As L/λ reac decreases, the system leaves the reaction-dominated regime and enters the di usion-dominated regime (for which λ reac λ cond holds) and σ and j σ become increasingly heterogeneous.
Imposing a vanishing net current at every position x, according to Eq. (16), we nd the streaming electric eld E, with ∂ x p a short-hand notation for ∂p ∂x . In this limit, we recover the standard Helmholtz-Smoluchowski equation (1) with a spatially constant electric eld consistent with our assumptions. However, just like j σ and σ, E is not constant close to the edges and Eq. (23) only holds several λ reac from x = ±L. e contributions of the edge e ect to the streaming potential ∆Φ = L −L dxE are negligible (note the antisymmetric nature of the edge e ect). Furthermore, we can identify using Eq. (23) the Stern-layer contribution to the surface conduction G s = D s σ eq βe 2 = G S s . e la er is proportional to the charge carrier density σ in the Stern layer and a 2D-analogue of G b = Dρβe 2 . Note that in our calculations, the EDL surface conductivity G EDL s , which originates from the increased ion density in the EDL, does not appear since we assumed that ρ is spatially constant.

C. Di usion-dominated regime
In the second dynamic regime, the di usion time τ dif is the smallest time scale, and the channel length L is the smallest length scale, L λ cond , λ reac . Consequently, this implies that the dimensionless streaming potential βe∆Φ ∼ EL = L/λ cond is small. Analogously to the Helmholtz-Smoluchowski regime, the ratio between the other two lengths, λ cond /λ reac , will have no signi cant impact on the results. We can use Eq. (20), which was derived assuming only λ reac λ cond , but now with L λ reac to write the surface charge, ux and chemical production rate as In this parameter regime, we thus recover a linear pro le for the surface charge density σ(x) also found numerically [29] and shown in Fig. 1(a). We can intuitively understand this linear pro le by realising that the system is di usion dominated. For a surface with translation invariance in one direction, the steady state would then be given by a linear pro le. e liner pro le is maintained because the chemical reaction is not fast enough to force σ to the equilibrium value (τ reac τ dif ). Additionally, since j σ ∝ τ dif /τ reac , the surface ux is very small, as can be observed in Fig. 1(b). At steady state, the surface charge pro le is therefore determined by a balance between the conduction caused by E and di usion is the opposite direction, which explains why the slope of σ(x) is given by the electric eld E. is result is analogous to a di usion-dominated convection-di usion problem. Within the di usion-dominated regime, we obtain from Eq. (16) a new expression for the streaming electric eld, Here we introduced E HS = |∂xp||ζeq| ηG b , the (magnitude of the) streaming electric eld as predicted by the Helmholtz-Smoluchowski equation (1) without Stern-layer conduction. Since we have de ned ζ eq as the dimensionless equilibrium zeta potential, E HS has dimensions of inverse length. Note that the solution does not depend on D s (except for the restriction that τ reac τ dif , τ cond ∝ D −1 s ) due to the vanishing j σ . e streaming electric eld is, similar to σ(x), heterogeneous, with a smaller value than E HS at the inlet and a larger value at the outlet. e impact of Stern-layer conduction is thus indirect in this case, and not direct via a charge current 'leaking' through the Stern layer as in the reaction-dominated regime. e analytical solution does, however, overestimate the heterogeneity of σ(x) and E(x). e reason for this is that we assumed that the counter ion concentration at the surface, ρ C,s , does not depend on the surface charge. In reality, of course, it does, since an increased surface charge a racts more counter ions. e lack of this regulation mechanism explains the overestimation of the surface charge and electric eld.
Eq. (25) additionally allows us to derive an expression for the streaming potential, where βe∆Φ HS S ≡ −z σ 2LE HS is the streaming potential predicted by the Helmholtz-Smoluchowski equation (1). For small ∆Φ, which, as we have discussed above, is always the case in the di usion-dominated regime, we nd that the streaming potential can be accurately estimated using the standard Helmholtz-Smoluchowski expression, ∆Φ ≈ ∆Φ HS . e lateral heterogeneity therefore has no signi cant e ect on ∆Φ S . is can be explained by the quasi-antisymmetric pro le of σ(x), as the streaming potential is a laterally integrated quantity.
is, combined with the tendency to measure at the center of the channel, might explain why such lateral heterogeneities have not been observed yet. Note that Eq. (26) breaks down for |∆Φ| → 1. However, we are in the regime where L λ cond , which implies that βe∆Φ ≈ 2EL 1. e streaming potential will therefore never diverge, but the system will change to a new regime as ∆Φ increases. Lastly, we can derive a surprisingly simple expression for the surface charge di erence between inlet and outlet in the di usion-dominated regime, e streaming potential therefore gives the fractional di erence in surface charge between the inlet and outlet, providing a good measure of the heterogeneity of the system. Given that here the streaming potential is approximately equal to the Helmholtz-Smoluchowski result, Eq. (27) gives a priori a measure of the heterogeneity to be expected, although one should keep in mind that in general our results overestimate the actual heterogeneity. D.
e conduction-dominated regime e third regime, the conduction-dominated regime, is reached when τ cond (λ cond ) is the smallest time (length) scale, τ cond τ dif , τ reac (and hence λ cond L, λ reac ), such that the wavenumbers can be approximated as k + ≈ E s and k − ≈ 0. Also here, the ratio of the remaining lengths does not impact the results. is allows us to simplify Eq. (12) signi cantly. e surface charge pro le is no longer linear or anti-symmetric compared to the equilibrium value σ eq , but rather exponential, while j σ vanishes It should be noticed that j σ only vanishes because we assumed a constant E. However, we have already seen that this is no longer generally the case, and j σ will in fact not exactly vanish in a fully self-consistent analysis, but our analysis does show that j σ is small. Our numerical calculations con rm that in this regime, as well as in the di usion-dominated regime, j σ is negligible compared to the bulk uxes [29] as the chemical reaction rates are too small for a signi cant surface ux to develop. Note that Eq. (28) shows that the density pro le is exponential, which is analogous to a convection-dominated convectiondi usion problem. If we take L λ cond (EL 1), we recover the same linear pro le as in the di usion-dominated regime Eq. (24).
To determine E in the conduction-limited regime we again impose a vanishing net charge current for every x, Eq. (16), from which we obtain the condition Eq. (29) can be solved numerically for E(x), and we obtain qualitatively similar behaviour as in the full numerical calculations [29] which is of course inconsistent with our assumptions that E is spatially constant. In Fig. 2 we plot the resulting surface charge pro le according to Eqs. (28) and (29) for several values of the pressure drop across the channel. For small pressure drops, the surface charge pro le is roughly linear, corresponding to the di usion-dominated regime. For larger pressure drops, the pro le becomes highly nonlinear, increasing exponentially as x approaches the outlet position at x = L. e point where the pro le crosses the equilibrium surface charge is, contrary to the di usion-dominated regime, no longer at x = 0. Our semi-analytical results agree qualitatively but not quantitatively with the numerical calculations. In particular, the surface charge close to x = L is much larger in the semi-analytical results. e reason for this must again reside in the assumptions that ρ C,s is constant, similar to the di usion-dominated regime. ere is, in fact, also some qualitative discrepancy between the numerical calculations and the current analysis. In the numerical calculations the average surface charge decreases with an increasing pressure drop, while Fig. 2 shows that the average surface increases with an increasing pressure drops. e reason for this probably lies again in the missing charge regulation mechanism discussed above, and the exponential increase of σ(x) greatly overestimates the average surface charge. Consequently, the predicted streaming potential is also overestimated, and cannot be accurately determined in the current analytical approach. e breakdown of the theory is not unexpected in this respect, as our linearised theory and lack of regulation work best for small driving forces. Lastly, we note that multiplying Eq. (29) with L, we see that the solution is given in terms of EL rather than E, if ∆p is xed. Hence, in this regime, σ is invariant under changes in L if ∆p is xed, as was in fact also suggested by the numerical solutions [29].

IV. SUMMARY & CONCLUSION
In this work we revealed some consequences of a chemically and physically dynamic Stern layer on electrokinetic phenomena. By allowing the surface charges to di use across the Stern layer and by assigning a nite rate to the adsorption and desorption reactions, we showed that a simple electrokinetic system develops novel properties. We identi ed three dynamical regimes, schematically represented in Fig. 3. ese regimes can be identi ed via three time scales (or equivalent lengths): the chemical reaction time scale τ reac = (k des + k ads ρ C,s ) −1 , the di usive time scale τ dif = L 2 /D s and the conductive time scale τ cond = 1/(D s E 2 ). Here, E = βe∂ x ψ s is the dimensionless streaming electric eld and has dimensions of inverse length. e particular regime is determined by which of the three time scales is the smaller one. ere are three basic parameters to Figure 3: Schematic representation of the three regimes and how to transition between them. Each regime is given for each of the 6 possible orderings of the three governing time scales τ cond , τ dif and τ reac .
change the regime of the electrokinetic system: the pressure gradient ∂ x p changes the conduction time τ cond and λ cond , the bulk counter ion density (for example, the pH of the solution) alters τ reac and λ reac , while τ dif varies with L. Note that the bulk ionic strength ρ is also an experimentally tunable parameter, which has a similar e ect as the pressure gradient; increasing ρ decreases E and thus increases τ cond and λ cond . In order to transition between the three regimes, the ordering of the time scales must be changed, which is schematically represented in Fig. 3. e rst regime, the reaction-dominated regime, emerges when the chemical reaction rates are the fastest process, τ reac τ dif , τ cond . Due to the high chemical rates we nd a mostly constant surface charge pro le, except for a region of size λ reac away from the ends of the surface. Moreover, we recover the standard Helmholtz-Smoluchowski equation (1) for the streaming electric eld, from which we can read o an explicit expression for the Stern-layer conductivity. e system is in the second regime, the di usion-dominated regime, if L is the smallest length scale, or equivalently if τ dif is the smallest time scale. In this regime, the system develops a linear surface charge pro le equal to the equilibrium value at the middle of the surface. Although the surface charge and thus the streaming electric eld are laterally heterogeneous, the resulting streaming potential is within a good approximation equal to the Helmholtz-Smoluchowski expression Eq. (1) with zero Stern-layer conductance. For large streaming potential and slow reaction rates, τ cond τ dif , τ reac , the system reaches the nal regime, the conduction-dominated regime.
e surface charge pro le is then exponential, and is no longer equal to the equilibrium value in the center of the channel. Consequently, the streaming electric eld is also exponential, and must be found numerically. In the conductiondominated regime, the electric eld and thus the streaming potential di er signi cantly from the Helmholtz-Smoluchowski expression.
We believe that this theoretical framework provides both a deeper physical understanding of the processes at work for a dynamic Stern layer, and is able to predict what properties to expect from an experimental setup. A di culty is the chemical rates, which seem to be unknown for many surface-electrolyte combinations. Perhaps, by exploring the behaviour of the surface charge by altering the pressure drop or system size, our framework might actually provide insight in the numerical values of the chemical rates.
is work is part of the D-ITP consortium, a program of the Netherlands Organisation for Scienti c Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW). is work is a part of an Industrial Partnership Program of the Netherlands Organization for Scienti c Research (NWO) through FOM Concept agreement FOM-15-0521. Financial support was provided through the Exploratory Research (ExploRe) programme of BP plc.
where R is the production rate of surface charges and j σ the (2D) ux of surface charges and ∇ S is the (2D) divergence in S. For example, for a at plate in the xy-plane, we have ∇ S = (∂ x , ∂ y ). We have implemented a type B dynamic model because the total number of ions (on surface plus in the water) is conserved. From Eq. (A.3) we can write down R(r s ) R(r s , t) = −k ads (Γ − σ(r s , t)) + k des ρ C (r s , t)σ(r s , t). (A.7) Analogously to the bulk equation, the surface current is given by, where λ B = βe 2 4π is the Bjerrum length with the permi ivity, and the function f S (r) (with dimension inverse length) encodes the location of the chargeable surface S. In the case of a chargeable plate parallel to the xy-plane at z = 0, we have f S (r) = δ(z). We couple the bulk ions to the surface charges via a Robin boundary condition in R, such that the total in ux of counter ions is equal to the destruction rate of the surface charges −R. In addition to the standard electric boundary condition and the no-slip boundary condition for u, we obtain the boundary conditions −n s · J C (r s , t) = −R(r s , t), n s · ∇φ(r s , t) = −4πλ B σ(r s , t) u(r s , t) = 0 (A. 10) with n s an inward pointing normal vector (into the uid) and J C the counter ion ux. All other, non-charged surfaces are impermeable for all ions. e charged surface at S is impermeable for all ions expect the counter ion. In order we check the consistency of the governing equations, we consider equilibrium condition where all time derivatives and uxes vanish. erefore, we nd that in equilibrium R(r s , t) = 0, and Eq. (A.7) reduces to the Langmuir adsorption isotherm, σ(r s , t → ∞) = Γ 1 + ρ C (r s , t → ∞) K In order for all bulk uxes to vanish, the functional derivative of F must reduce to a constant, δβF b δρ i (r) = constant = µ R , (A. 12) with µ R the chemical potential of the system (there is no external potential except for the hard wall potential) xed by ρ b,i , the bulk concentrations of the ions. Solving for the densities ρ i we get the Boltzmann distribution for all ionic species, and therefore represent a valid equilibrium condition. Alternatively, we could also have enforced j σ (r s ) = 0, and thus a constant rst derivative of F s , which can then be solved for the equilibrium surface charge σ(r s ) = Γ 1 + C σ e zσφ(rs,t→∞) −1 , r s ∈ S, (A. 15) with C σ an integration constant. We can then combine this with the condition R = 0 to obtain once again the Langmuir adsorption isotherm for the surface charge and the Boltzmann weight for the dissolved ions. is shows the internal consistency of the theory, for if we set any 2 of J i , j σ or R to zero, it follows that the third vanishes.
system we thus demand that z σ j σ (x) + where z σ j σ (x) is the net charge current through the Stern layer. Using Eq. (8), and the solution to ρ e , Eq. (9), we nd the net charge current through the liquid where have used that H λ D and included a factor z σ in the second term on the right hand side since E was de ned as a positive quantity. Lastly, we use Eq. (12) and (5)