Accounting for adsorption and desorption in lattice Boltzmann simulations

We report a Lattice-Boltzmann scheme that accounts for adsorption and desorption in the calculation of mesoscale dynamical properties of tracers in media of arbitrary complexity. Lattice Boltzmann simulations made it possible to solve numerically the coupled Navier-Stokes equations of fluid dynamics and Nernst-Planck equations of electrokinetics in complex, heterogeneous media. Associated to the moment propagation scheme, it became possible to extract the effective diffusion and dispersion coefficients of tracers, or solutes, of any charge, e.g. in porous media. Nevertheless, the dynamical properties of tracers depend on the tracer-surface affinity, which is not purely electrostatic, but also includes a species-specific contribution. In order to capture this important feature, we introduce specific adsorption and desorption processes in a Lattice-Boltzmann scheme through a modified moment propagation algorithm, in which tracers may adsorb and desorb from surfaces through kinetic reaction rates. The method is validated on exact results for pure diffusion and diffusion-advection in Poiseuille flows in a simple geometry. We finally illustrate the importance of taking such processes into account on the time-dependent diffusion coefficient in a more complex porous medium.


I. INTRODUCTION
The dynamical properties of fluids in heterogeneous materials offers a great challenge and have implications in many technological and environmental contexts. Inherently multi-scale in time and space, mesoscale properties such as the diffusion or dispersion coefficient reflect the nano-to-microscopic geometry of the media and the inter-atomic interactions between flowing particles, or tracers, and surface atoms. Experimentally, information about the microstructure of porous media can be extracted from diffusion measurements by pulsed gradient spin echo nuclear magnetic resonance (PGSE-NMR) [1][2][3]. At short times, the dynamics of a pulse of tracers is connected to the geometry of the porous medium at the pore scale. At longer times, macroscale properties such as porosity and tortuosity come into play. Theoretically, stochastic approaches have been an important support to the understanding of the underlying phenomena [2][3][4], and have been used recently to show that adsorption and desorption processes may strongly modify the short and long time dynamics of the tracers [5,6].
In numerous if not all practical situations involving particle diffusion and advection, the carrier fluid is in contact with confining walls where adsorption may occur. These processes depend on the chemical nature of the so- * maximilien.levesque@gmail.com lute, which explain why particles with the same charge may diffuse in the same medium with different effective diffusion coefficients. This species-dependent affinity is at the heart of all chromatographic techniques used in analytic and separation chemistry [7]. It also plays a crucial role in the dissemination of toxic or radioactive pollutants in the environment, and conversely in remediation strategies. Recently, the great interest for nanofluidic devices and for the transport in heterogeneous porous media has also raised the issue of the relevance of models which do not take into account these sorption processes. Moreover, it was recently shown that stochastic resonance between these processes and some external field may be of practical importance, e.g. for molecular sorting [5,8].
At the mesoscale, the dynamics of particles in a fluid can be described by the continuity equation where ρ(r, t) is the one-particle density at position r and time t, ∂ t ≡ ∂/∂t and J is the particle flux, which is a function of the velocity field of the carrier fluid, the bulk diffusion coefficient D b of the particles and, if any, their charge and the local electric field arising from their environment. If adsorption is taken into account, solidliquid interfaces located at r have a surface concentration Γ (r, t) (length −2 ) that evolves with time according to: where k a (length·time −1 ) and k d (time −1 ) are kinetic adsorption and desorption rates. For molecules, the rates can vary widely. As an example, the dissociation rate of DNA double strands on a surface grafted with singlestrand DNA ranges from 10 −5 to 10 −3 s −1 for a few tens of base pairs [9]. Their adsorption rate can be adjusted by changing the grafting density. Finally, we assume that the tracers (the solutes) neither diffuse into the solid phase (even though that process can easily be accounted for using our algorithm) nor dissolve the surfaces [10]. The time-dependent diffusion coefficient D(t) and the dispersion coefficient K of the tracers can be investigated by following the spreading of a tracer pulse in the fluid. This would amount to solving Eqs. 1 and 2, e.g. with a finite element method, for all possible initial conditions, which is computationally intractable for complex systems such as heterogeneous porous media. An alternative is to deduce K and D(t) from the tracer velocity auto-correlation function (VACF) following [11][12][13][14]: where the VACF in the direction γ ∈ {x, y, z} is At long times Z γ (∞) =v 2 γ withv γ the average velocity of the flow. The issue of averaging over initial conditions in Eq. 5 can be handled elegantly and efficiently using the moment propagation method [12,15,16], which was recently extended to charged tracers [14,[17][18][19][20][21].
In order to compute the VACF of the tracers from Eq. 5, one has to keep track of their velocity. For this purpose, we use the underlying dynamics of the fluid given by Eq. 1, which does not rely on the velocity of individual particles, but on the one-particle solvent density. Moreover, the simulation of heterogeneous multiscale media requires a numerically efficient method. The lattice-Boltzmann (LB) method [22][23][24][25][26][27][28] offers a convenient framework to deal with such situations. In the LB approach, the fundamental quantity is a one-particle velocity distribution function f i (r, t) that describes the density of particles with velocity c i , typically discretized over 19 values for three-dimensional LB, at a node r, either fluid or solid, of a lattice of spacing ∆x, and at a time t discretized by steps of ∆t. The dynamics of the fluid are governed by transition probabilities of a particle moving in the fluid from one node to the neighboring ones: where ∆ i , the so-called collision operator, is the change in f i due to collisions at lattice nodes. This LB equation recovers the fluid dynamics of a liquid and the moments of the distribution function are related to the relevant hydrodynamic variables. The reader is referred to Refs. [26] and [28] for reviews of the method.

II. ALGORITHM
In order to compute the dynamical properties of tracers evolving in a fluid described by the LB algorithm, i.e. to solve Eqs. 3-5, we use the moment propagation (MP) method [11,12]. Other methods could have been used, such as the numerical resolution of the macroscopic equations or Brownian Dynamics to simulate the random walk of tracers biased by the LB flow [29]. The latter method has often been successfully used, e.g. by Boek and Venturoli [30]. Nevertheless, the MP method offers many advantages. First, MP relies on the same ground as LB. It is based on the propagation of a position and velocity distribution function [12,15], therefore offering an elegant unified approach. Secondly, MP allows for the propagation of any moment of the distribution function f i (r, t), which offers great opportunities. For example, Lowe, Frenkel and van der Hoef exploited these higher moments to compute self-dynamic structure factors [31]. The LB-MP method has been thoroughly validated by Merks for low Péclet and Reynolds numbers [16].
In the moment propagation algorithm, any quantity P (r, t) can be propagated between fluid nodes. This quantity will be modified by adsorption and desorption processes. In their absence, P (r, t + ∆t) = P ⋆ (r, t + ∆t) with: where the first sum runs over all discrete velocities connecting adjacent nodes. The probability of leaving node r along the direction c i is noted p i (r, t). The last term in Eq. 7 represents the fraction of particles that did not move from r at the previous time-step. The expression for p i , which is central in the algorithm, depends on the nature of the tracers. It is given by: where the first two terms account for advection and are obtained by coupling the tracer dynamics to that of the fluid evolving according to the LB scheme. The weigths ω i are constants depending upon the underlying LB lattice. The last term describes diffusive mass transfers. The dimensionless parameter λ determines the bulk diffusion coefficient D b = λc 2 s ∆t/4, with c s = k B T /m the sound velocity in the fluid. It also determines the mobility of tracers under the influence of chemical potential gradients (including the electrostatic contribution) which are accounted for in the Q term as described in Ref. [17]. For neutral tracers, Q = 1.
We now introduce a new propagation scheme in order to account for adsorption and desorption at the solidliquid interface. While Eq. 7 still holds for nodes r which are in the fluid but not at the interface, for the fluid interfacial nodes we define a new propagated quantity P ads (r, t) associated with adsorbed particles: P ads (r, t + ∆t) = P (r, t) p a + P ads (r, t) (1 − p d ) , (9) where p a = k a ∆t/∆x is the probability for a tracer lying at an interface to adsorb, and p d = k d ∆t is the probability for an adsorbed, immobile tracer to desorb. There is no restriction in the definition of k a and k d so that they may depend on geometrical considerations and on the local tracer density. Finally, the evolution of the propagated quantity associated with free tracers now includes a term accounting for the desorption of adsorbed particles: P (r, t + ∆t) = P ⋆ (r, t + ∆t) + P ads (r, t) p d . (10) where P ⋆ is still given by Eq. 7. In order to compute the VACF of the tracers, one propagates as P the probability to arrive at position r at time t, weighted by the initial velocity of tracers. Thus one needs to initialize, for each direction γ, a propagated quantity according to the Maxwell-Boltzmann distribution. The Boltzmann weights for solid (S), fluid (F ) and interfacial (I ⊂ F ) nodes read: where e −β∆µ ads (r) = k a / (k d ∆x) corresponds to the sorption free energy for interfacial tracers, β ≡ 1/k B T with k B the Boltzmann's constant and T the temperature, and Z is the partition function of the tracers. The excess chemical potential µ ex (r) includes in the case of tracers with charge q a mean-field electrostatic contribution qψ (r) with ψ the local electrostatic potential. The VACF is then simply given as in the no sorption case [17] by:

III. VALIDATION
In order to validate this new scheme, we compare numerical results and exact theoretical solutions of Eqs. 1 and 2. This allows to assess the validity of our method independently of experimental results. We consider the diffusion and dispersion of tracers in a slit pore, i.e. between two walls at positions x = 0 and x = L. The time-dependent diffusion coefficient D (t) in the direction normal to the wall is given by [6]: where s is the Laplace conjugate of time t, L −1 is the inverse Laplace transform, χ = s/D b and κ = χL/2. In Fig. 1, we compare the exact time-dependent diffusion coefficient D(t) of Eq. 13 with the one extracted from Lattice-Boltzmann simulations for different sorption strength. This last quantity is defined by the fraction of adsorbed tracers. Unless otherwise stated, all simulations are performed within a slit pore of width L = 100 ∆x and a bulk diffusion coefficient D b = 10 −2 ∆x 2 /∆t. These values are chosen very conservative, since the LB method is known to be efficient even in narrow slits (even for L < 10 ∆x) and for a wide range of magnitudes in the diffusion coefficients [16,25]. L and D b therefore account for a negligible part of the difference with exact results. We can thus purposely assess the effect of the new algorithm only. In Fig. 1, we report time-dependent diffusion coefficients calculated by our method for a fixed sorption rate k a = 10 −1 ∆x/∆t and decreasing desorption rates k d ∆t = 10 −2 , 10 −3 and 10 −4 resulting in an increasing fraction of adsorbed tracers f a of approximately 16 %, 66 % and 95 %. This fraction is given by [6]: of D(t), as the partial immobilization at the surface slows down the exploration of the pore. We have shown recently that for this range of parameters, this slope is given by k d /(1 + k a L/2D b ) [6]. In this illustration of diffusion between two walls, the confinement is total so that for sufficiently long times, the effective diffusion coefficient decreases exponentially with time and tends to zero. In order to quantify the error on D(t) with respect to the exact result of Eq. 13, we plot in figure 2 the relative error on the slope and origin of the linear fit of log 10 D(t) for long times, i.e. for D b t/L 2 > 0.5, as a function of k a and k d . In the whole range of k a and k d , the relative errors on the slope and origin remain under 5 % and 2 %, respectively, which is highly satisfactory. The effect of a pressure gradient has also been studied on the same system. The resulting Poiseuille flow induces Taylor-Aris dispersion [32,33] of the tracers with a dispersion coefficient K, which is known exactly in the presence of adsorption and desorption in the simple slit geometry [5]: where P e = Lv/D b is the Péclet number and y = k a /k d L.
In Fig. 3, we compare the dispersion coefficient as calculated by LB with the exact results of Eq. 15 for various sorption strengths, as a function of the Péclet number. Adsorption significantly increases the dispersion, as it slows down part of the tracers. The agreement between our scheme and exact results is excellent, even for strong adsorption and Péclet numbers above 100.
As mentioned above, the time-dependence of the diffusion coefficient is a signature of the intrinsic geometric properties of a porous medium. The simplest model of such media consists of a compact face centered cubic (fcc) lattice of spheres of radius R [34]. The porosity, i.e. the fraction of empty (or fluid) space, is 1−π/ 3 √ 2 ≈ 26%. The unit cell contains 4 octahedral cavities of radius ≈ 0.41R connected by 8 smaller tetrahedral cavities of radius ≈ 0.22R by a small channels of radius ≈ 0.15R. This fcc lattice is illustrated in Fig. 4 for a lattice parameter L = 100∆x.
We report the time-dependent diffusion coefficient for this model porous medium in Fig. 5. At t = 0, the dif-  Fig. 4, as a function of the reduced time. Several fractions of adsorbed tracers, or sorption strength, fa defined in Eq. 14, are presented: black circles, 0%, i.e., without adsorption; red squares, 30%; green upward triangles, 79%; and blue downward triangles, 97%.
fusion coefficient is again given by the fraction of free particles times their bulk diffusion coefficient. After a reduced time D b t/R 2 = 0.5, tracers have explored the whole porosity and the diffusion coefficient tends toward the effective diffusion coefficient. The time dependence is strongly influenced by adsorption/desorption, as in the slit pore case. It is thus essential to consider these phenomena when interpreting experimental measurement of effective and time-dependent diffusion coefficients.

IV. CONCLUSION
In summary, we have proposed a new scheme that accounts for adsorption and desorption in a generic Lattice-Boltzmann scheme, allowing for the calculation of mesoscale dynamical properties of tracers in media of arbitrary complexity. These processes are modelled by kinetic rates of adsorption and desorption taking place at interfacial, fluid nodes. The algorithm has been validated over a wide range of adsorption and desorption rates and Péclet numbers in the slit pore geometry where exact results are available. Finally, we have shown on a more complex porous medium that adsorption and desorption processes may not be neglected, as they strongly modify the short, intermediate and long time behaviors of the diffusion coefficient as well as the dispersion coefficient. In turn, this demonstrates that neglecting interactions with the surface in the interpretation of D(t) as a probe of the geometry of the porous medium, as measured experimentally e.g. by PSGE-NMR, may lead to incorrect conclusions. This scheme may now be used in two ways. First, one could predict the effective diffusion coefficient in complex heterogeneous media for species with known adsorption and desorption rates (from experiments or molecular simulations). Conversely, from reference measurements of the time-dependent diffusion coefficient in controlled geometries, one could extract the adsorption and desorption rates k a and k d . Moreover, in the case of diffusion in the solid stationary phase, the method would allow us to relate the relevant diffusion constant to the shape of the elution profile.
While the present method is very general and also applies in principle to the case of irreversible adsorption (k d = 0, i.e. p d = 0 in Eq. 10), such a situation is of interest only outside of equilibrium. Indeed, in that case at equilibrium all the solute is adsorbed on the surface and its VACF corresponds to the sorbed species only. The propagation scheme (Eqs. 7-10) could nevertheless be used to investigate irreversible adsorption out of equilibrium by considering the density as the propagated quantity P (instead of the one described here to compute the VACF), as was done e.g. by Warren to simulate electrokinetic phenomena [35]. As an example of practical application where (possibly irreversible) sorption is coupled to electrokinetic phenomena, we can for example mention the case of ion adsorption onto charged minerals such as clays.