Adelchi J Asta

Maximilien Levesque

Benjamin Rotenberg

Moment propagation method for the dynamics of charged adsorbing/desorbing species at solid-liquid interfaces

Moment propagation method for the dynamics of charged adsorbing/desorbing species at solid-liquid interfaces

INTRODUCTION

The dynamics of ions at charged solid/liquid interfaces, in particular in porous media, results from the complex interplay between diffusion, advection by the fluid flow and electrostatic interactions. In turn, the motion of ions also influences the local charge distribution, hence the electrical force acting on the fluid. Such electrokinetic effects, which are encountered from energy conversion exploiting couplings on the scale of interfaces [1,2] to large scales applications in Earth Sciences [3] have been the subject of extensive theoretical and numerical studies (see e.g. Refs [4][5][6][7][8][9][10]) and a number of simulation tools have been proposed to model them at various levels, from the molecular to the macroscopic ones (see e.g. Refs. [11,12] for reviews).

It has also become obvious recently that the behaviour of charged solid/liquid interfaces does not only depend on the electrostatic interactions between the ions and the wall, but also on their chemical nature, which modifies their short-range interactions, e.g. via the solvation properties of the ions. [13,14] Molecular simulations have clarified such ion-specific effects, as well as the consequences on the electrokinetic effects on both hydrophilic and hydrophobic surfaces. [15] Details of the electronic properties of the walls also have strong effects on the interaction with solutes: Density Functional Theory (DFT) calculations suggest for example that differences in the interactions between the hydroxide ions and either carbon or boron nitride nanotubes, which have the same geometry but different electronic structures [16], may be at the origin of the dramatically different hydrodynamic behaviour observed on larger scale for the flow of water through these tubes. [17,18] At a more coarse-grained level, specific effects can be introduced in mean-field theories such as Poisson-Boltzmann via Potentials of Mean Force (PMF) which can be computed using molecular simulations. [19] Even without electrokinetic couplings, specific effects play an important role on the dynamics of charged species. For example, the diffusion of ions in porous materials depends not only on their charge but also on their chemical nature. An illustration can be found e.g. with the different properties of Na + and Cs + tracers in clays, which go beyond their different behaviour in the bulk. [20] Of course, charged tracers also experience the effects of hydrodynamics and electrokinetic couplings in the fluid as a whole, in addition to their own dynamics and interactions with the solid surfaces. This is exploited in practice in analytical chemistry e.g. in chromatography or electrophoresis experiments. [21] It is well known e.g. in the chemical engineering community that even in the case of uncharged solids and solutes, the coupling between pore-scale motion and surface adsorption may result in intricate macroscopic transport properties [START_REF] Brenner | Macrotransport Processes Butterworth-Heinemann Series I[END_REF]. Surface charge and charged solutes only increase the difficulty to understand and predict the emerging behaviour.

At the coarse-grained or macroscopic levels, specific interactions with the surfaces are usually introduced via adsorption and desorption reactions, with corresponding rates. Such reactions are usually assumed to be of firstorder, even though in principle the description can be improved to account for more complex features such as the saturation of surface sites. In the present work, we describe a numerical scheme which allows the simulation of mobile tracers evolving under the combined effects of diffusion, advection by the flow, migration under local electric field and adsorption/desorption with the surface.

We combine previous developments for the separate description of these phenomena, in particular taking advantage of the Lattice Boltzmann Electrokinetics (LBE) method to capture electrokinetic effects in the underlying fluid. [START_REF] Capuani | [END_REF][24][25][26][27][28][29] The Moment Propagation (MP) approach [30][31][32][33][34][35][36] is not limited to the mere computation of average fluxes but allows the computation of observables which reflect the statistical properties, averaged over all possible trajectories, of their dynamics, such as their velocity auto-correlation function (VACF).

Section presents the theoretical basis for the description of the transport of solvent and ions as well as their adsorption/desorption at solid/liquid interfaces. The lattice-based algorithms are then presented in Section . Section finally considers the case of dispersion by an electro-osmotic flow in a slit-pore with charged walls and counterions in the absence of added salt, for which some properties can be determined analytically. The algorithm is validated by examining the average velocity of tracers as a function of surface charge, of the tracer valency and of the adsorption/desorption rates. We then illustrate the ability of the MP method to estimate properties for which no analytical results are available, by considering their dispersion coefficient.

ELECTROKINETIC AND ADSORPTION/DESORPTION EQUATIONS

The dynamics of ions in a solvent is governed by several processes, such as diffusion, migration under the effect of the local or an external electric field and advection by the flowing fluid. These phenomena can be encompassed in a conservation equation for the ionic densities ρ k as a function of position r and time t:

∂ t ρ k (r, t)+∇•[ρ k (r, t)u(r, t) + βD k q k eρ k (r, t)E + j k (r, t)] = 0
(1) where u is the local velocity of the fluid, β = 1/k B T with k B Boltzmann's constant and T the temperature, D k is the diffusion coefficient of species k, q k its valency, e the elementary charge and E the applied electric field. Finally, j k is the flux in response to local thermodynamic forces, which can be written for dilute enough solutions as:

j k (r, t) = -βD k ρ k (r, t)∇µ k (r, t) (2)
where µ k the local chemical potential. Separating the ideal and excess contributions,

µ k = k B T ln(ρ k /ρ • k )+µ ex k , the fluxes simplify to j k (r, t) = -D k ∇ρ k (r, t) -βD k ρ k (r, t)∇µ ex k (r, t) (3)
Finally, if solutes interact only via electrostatic interactions, at the mean-field level the excess chemical potential is simply the product of the charge q k e and ψ the local electrostatic potential, solution of the Poisson equation:

∆ψ(r, t) = - e 0 r k q k ρ k (r, t) (4)
with 0 the vacuum permittivity and r the dielectric constant of the solvent. Summarizing, and dropping the explicit space-and time-dependence, this leads to the Nernst-Planck equation:

∂ t ρ k +∇•[ρ k u -D k ∇ρ k + βD k q k eρ k (E -∇ψ)] = 0 (5)
Finally, the dynamics of the fluid is governed by the Navier-Stokes equation for the velocity u(r, t):

∂ t u + u • ∇u = η∆u - k ρ k ∇µ k + f ext V (6)
with η the fluid viscosity, f ext V the external force density and where the first term in the right hand side corresponds to the thermodynamic force, which can be expressed as a pressure gradient using the Gibbs-Duhem equation. Separating again the ideal and excess parts in the chemical potential, in the case where µ ex k = q k eψ, and collecting the latter with the contribution of the external force, we recover the usual electric force density e ( k q k ρ k ) (E -∇ψ).

The combination of the Poisson (4), Nernst-Planck (5) and Navier-Stokes ( 6) equations provides a simple yet reasonably accurate description of the coupled ionic and solvent dynamics, provided that the above-mentioned assumptions apply. They must be supplemented by electrostatic and hydrodynamic boundary conditions. We will focus here on fixed surface charge density σe (with σ per unit surface), which may in principle depend on the position on the surface, and no-slip (u = 0) at the fluid-solid interface.

We will further investigate the dynamics of charged tracers within such an ionic solution. Tracers follow the same evolution equation (5) as the major ions, i.e. experience diffusion, advection and migration, but they do not influence in return the electrostatic potential ψ or the fluid velocity u (or the other ionic densities). However, we will also consider the adsorption and desorption of such tracers at the surface of solid walls. Introducing Γ k the surface concentration of species k (length -2 ) and assuming first-order kinetics for both processes, the evolution of the adsorbed and mobile tracers at the surfaces follow:

∂ t Γ k (r, t) = -k d Γ k (r, t) + k a ρ k (r, t) (7)
with k d (time -1 ) and k a (length•time -1 ) are the desorption and adsorption rates, respectively. In equilibrium, this results in the Henry law adsorption isotherm. These coupled non-linear equations, namely (4-6) for the dynamics of the major species and of the fluid and ( 5)+ (7) for the adsorbing/desorbing charged tracers should then be solved, in general numerically, to predict the evolution of the system. At this mean-field level, however, such a solution does not provide information on the dynamics of individual tracers and how they explore the fluid and the interface under the combine effects of all microscopic processes. Here we extend the Moment Propagation method, which allows the computation of observables which reflect the statistical properties, averaged over all possible trajectories, of their dynamics, such as their velocity auto-correlation function (VACF).

ALGORITHMS Lattice-Boltzmann Electrokinetics

Previous work, cited in the introduction and below, have demonstrated the power of lattice-based method to achieve both goals, namely solving the electrokinetic equations as well as computing statistical properties over individual trajectories. Central to these descriptions are: (a) the probability distribution f (r, v, t) for finding a particle with a given velocity at a given position and a given time; (b) the spatio-temporal discretization of the evolution equation for this quantity, which is also associated with a discretization of the velocity-space. In a nutshell, this requires introducing a lattice of spacing ∆x and a time step ∆t together with a set of discrete velocities {c i } connecting nodes of the lattice. The probability distribution is discretized accordingly into populations f i (r, t), which evolve following a discretized kinetic equation of the form:

f i (r + c i r, t + ∆t) -f i (r, t) = ∆ i (r, t) . (8)
In this equation, the left-hand side corresponds to the transport of particles to neighbouring nodes according to their discrete velocity, while the right-hand side (collision operator) accounts for the redistribution of their velocities due to interactions within the fluid. This algorithm allows to propagate the distribution function f , from which the hydrodynamic observables (local mass density, momentum and stress tensor) are computed as its moments in velocity-space. The collision operator captures interactions within the fluids and the effect of the local force on the fluid, including the thermodynamic forces in Eq. ( 6) which depends on the local concentration of each species. The latter are updated on the same lattice using the so-called link-flux method of Capuani et al. [START_REF] Capuani | [END_REF], resulting in an evolution of the concentrations according to Eq. ( 5), coupled with that of the momentum of the fluid according to Eq. ( 6) via Eq. ( 8). Finally, the Poisson equation ( 4) is solved numerically on the same lattice using the Successive Over Relaxation method [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF]. Finally, no-slip boundary conditions at the fluid-solid interface are implemented by the bounce-back rule [START_REF] Ladd | [END_REF] and solute fluxes through links between fluid and solid nodes are set to zero to enforce the no-flux boundary condition. This combination of LB and link-flux methods has already been used succesfully to simulated coupled electrokinetic effects in colloidal systems [START_REF] Capuani | [END_REF][24][25][26], channels and porous media [36,39,40] or dielectric droplets under an electric field [41]. In the following we will consider specifically the simple case of the electro-osmotic flow in a slit pore with charged walls, where the fluid is driven by an applied electric field, for which some properties can be computed analytically, offering the possibility to validate the numerical results.

Moment Propagation

In order to compute the VACF of tracers within such electrokinetic flows, we take advantage of the probabilistic description underlying the LB method via the Moment Propagation (MP) approch introduced by Lowe and Frenkel [30,42]. Several descriptions and extensions of this method have since been proposed for various applications, see e.g. [33,35,36,[43][44][45][46]. In the present work, we show that this method can be used for charged mobile (hence experiencing diffusion, advection and migration) and adsorbing/desorbing species -a combination of features which had to date not been investigated previously despite its relevance in the many contexts described in the introduction. We first show how the MP method can be implemented to propagate any quantity related to the transport of tracers under the combined effect of all the above-mentioned processes. Then, we show how a particular choice of propagated quantity can be made to compute the VACF, before finally expressing averaged quantities such as the average velocity or dispersion coefficient of tracers from their VACF.

We consider here the propagation of an arbitrary quantity P (r, t) defined on the same lattice as the one used for the LB/link-flux simulations. In the absence of adsorption/desorption processes, it is updated in steps according to P (r, t + ∆t) = P * (r, t + ∆t) with:

P * (r, t + ∆t) = P (r, t) 1 - i p i (r, t) + i P (r -c i ∆t, t) p i (r -c i ∆t, t) (9)
where the sums run over discrete velocities and p i (r, t) stands for the probability of leaving a node r with velocity c i . The first term therfore corresponds to particles that have not left r between t and t + ∆t. All the dynamics is then encompassed in the definition of these probabilities, which are designed to capture the relevant physical processes [36]: (10) with w i the weight associated with velocity c i in the underlying LB lattice [47] and ∆µ ex i (r, t) = µ ex (r + c i ∆t, t)µ ex (r, t) if r + c i ∆t belongs to the fluid phase, while p i (r, t) = 0 otherwise. The first term accounts for advection by the fluid, whereas the term between curly brackets is specific to the tracer and includes both the effect of diffusion and migration under the influence of the external and internal electric fields (at the present level of description, µ ex = qψ with q the charge of the tracer and ψ the potential arising from the surface charges and major ions). More precisely, the advective contribution is computed from the LB populations as:

p i (r, t) = p adv i (r, t) + λw i 1 4 βqE • c i ∆t + 1 1 + e β∆µ ex i (r,t)
p adv i (r, t) = f i (r, t) ρ f (r, t) -w i (11)
where ρ f (r, t) = i w i f i (r, t) is the local density of the fluid, and the λ parameter in Eq. ( 10) is related to the diffusion coefficient D of the tracer via:

λ = 4D c 2 s ∆t (12)
with c s the speed of sound associated with the LB lattice. As a result, any property P will be propagated following the same dynamics as Eq. ( 5). We now complement this algorithm with the approach of Ref. [43] for adsorption/desorption at the solid-liquid interface. To that end, we introduce for interfacial fluid nodes a propagated quantity P ads (r, t) associated with adsorbed particles. It is updated according to: P ads (r, t + ∆t) = P ads (r, t)(1p d ) + P (r, t)p a (13) where p a = k a ∆t/∆x and p d = k d ∆t, with k a and k d the rates defined in Eq. (7). Conversely, the evolution of the propagated quantity for mobile species is modified to: P (r, t+∆t) = P * (r, t+∆t)-P (r, t)p a +P ads (r, t)p d (14) with P * defined in Eq. (9).

The combination of Eqs. (9-13) therefore corresponds to the evolution of particles according to the coupled diffusion-advection-migration and adsorption/desorption equations ( 5) and (7) for tracers. However, one can propagate properties beyond the mere density of particles. Indeed, a proper choice of P (and corresponding P ads ), defined by the initialization discussed below, allows the computation of the VACF. We briefly recall here the derivation of Ref. [36]. The starting point is the definition of the VACF for the γ ∈ {x, y, z} component of the velocity, which in the present case where the velocities of particles can only assume discrete values:

Z γ (t) = v 0 γ v t γ = r 0 ,c 0 γ π(r 0 , c 0 γ )π(c t γ |r 0 , c 0 γ )c t γ c 0 γ = r 0 ,c 0 γ π(r 0 , c 0 γ )c 0 γ r t π(r t , c t γ |r 0 , c 0 γ )c t γ ( 15
)
where the superscripts refer to times 0 and t and the discrete sums run over nodes and set of velocities associated with the chosen lattice. Moreover, π(r 0 , c 0 γ ) is the probability of being at node r 0 , with velocity c 0 γ and π(c t γ |r 0 , c 0 γ ) the probability to have a velocity c t γ , given that the particle was initially at node r 0 with velocity c 0 γ (and similarly for the joint distribution in the sum over r t ). Eq. 15 can be rewritten by replacing π(r t , c t γ |r 0 , c 0 γ ) by the post-collisional average π(r t , c t γ |r 0 , c 0 γ ) i p i (r, t)c iγ = π(r t , c t γ |r 0 , c 0 γ )u * γ (r, t), with the local average tracer velocity defined by:

u * (r, t) = i p i (r, t)c i ( 16
)
which to first order in β||c i • ∇µ ex ||∆t reduces to u * ≈ u + βD(qeE -∇µ ex ) = u + βDqe(E -∇ψ). We can then rewrite the VACF as a sum over all lattice nodes at time t:

Z γ (t) = r t   r 0 ,c 0 γ π(r 0 , c 0 γ )c 0 γ π(r t |r 0 , c 0 γ )   u * γ (r, t) = r P (r, t, γ)u * γ (r, t) (17)
which also defines the relevant probabilty P (r, t, γ) to be propagated for each component γ of the VACF, namely the probability to arrive at node r at time t, weighted by the initial velocity of the particle. Since the particles adsorbed at the solid-liquid interface are considered as immobile (i.e. we neglect here surface diffusion), they do not enter directly in u * . However, the adsorption and desorption processes do contribute to the VACF via the evolution of P in Eq. 17, which is coupled to that of the corresponding P ads -see Eqs. ( 13) and ( 14). The crucial step is then the definition of the initial values of both propagated quantities. Specifically, this is achieved by the following choice:

P (r, 1, γ) = i e -βµ ex (r-ci∆t) Q p i (r -c i ∆t, 0)c iγ (18)
for all fluid (F ) nodes, including interfacial (I) ones, with Q the partition function of the tracers, which also includes the adsorbed ones:

Q = r∈F \I e -βµ ex (r) + r∈I e -βµ ex (r) 1 + e -β∆µ ads (r) (19)
where e -β∆µ ads = k a /k d ∆x defines the tracer sorption free energy. We further initialize the corresponding quantity for adsorbed tracers to 0, because of their vanishing velocity.

The initial value of the VACF is simply given by

Z γ (0) = i e -βµ ex (r) Q p i (r, 0)c iγ (20)
and subsequent values are computed with Eq. ( 17) where the propagated quantity is initialized via (18) and evolved according to Eqs. ( 9), ( 13) and ( 14). This completes the description of the algorithm. Note that this should be applied separately for each direction γ and of course for each tracer, defined in the present case by its valency q (which enters in µ ex ), its diffusion coefficient D (which enters in the transition probabilities p i ) and the adsorption/desorption rates k a and k d .

All the above algorithms have been implemented in the LABOETIE code.[48]

RESULTS AND DISCUSSION

Dispersion of charged adsorbing tracers by an electro-osmotic flow

In order to illustrate the ability of the moment propagation method to capture the combined effects of adsorption/desorption, in addition to that of advection, electromigration and diffusion, we consider here the simple case on which it had previously been tested in the absence of adsorption/desorption. Specifically, we investigate the dispersion of tracers with valency q ∈ {-1, 0, +1} by an electro-osmotic flow in a slit pore with parallel walls of surface charge density σe < 0 separated by a distance L, with monovalent counterions, as illustrated in Figure 1. This system has the advan-FIG. 1. Electro-osmotic flow in a slit pore between negatively charged walls.

tage, for validation purposes, that a number of properties can be computed analytically in the absence of added salt. In particular, the electrostatic potential in the direction perpendicular to the surfaces is given by: ψ(x) = ψ(L/2) + k B T e ln cos 2 (αx)/ cos 2 (αL/2) , with x ∈ [-L/2, +L/2] and where the characteristic length α -1 is the solution of αL 2 tan αL 2 = πσLl B with the Bjerrum length l B = e 2 /4π 0 r k B T . In the presence of an applied electric field E y parallel to the surfaces, the steadystate profile is governed by the balance between the electrostatic force eρ + (x)E y (since only monovalent cations are present in the fluid) and viscous force ηu y (x), resulting in the electro-osmotic flow:

u y (x) = eE y 2πηl B ln cos(αx) cos(αL/2) ≡ u ref h(x) (21)
which also defines the reference velocity u ref = eE y /2πηl B and the scaling function h. For small surface charge density (αL → 0), the velocity profile is almost parabolic, u y (x) ∼ u ref

(αL) 2 8
1 -4x 2 L 2 , as in the case of a Poiseuille flow (applied pressure gradient instead of electric field, which corresponds of course to a different prefactor).

Note that in the above the counterions are not subject to adsorption/desorption, which only applies to the tracers. The fraction of adsorbed tracers f ads (or equivalently of mobile tracers f mob = 1f ads ) depends on the surface charge and distance between surfaces via the product αL, as well as on their valency q and the ratio k a /k d L. More generally, the equilibrium distribution of tracers within the pore is given by the normalized Boltzmann weights:

B q (ξ) = B * q (ξ) 1/2 -1/2 B * q (ξ) dξ + ka k d L B * q ( 1 2 ) + B * q (-1 2 ) ( 22
)
with ξ = x/L ∈ [-1/2, +1/2] and

B * q (ξ) = e -qeψ(ξ)/k B T = cos(αLξ) cos(αL/2) -2q (23)
The integral in the denominator of Eq. ( 22) can be performed analytically for q ∈ {-1, 0, +1}, with the results 1 + sin αL αL / (1 + cos αL), 1 and sin αL αL respectively. The variation of the fraction of adsorbed tracers with αL and k a /k d L will be discussed below.

From the dynamical point of view, we will analyze two properties relevant in practice for the transport of tracers, namely the average velocity and the dispersion coefficient, in the direction of the flow:

vy,q = 1/2 -1/2 B q (ξ) [βDqE y + u ref h(ξ)] dξ (24)
and

D y,q = ∞ 0 [v y,q (
0)vy,q ] [v y,q (t)vy,q ] dt (25)

Both can be determined from the VACF (computed in the following with the moment propagation method) as: vy,q = Z y (∞) and:

D y,q = ∞ 0 [Z y,q (t) -Z y,q (∞)] dt (26)
with Z y,q (∞) = lim t→∞ Z y,q (t).

Simulation parameters

We use the standard D3Q19 lattice, a cubic lattice with 19 speeds connecting each node to the nearest and next-nearest neighbouring nodes (and zero velocity). The associated speed of sound is c 2 s = 1 ∆t . The diffusion coefficient of the counterions and all the tracers is D = 0.05 ∆x 2 ∆t , which ensures that the Schmidt number Sc = ν/D is large, as in the case for small ions in water. In order to resolve the variations of the electrostatic potential, ionic concentrations and velocities over a distance α -1 , we use a lattice spacing ∆x = 2.5 l B . Since the Bjerrum length in water at room temperature is l B ≈ 7 Å, the distance between walls is L = 100∆x ≈ 175 nm. Simulations are performed for 4 surface charge densities, namely 2N x N y σ∆x 2 = -0.1, -1.0, -2.0 and -5.0 corresponding to αL = 0.96, 2.29, 2.63 and 2.91. For each case, the electric field is applied in the y direction with magnitude βeE y ∆x = 0.0 to 0.15 in reduced units.

Finally, once the steady-state is reached in each case with LBE, the populations f i are used in the moment propagation for 3 different tracers with valency q ∈ {-1, 0, +1} to obtain the corresponding VACFs, from which the average velocity vy,q and dispersion coefficient D y,q are obtained as explained above. For each tracer, 4 different simulations to investigate the effect of adsorption/desorption. In addition to the reference case without adsorption (k a ∆t/∆x = 0), we consider a finite adsorption rate k a ∆t/∆x = 10 -1 and three desorption rates k d ∆t = 10 -2 , 10 -3 and 10 -4 .

Fraction of adsorbed tracers

Before discussing the dynamical properties, we first summarize the equlibrium fraction of adsorbed tracers for the considered systems in Figure 2, which is calculated from the analytical expression Eq. [START_REF] Brenner | Macrotransport Processes Butterworth-Heinemann Series I[END_REF]. For given FIG. 2. Fraction of adsorbed tracers f ads with valency q = +1 (top), 0 (middle) and -1 (bottom), as a function of αL which quantifies the strength of electrostatic interactions with the walls (see section ). In each panel, we consider a finite adsorption rate ka∆t/∆x = 10 -1 and three desorption rates: k d ∆t = 10 -2 (red), 10 -3 (blue) and 10 -4 (green). In the last case simulated in the present work (without adsorption), obviously f ads = 0 regardless of the tracer valency (not shown).

The lines are computed from Eq. ( 22) and the symbols indicate the values of αL corresponding to the simulated systems.

electrostatic conditions (fixed αL), f ads increases with increasing k a /k d for all tracers. For fixed k a /k d , the variation with of f ads with the electrostatic conditions depends on the valency of the tracer: f ads does not depend on αL for neutral tracers, while for the present case of a negatively charged surface it increases (resp. decreases) with αL for cations (resp. anions). This is a direct consequence of the effect of the surface charge density on the concentration of mobile tracers at the surface, which are in equilibrium with the adsorbed ones.

Electro-osmotic flow profile

Figure 3a reports the electro-osmotic flow profile for an applied electric field βeE y ∆x = 0.1 and a surface charge density 2N x N y σ∆x 2 = -2.0. Together with the above simulation parameters, this corresponds to αL = 2.63. As expected from previous work on a nearly identical system [36], the LBE results are in excellent agreement with the analytical solution Eq. ( 21). The fluid velocity is maximal near the center and more inhomogeneous near the walls (where it vanishes). Figure 3b then shows the variation of the average velocity:

ū = 1 L L/2 -L/2 u y (x) dx ( 27
)
with u y (x) given by Eq. 21, as a function of αL. As expected, the average electro-osmotic flow increases with increasing surface charge density. In particular, it vanishes as u ref

(αL) 2 12
in the limit of small surface charges (αL → 0).

Average tracer velocity

Figure 4 indicates the average velocity vy,q of tracers with valency q ∈ {-1, 0, 1} as a function of the Péclet number P e = ūL/D with ū the average fluid velocity, for a surface charge density corresponding to αL = 2.63 and adsorption/desorption rates (k a ∆t/∆x, k d ∆t) = (10 -1 , 10 -3 ). We first note that the excellent agreement between the simulation (LBE/MP) results and the analytical solution Eq. ( 24) for all tracers. Even though the average velocity depends only on the limit of the VACF at long times, this provides a first validation of the proposed MP algorithm combining electrokinetics and adsorption/desorption. We further note that the latter decreases the average velocity, by an amount which depends on the charge of the tracer. More precisely, it follows from Eqs. ( 22) and ( 24) that the average velocity is simply equal to the product of the average velocity without adsorption/desorption (previously studied in Ref. [36]) and the fraction of mobile tracers, f mob = 1f ads . The proportionality between the average tracer velocity and average fluid velocity suggests to investigate how the ratio vy,q /ū depends on the surface charge and the adsorption properties. This is illustrated in Figure 5. The good agreement with analytical results in all cases (varying αL, q, k a and k d ) further validates the present MP scheme.

In the simpler case of neutral tracers, the average tracer velocity without adsorption is equal to the average fluid velocity. In the presence of adsorption, the ratio vy,0 /ū is independent of the surface charge and simply equal to the fraction of mobile tracers, which decreases as k d increases for fixed k a . Similarly, the ratio vy,q /ū also decreases with increasing adsorption in the case of positive and negative tracers. However, the behaviour is more complicated since the mobile fraction as well as the distribution of charged tracers within the pores (hence the flow) depends on the surface charge. For positive tracers, which are located closer to the surface where the fluid is slower than near the center of the pore, the ratio vy,+1 /ū decreases with increasing surface charge (αL). This ratio is larger than 1 for small surface charge, because the motion of the counterions (q = +1) is dominated by the direct effect of the electric field. Specifically, in this regime vy,+1 ∼ βDeE y , while ū ∼ eEy 2πηl B (αL) 2 so that their ratio diverges as (αL) -2 . It gradually vanishes as the surface charge increases, because the motion becomes dominated by the advection by the electro-osmotic flow, which is more efficient for the fluid as a whole than for the cations located on average closer to the surface. ber P e = ūL/D with ū the average fluid velocity, for tracers with valency q ∈ {-1, 0, 1}, for a surface charge density corresponding to αL = 2.63 and adsorption/desorption rates (ka∆t/∆x, k d ∆t) = (10 -1 , 10 -3 ). The simulation results (symbols) are compared to the analytical results Eq. ( 24) (solid line). In each case, the analytical solution without adsorption (dashed line).

In contrast, co-ions (q = -1) are located on average closer to the center of the pore and the direct effect of the electric field is to drive them in the direction opposite to the electro-osmotic flow. As a result, the ratio vy,-1 /ū is negative for small surface charge and diverges as -(αL) -2 for αL → 0, but is positive for large αL as the motion of co-ions becomes dominated by advection by the electro-osmotic flow. This ratio tends to 1 as αL → π regardless of adsorption, since the fraction of mobile co-ions also goes to 1 in this limit where both vy,-1 and ū are dominated by the flow near the center of the pore.

Dispersion coefficient

Due to their diffusion in the direction perpendicular to the surface, tracers experience various streamlines with different velocities. The adsorption/desorption processes also participate in the dispersion since they offer other possibilities for tracers to adopt different states (here the adsorbed species are considered as immobile, but surface diffusion would also contribute). The dispersion coefficient in the direction of the flow, D y , normalized by the diffusion coefficient D, as a function of the Péclet number ūL/D, is shown for tracers with q ∈ {-1, 0, 1}, for a surface charge density corresponding to αL = 2.63.

In this geometry, for all cases the dispersion coefficient increases quadratically with the Péclet number ūL/D, In each case, simulation results (symbols) for various adsorption/desorption rates are compared with the anaytical results Eq. ( 24) (lines).

more precisely as:

D y D = f mob + f (αL, q) × P e 2 (28)
which defines the prefactor f (αL, q) which also depends on the adsorption/desorption rates k a and k d . While in the absence of electric field (E y = 0) the dispersion coefficient is simply D y = f mob D, this is not the case in the presence of an applied field. Following De Leebeck and Sinton for the case of a cylindrical channel [START_REF] De Leebeeck | [END_REF], we obtained in Ref. [36] the expression for the dispersion coefficient for a slit pore in the no-adsorption case:

f (αL, q, k a = 0) = - 1/2 -1/2 dξ B q (ξ)g q (ξ)× ξ 0 dξ 1 B q (ξ ) ξ 0 dξ B q (ξ )g q (ξ ) (29)
where g q (ξ) = [u y (ξ) + βDqeE yvy,q ]/ū measures the local deviation from the average velocity. The simulation results reported as symbols in Figure 6 are in good agreement with this analytical result, reported as the dashed FIG. 6. Dispersion coefficient, normalized by the diffusion coefficient, Dy,qL/D as a function of the Péclet number ūL/D, for tracers with valency q ∈ {-1, 0, 1}, for a surface charge density corresponding to αL = 2.63, without adsorption (dashed lines) and with adsorption/desorption rates (ka∆t/∆x, k d ∆t) = (10 -1 , 10 -3 ) (solid lines). The simulation results are shown as symbols, while the lines correspond to the quadratic form Eq. (28). In this equation, the value at zerofield f mob = 1f ads corresponds to Figure 2; the curvature is known analytically in the no-adsorption case, see Eq. ( 29) and fitted numerically in the presence of adsorption/desorption. Results of this fitting are discussed in Figure 7.

lines. In this no-adsorption case, counterions (q = +1) are more dispersed than neutral tracers (q = 0) because they are mainly located near the walls, where the velocity profile is more inhomogeneous (see Figure 3). The opposite behaviour is observed for co-ions (q = -1), which are mainly located near the center where the flow is more homogeneous.

In the general case with adsorption, Eq. ( 29) does not apply, and the "curvature" f (αL, q) must be determined numerically by fitting the simulation results to Eq. (28). Note that in each case the only fitting parameter is f (αL, q) since f mob = 1f ads is known independently (see Figure 2). The solid lines in Figure 6 illustrate in this particular case that the behaviour remains indeed quadratic with applied field (hence P e) and that the effect of adsorption depends strongly on the valency of the tracer. For example, in this case dispersion is increased by adsorption/desorption for neutral and negative tracers, while it is decreased for positive ones. This dependence is shown in more detail in Figure 7 which reports f (αL, q) for the same adsorption/desorption rates as for the discussion of the average velocity (Figure 5).

f (αL, -1)
Analytical (no ads.)

LBE/MP (k a ∆t/∆x, k d ∆t) (0, 0) (10 -1 , 10 -2 ) (10 -1 , 10 -3 ) (10 -1 , 10 -4 )
FIG. 7. Prefactor of the dispersion vs Péclet number, f (αL, q) defined by Eq. ( 28), for q = +1 (top), 0 (middle) and -1 (bottom). In the no-adsorption case (black), simulation results are compared to the analytical expression Eq. 29 (solid line). The other simulation results are for finite adsorption rate ka∆t/∆x = 10 -1 and three desorption rates: k d ∆t = 10 -2 (red), 10 -3 (blue) and 10 -4 (green). The dashed-lines are only guide for the eyes.

It is less easy to understand than the average velocity, which boils down to the combined effects of the fraction of mobile species and their individual mobility.

In the limit of small surface charges (αL → 0) where the flow profile is almost parabolic and where the tracer distribution is (almost) flat even for charged tracers, one recovers the result for a Poiseuille flow, namely f = 1 210 . Already in this regime the effect of adsorption/desorption is not simple: Consistently with the analytical result for Taylor-Aris dispersion by a Poiseuille flow [51][52][53], f does not vary monotonically with e.g. decreasing desorption rate k d at fixed adsorption rate k a . In this case (without any effect of charge), it further depends separately on the ratios k a /k d L and D/k d L 2 [START_REF] Brenner | Macrotransport Processes Butterworth-Heinemann Series I[END_REF]. In the present case of charged tracers and electro-osmotic flows, the situation is even more complicated, even though one can note similar trends for all tracers in the limit αL → 0.

For very large surface charge densities (αL → π), f tends to decrease with increasing αL for both co-and counterions, because their distribution becomes increasingly concentrated near or away from the surface, respectively, while the electro-osmotic velocity profiles becomes flatter near the center (which also explains the less pronounced decrease of f for neutral tracers). For counterions (resp. co-ions), this is likely due to the fact that the fraction of mobile tracers then becomes very small (resp. large), see Figure 2, so that the average velocity becomes dominated by adsorbed (resp. mobile) tracers, which have the same vanishing (resp. large) velocity. In the cases with smaller desorption rates k d , this decrease with αL is observed for relatively small values in the case of co-ions (q = +1), it is preceded for counter-ions (q = -1) by an increase with αL and is almost not visible for neutral tracers. For the largest desorption rate, one observes for all tracers a maximum as a function of αL. Overall, the interplay between adsorption/desorption, migration and advection by the electro-osmotic flow results in a very rich behaviour of the dispersion coefficient -more difficult to interpret than the average mobility. This underlines the usefulness of a numerical simulation tool to investigate the combined effects of all processes.

CONCLUSION

We have extended the Moment Propagation (MP) method to capture the combined effects of adsorption/desorption of charged tracers, their migration under local and applied electric fields, as well as their advection by the local velocity of the fluid. This was achieved by combining previous developments for the separate description of these phenomena, in particular taking advantage of the Lattice Boltzmann Electrokinetics method to capture electrokinetic effects in the underlying fluid. For validation purposes, we examined the simple case of dispersion by an electro-osmotic flow in a slit-pore with charged walls and counterions in the absence of added salt, for which some properties can be determined analytically. We computed with MP the velocity autocorrelation function (VACF) of charged and neutral tracers, from which we extracted their average mobility and dispersion coefficient. Comparison with analytical results for the former allowed to validate the algorithm; the dispersion coefficient then illustrates an example of property which can be provided by the MP method when no analytical results are available. For both properties, we discussed the combined effects of the surface charge, of the tracer valency and of the adsorption/desorption rates.

Beyond the average velocity and dispersion coefficent, the VACF is also related to the time-dependent diffusion coefficient D(t). This is of particular interest in porous materials since it reflects how each tracer explores the porosity accessible to it (which depends on its valency in charged porous materials) and can be measured experimentally using NMR [START_REF] Sen | Concepts in Magnetic Resonance Part A[END_REF]. As also mentioned previously, the MP method is not limited to the VACF and can be used e.g. to compute NMR spectra of species diffusing in porous materials [START_REF] Merlet | [END_REF]. The present algorithm allows the computation of such properties for charged species undergoing adsorption/desorption, migration and advection, in simple geometries as presented here but also in more complex media. Finally, we note that we have considered here the coupling of adsorption/desorption with transport of charged tracers only, so that these processes do not influence the dynamics of the underlying fluid. While in the present case of steady-state in an infinite slit pore this would only result in the renormalization of the surface charge density, this is not true in general and a richer behaviour is expected, especially in the transient regime. In order to do this and hence describe the interplay between adsorption/desorption and electrokinetic effects, the present development of the MP method would need to be also extended to the Lattice-Boltzmann Electrokinetics scheme.

3 ∆x ∆t 2 .

 32 [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF] We use N x × N y × N z = 106 × 5 × 5 lattice points, with three layers of solid nodes on each side. The distance between the solid-liquid interfaces is thus L = (N x -6)∆x = 100∆x. Periodic boundary conditions are used in all directions. The relaxation time in the collision operator of the LB scheme is chosen as τ = ∆t (with ∆t the time step); this corresponds to a kinematic viscosity ν = 1 6 ∆x 2

  ∆t/∆x, k d ∆t) (10 -1 , 10 -2 ) (10 -1 , 10 -3 ) (10 -1 , 10 -4 )

FIG. 3 .

 3 FIG. 3. (a) Electro-osmotic flow profile in a slit pore with uniformly charged walls and no added salt (counterions only, in the presence of an applied electric field along the surfaces (see text for details). The Lattice-Boltzmann Electrokinetics simulations, normalized by the reference velocity u ref = eEy/2πηlB, are compared with the analytical solution Eq. (21). (b) Average fluid velocity ū, normalized by the reference velocity u ref , as a function of αL.

FIG. 4 .

 4 FIG.4. Average tracer velocity, adimensionalized as a Péclet number vy,qL/D as a function of the actual Péclet number P e = ūL/D with ū the average fluid velocity, for tracers with valency q ∈ {-1, 0, 1}, for a surface charge density corresponding to αL = 2.63 and adsorption/desorption rates (ka∆t/∆x, k d ∆t) = (10 -1 , 10 -3 ). The simulation results (symbols) are compared to the analytical results Eq. (24) (solid line). In each case, the analytical solution without adsorption (dashed line).

FIG. 5 .

 5 FIG.5. Average tracer velocity vy,q, normalized by the average flow velocity ū, for charged adsorbing tracers with valency q = +1 (top), 0 (middle) and -1 (bottom), as a function of the strength of electrostatic interactions with the walls (αL). In each case, simulation results (symbols) for various adsorption/desorption rates are compared with the anaytical results Eq. (24) (lines).

ACKNOWLEDGEMENTS

The authors acknowledge financial support from the French Agence Nationale de la Recherche (ANR) under grant ANR-15-CE09-0013-01. This work follows from previous thrilling collaborations with Daan Frenkel over the last decade and it is a great pleasure to dedicate the present article to him on the occasion of his 70th birthday.