Jean-Mathieu Vanson

Anne Boutin
email: anne.boutin@ens.fr

Michaela Klotz

François-Xavier Coudert
email: fx.coudert@chimie-paristech.fr

Transport and adsorption under liquid flow: the role of pore geometry

We study here the interplay between transport and adsorption in porous systems with complex geometries under a fluid flow. Using a Lattice Boltzmann scheme extended to take into account adsorption at solid/fluid interfaces, we investigate the influence of pore geometry and internal surface roughness on the efficiency of fluid flow and the adsorption of molecular species inside the pore space. We show how the occurrence of roughness on pore walls acts effectively as a modification of the solid/fluid boundary conditions, introducing slippage at the interface. We then compare three common pore geometries, namely honeycomb pores, inverse opal, and materials produced by spinodal decomposition. Finally, we quantify the influence of those three geometries on fluid transport and tracer adsorption. This opens perspectives for the optimization of materials geometries for applications in dynamic adsorption under fluid flow.

INTRODUCTION

Due to their high specific surface area, porous materials are widely used in industrial-scale processes for a broad range of applications involving surface interactions such as phase separation, gas mixture separation, or ions exchange and capture. In the liquid phase, practical applications at large scale include, for example, water decontamination and removal of pollutants such as heavy metals or radioactive ions. [START_REF] Hua | nanotechnologies for the Treatment of Water[END_REF][START_REF] Fu | [END_REF][3] Understanding at the microscopic scale the physical phenomena occurring in these materials is key to improve and optimize their working capacity. The adsorption capacity itself, namely the density of adsorption sites and their activity, is the first parameter to consider. Nevertheless, the best adsorbent would be completely useless if the topology of the material's pore space does not allow the species to move freely to the active adsorption sites, and thus transport of molecular and ionic species are also of paramount importance. Both transport and adsorption properties of porous material directly depend on the internal pore geometry. As a consequence, understanding how the geometry of porous materials impacts both transport and adsorption represents a great stake to design more and more efficient systems.

The topics of fluid transport [START_REF] Sahimi | Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approaches[END_REF] and physical adsorption [START_REF] Roque-Malherbe | Adsorption and diffusion in nanoporous materials[END_REF] in porous materials have been thoroughly investigated using computational methods in the literature. However, there exist relatively few studies demonstrating how to use numerical methods to study the coupling of fluid transport and adsorption, especially in complex or "realistic" porous materials. Of the examples available in the recent literature, some use atomistic-scale modelling, studying for example the adsorption and diffusion in mesoporous silica through molecular dynamics and Monte Carlo methods. [START_REF] Coasne | [END_REF]7] Another approach, at the other end of the scale, is to perform three-dimensional numerical studies based on stochastic models, for example to shed light on the adsorption kinetics of chromatographic packed beds. [8,9] More recently, Botan et al. proposed a bottom-up model, rooted on statistical mechanics, to upscale molecular simulation and describe adsorption and transport at larger time and length scales. [10] It is important to note that, when dealing with hierarchical porous materials which present complex pore geometries presenting multiple length scales, atomistic molecular simulation methods become computational prohibitive, because of the large system sizes necessary for an accurate representation of the pore space. While those methods perform well for nano-sized systems and can describe in detail the local phenomenon of adsorption, they do not allow to efficiently compute solute properties at a macroscopic level and in the time scale required for fluid dynamics. On the other hand, at a macroscopic level computational fluid dynamics is the focus of an entire field of research, and simulates very well the behavior of the fluid. However, they are difficult to adapt to multi-phase systems and to take into account the adsorption process in heterogeneous systems. There has thus been in recent years a focus on the development of methods aimed at modeling transport and adsorption in hierarchical porous materials of various nature and pore sizes. For an introduction to those, we refer the reader to the recent review of Coasne. [11] In this work, we use a lattice-based mesoscopic fluid simulation method, namely a Lattice Boltzmann model [START_REF] Succi | The Lattice-Boltzmann equation for fluid dynamics and bey[END_REF][START_REF] Ladd | [END_REF] recently expanded to take into account adsorption, [14,15] to investigate the effects of pore geometry on the adsorption and transport of species in a fluid flow. In the following sections, we first describe the Lattice Boltzmann model used in this study. We then investigate the effect of random roughness on transport and adsorption and finally study the influence of three different geometries on transport and adsorption.

I. METHODS

A. The Lattice Boltzmann method

The Lattice Boltzmann (LB) simulation method finds its origins in the 1980s and comes from bringing together the idea behind Lattice Gas Cellular Automata and concepts of statistical physics, through the Boltzmann equation. [16][17][18] As a lattice-based technique governed by local time-evolution equations, it is relatively simple to implement and to parallelize on multi-core systems. Moreover, local microscopic interactions can be readily implemented in the model, which is of high interest in our case for modeling fluid behavior in porous media. [START_REF] Succi | The Lattice-Boltzmann equation for fluid dynamics and bey[END_REF] Unlike classical computational fluid dynamic methods, the Lattice Boltzmann method does not solve explicitly the Navier-Stokes equation; however, by numerical integration of the Boltzmann equation it can be shown to satisfy the incompressible Navier-Stokes equation. [19] At the center of the Lattice Boltzmann method is the propagation of the one-particle velocity distribution function f (r, c, t) equivalent to the probability of a particle to be at node r of the underlying lattice, with velocity c at a given time t. Time, space and velocities are all discrete quantities in this scheme. Space is discretized by adopting a cubic mesh (or lattice) as a basis for the simulation. Velocities are discretized by projecting them on a finite number of lattice vectors. In three dimensions, several different models of discretized velocities exists, such as D3Q15, D3Q19, D3Q27 (featuring 15, 19 and 27 lattice vectors respectively). [START_REF] Succi | The Lattice-Boltzmann equation for fluid dynamics and bey[END_REF]18] Here we chose to use the D3Q19 model (see Fig. 1b), as a best compromise between precision and computational speed. [20] Time is discretized by integrating the propagation equation numerically, by finite time steps ∆t. The dynamics of the fluid on the lattice are governed by the following propagation equation:

f i (r + c i ∆t, t + ∆t) = f i (r, t) + (f e i (r, t) -f i (r, t)) τ + F ext i (1)
where f i is the component of f on velocity vector i, i.e. f i (r, t) = f (r, c i , t). The field f e i corresponds to the local Maxwell-Boltzmann equilibrium distribution, and τ is the relaxation time. The term F ext i accounts for external forces acting on the fluid and creating the fluid flow; in our case, they will correspond to a unidirectional pressure gradient throughout the system. This equation is implemented in our simulations following the method of Ladd and Verberg, [START_REF] Ladd | [END_REF] relevant for simulations of fluid dynamics in porous materials. We assume, in this type of materials, a laminar flow regime. The permeability of the fluid can thus be computed using the Darcy law:

K j Φ = νρ v j F j ext , (2)
where j = x, y, z corresponds to one of the three directions of space, v j , to the mean velocity of the fluid, F ext to the external forces, ν to the kinematic viscosity of the fluid, and ρ the volumetric mass density of the fluid.

B. The moment propagation method

To simulate the dynamical properties of solute dispersed in the fluid we use the moment propagation method proposed by Lowe and Frenkel [21,22] and further validated by Merks et al. [23] In this method, a propagated quantity P (r, t) is defined on the lattice which evolves following:

P (r, t + ∆t) = i P (r -c i ∆t, t)p i (r -c i ∆t, t) + P (r, t) 1 - i p i (r, t) (3)
where p i (r, t) corresponds to the probability of leaving node r with speed c i :

p i (r, t) = f i (r, t) ρ(r, t) -w i + w i λ with λ = 2D b v 2 T ∆t (4)
Here ρ is the fluid density, w i are constant weights of the speed model (D3Q19 in this case), D b is the diffusion coefficient of the tracers in the fluid in bulk phase, and v T is the fluid's speed of sound (v 2 T = 1 3 ∆x 2 /∆t 2 , with ∆x the lattice spacing). The propagated quantity is not a physically understandable parameter but, by its mathematical construction, provides access to the behavior of tracers inside the fluid. For a particular choice of the propagated quantity, namely the probability to arrive at position r at time t, weighted by the initial velocity of the tracers (in practice, one quantity is propagated for each component of the velocity), the velocity auto-correlation function Z is then computed by:

Z(t) = r P (r, t) i p i (r, t)c i (5)
The dispersion coefficient K is an interesting dynamical property of the tracers inside the fluid. [24][25][26] It quantifies the spreading of particles inside the fluid and may be defined from the standard deviation of the position of tracers over long times:

K = lim t→∞ σ 2 2t where σ 2 = r -r 2 (6)
with r the average position of the particles at the considered time. In practice, we compute it from the offsetted integration of the velocity auto-correlation function: [8,31]With the growing interest for the shale gas, we find also some studies about the transport and adsorption in kerogen pores. [32,33]More recently Long et al. [34] developped a new scheme to introduce all the IUPAC adsorption isotherms in Lattice Boltzmann scheme.

K = ∞ 0 [Z(t) -Z(∞)] dt (7)
We use in this work a novel Lattice Boltzmann model coupling transport of species and adsorption developed recently [14] and extended recently by the authors of the present paper to account for saturation and heterogeneity in the adsorbed density. [15] In this scheme adsorption takes place on the interfacial sites of the material it i.e. the fluid nodes having at least one neighboring solid node. The neighbors are detected following the D3Q19 speeds model. The adsorption process is described as an equilibrium between two populations: adsorbed and non adsorbed (free) species. The adsorption kinetic is set up using three parameters: the adsorption coefficient K a , the desorption coefficient K d and the saturation coefficient D max . The interplay between transport and adsorption is computed using the adsorbed and free densities:

D ads (r, t + ∆t) = 1 - D ads (r, t) D max D free (r, t)p a + D ads (r, t)(1 -p d ) D free (r, t + ∆t) =D free (r, t) 1 -p a + p a D ads (r, t) D max + D ads (r, t)p d (8)
where p a = k a ∆t/∆x and p d = k d ∆t. At t = 0. This two quantities are then equilibrated. This scheme corresponds to a Langmuir adsorption model. The adsorbed quantity n ads follows:

n ads (C ext ) = Q max m s κC ext 1 + κC ext (9)
where m s corresponds to the mass of material.

Q max = D max S s , C ext = C tot (1 -F a ) and κ = K a /(K d D max ).
The adsorbed fraction F a may be computed analytically with the relation:

F a = 1 + p d V p p a S s -1 (10)
Nevertheless this equation does not take account for saturation of the adsorption sites neither for eventual heterogeneities on the adsorbed density due to fluid flow. After this preliminary step to compute the adsorbed and free densities we propagate on the same scheme the two propagated quantities P and P ads to compute the dynamic of the tracers and their interactions with the adsorption sites to reach thermodynamic equilibrium:

P ads (r, t + ∆t) = 1 - D ads (r, t) D max P (r, t)p a + P ads (r, t)(1 -p d ) P (r, t + ∆t) =P (r, t) 1 -p a + p a D ads (r, t) D max + P ads (r, t)p d (11)
This last equation gives the framework to compute the Velocity auto-correlation function, the diffusion coefficient and the dispersion coefficient thanks to Eq. 5 and Eq 7.

D. Practical details

The simulation reported here are performed considering a laminar flow regime. We also use periodic boundary conditions on the three axis (x, y, z) and no slip boundary conditions at the liquid/solid interface for the Lattice Boltzmann scheme and the moment propagation method. We employed convergence criteria during simulations to ensure the convergence and verify the parameter rely on a the steady-state. We used a convergence criteria of 10 -14 (in relative step-to-step variation) for the average velocity of the fluid along the three directions of space, 10 -12 ∆x/∆t on the velocity autocorrelation function, 10 -11 for the step-to-step variations of the fraction adsorbed, and 10 -9 for the stepto-step variations of the dispersion coefficient. For the heterogeneity coefficient and its probability distribution function we did not use any convergence criteria but we performed several simulations at different number of time step to be sure to reach the steady state. The Lattice-Boltzmann scheme employed here works in reduced units. The results presented in this study are in scientific international (SI) units (except the mesh size x , y , z ). The method to switch between reduced units and SI units is available in supplementary information. Throughout the simulations we fixed: the bulk diffusion coefficient (D b = 6.04 10 -8 m 2 .s -1 ), the kinematic viscosity (ν = 10 -6 m 2 .s -1 ), the density of the fluid (ρ = 1000 kg.m -3 ) and the density of the solid (ρ s = 4970 kg.m -3 ).

II. IMPACT OF ROUGHNESS ON TRANSPORT AND ADSORPTION

From the published literature, many studies of adsorption and transport of fluids in porous media focus on simple geometrical model of the porosity, with "regular" or smooth surfaces. The impact of surface roughness at the local scale on the adsorption properties has been treated by some studies in gas phase [35,36] and on protein adsorption. [37,38] In an earlier Lattice Boltzmann study of adsorption and transport, we have seen an impact of local surface patterning, e.g. by comparing a smooth slit pore to one with grooves on the solid walls. [15] Herein we want to go further and investigate the effect of roughness on fluid flow and adsorption in a more realistic and geometrically complex model of pores with rough surfaces.

We focus here on roughness as a microscopic geometric heterogeneity on the internal surface of a pore of larger dimensions. The roughness thus represents a deviation -or the presence of defects -from an ideal geometry. It may have several origins, such as mechanical, chemical process or physical processes. It is omnipresent in real materials, but its scale and thus its impact depend drastically on the synthesis, activation, and chemical and physical history of each porous material. The effect of surface roughness has been studied in many research fields like biology [38,39], optic [40], coatings [41] or fluid dynamic. [42,43] For fluids in particular, in the case of hydrophobic interactions at the solid/liquid interface, the roughness can strongly affect the flow profile and in some cases it leads to a very low drop pressure due to slippage at the liquid/solid interface. [44,45] A. Generating rough surface models and measuring roughness

We describe here a simple model used to generate geometries of porous solids with roughness on their internal surface by a stochastic process of aggregation that mimicks the random deposition of nano-sized solid particles on the walls of an pre-existing pore system. To do so, we rely on the Lattice Boltzmann's underlying lattice vectors and definitions of neighboring nodes. Starting from an initial geometry (which we call skeleton), a fluid node is randomly selected. We evaluate its degree of connectivity with solid neighbors (see Fig. 1a) by computing an aggregation coefficient α:

α = (solid) w i (12)
α corresponds to a sum over all the node's solid neighbors, weighted by the coefficient of the D3Q19 speed model (see Fig. 1b). As an input parameter of the generation algorithm, we define the aggregation condition A c . If α > A c the node becomes solid, otherwise it remains fluid -this mimicks a process of aggregation of smaller particles, which are allowed to "stick" to the existing surface if the contact is large enough. Then, we repeat the aggregation process to another node chosen randomly again and again until we reach a convergence criterion on porosity or specific surface area.

After obtaining a new model of porous solid from this algorithm, we apply a filter to remove all nonconnected porosity (inaccessible cages) which may have been created during the aggregation process. This process ensures the connectivity of all the fluid nodes for the Lattice Boltzmann simulation and avoids artifacts in the moment propagation. [46] To do that we use a simple neighbor to neighbor propagation process. We initialize a quantity on one of the two sides of the simulation box orthogonal to F ext , Then we propagate the quantity from neighbor to neighbor. At the end of the propagation the fluid nodes where the quantity is not set up becomes solids. After that the process is repeated with an initialization on the opposite side of the simulation box.

This entire procedure allows us to create models of rough porous materials based on any given geometry defined on a cubic lattice, and through the parameter A c we can tune the extent of roughness. To quantify this, we define the following roughness coefficient:

R r = (|y i | -y i ) 2 (13)
corresponding to standard deviation of the minimal distance between the surface of the aggregated geometry and the original skeleton. There exist many other definitions of the roughness coefficient, mainly from the field of mechanical engineering, [START_REF] Degarmo | Materials and Processes in Manufacturing[END_REF] but as our goal here is merely to compare between different geometries a universal definition is not necessary. Figure 2 shows the evolution of the roughness coefficient as function of the number of aggregation steps. In this case, we chose as skeleton geometry a slit pore with a mesh size of x = 50∆x, y = 50∆x, z = 52∆x (size of the simulation box) and a convergence criterion on the porosity Φ = 0.7. Each point corresponds to the mean value of a set of 10 generations having the same input parameters. The errors bars corresponds to the standard deviation. We see that each aggregation condition A c gives rise to different values of roughness coefficient R r and its evolution as a function of number of steps in the generation algorithm. The error bars are small, showing that although we chose a stochastic procedure, the overall result is not very sensitive to the randomness. The value of R r stabilizes at high number of steps T max , and the parameter A c acts as a good control parameter to tune the roughness coefficient.

B. Impact of roughness and disorder

Before we set out to explore the influence of roughness on dynamical properties of transport and adsorption, in this section we quantify the impact of the randomness (or disorder) in the roughness on these properties. For this purpose, we created a series of 10 porous geometries for each value of A c (1/8, 1/16, 1/32) having otherwise the same values of inputs parameters [START_REF]Input parameters: ∆x = 0.5 µm, Ka = 6.04 m[END_REF] Looking at the standard deviation of the measured quantities related to adsorption and transport in the random pore spaces, we describe in Fig. 3 the roughness coefficient (R r ), the specific surface area (S s ), the permeability (K φ ), the average velocity of the fluid on y direction ( v y ), the fraction of tracers adsorbed (F a ), the mean pore size ( d p ), the porosity (Φ), the porous volume (V p ), the dispersion coefficient (K), and the Heterogeneity of adsorbed density ξ. Although there clearly are some variations on independent realizations of the rough geometry for a given value of roughness, this is relatively minor, with all the normalized standard deviations below 6%. The quantities most impacted (more than 2.5%) are related to heterogeneity of the adsorbed tracers at high roughness and the dispersion coefficient. Given the overall low sensitivity, for the two following sections, we will neglect the deviation caused by the random part of the roughness generation and describe the roughness of the surfaces simply by the R r coefficient.

C. Influence of adsorption on transport of tracers

We first study the impact of roughness on fluid properties, computing the velocity profile and the permeability coefficient in order to quantify the importance of roughness on geometries. Fig. 4 shows the flux profiles along the pressure drop for a non aggregated geometry (slit pore with an equivalent mean pore size of d p = 1.6 µm and plane perpendicular to z axis on top and down plane of the simulation box) and three aggregated geometries having A c = 1/8, 1/16, 1/32 on a slit pore (simulation box of size x = 50∆x, y = 50∆x, z = 52∆x, planes perpendicular to z axis on top and down of the simulation box and convergence criterion on porosity for roughness aggregation (Φ = 0.7).

The velocity of the fluid inside the pore decreases when the roughness is larger. The roughness obstructs the pore while the mean size of the pore d p (see Tab. I) remains constant. The difference in the flux profile comes from the deviations of the surface, namely the roughness. Computing the velocity profile for a slit pore with a pore size corresponding to the mean pore size of aggregated geometries highlights an unexpected effect. For low roughness value (A c = 1/8) the velocity is higher compared to the slit pore whereas we were expecting a value lower than the slit pore one. Investigating this unexpected effect, we showed that it appears as an artefact of the discretization in the simulation, and does not have physical meaning. In fact, close surface of the pore, the roughness creates some local pores only defined with one or two nodes -i.e., the length scale of the rugosity is close to the lattice spacing. Thus, the application of the bounce back rules (to ensure no-slip boundary conditions) combined with a single relaxation time may give a local dependence of the viscosity of the fluid on the local pore size. [START_REF] Ginzburg | [END_REF] To confirm this, we investigated the influence of mesh size on fluid behavior, by performing Table I. Values of the geometrical and permeability properties of the slit pore with and without roughness. simulations on the same geometry with several refined meshes. The results (detailed in supplementary information) show a dependence of the results of permeability coefficient on the discretisation of the little pores located in the roughness. The presence of this artefact in randomly generated surface of pore suggests to be careful when we generate geometries, especially those coming from tomography pictures where the surface roughness is poorly controlled. Fig. 4b shows the evolution of the permeability coefficient as a function of the roughness, and this evolution is the same as that of the fluid's velocity: it decreases with increasing roughness coefficient. In the regime studied, the evolution appears rather linear. The values of the permeability coefficient confirm the unexpected behavior seen previously on the flux profiles, namely that the roughness with A c = 1/8 gives a higher permeability than the slit pore. This counter-intuitive effect is very interesting because it offers the opportunity to improve the materials. Having the same mean pore size it is possible to decrease the drop pressure just with the introduction of some controlled roughness. Plus, the introduction of this roughness will increase the specific surface area, i.e. the adsorption capacity. We demonstrate here a way to increase the adsorption capacities and the permeability at the same time, through geometrical tuning of the inner pore surface.

We now turn to the dispersion coefficient, which is representative of the spreading of tracers in the fluid. In previous work we have studied the influence of an ordered roughness (a slit pore with crenelated pores on the walls) and have shown the ratio between the crenels' height and width r = h/w, have an influence on the dispersion coefficient in presence of adsorption: the dispersion coefficient increase with the ratio r. [15] As the crenels are analogous to an ordered roughness of the pore surface, we expect to observe something similar here.

And indeed, the random roughness generated here has also an influence on the dispersion coefficient, as is shown on Fig. 5a. For a given roughness value, the influence of F a on the dispersion coefficient is the same as the one observed previously with slit pore and crenelated pores. [15,50] For low adsorption strength, all the tracers are free and have a dispersion due to fluid flow. When the tracers start to adsorb, they create a disparity of position compared to the tracers free in the fluid: the dispersion coefficient increases. For high adsorption strength the majority of the tracers are immobile (adsorbed): the dispersion tends to zero. At intermediate regime we thus observe a maximum of dispersion. [15] Nevertheless, the influence of roughness on dispersion is inverted. With crenelated pores the dispersion coefficient increases when the roughness increases, whereas in the present case the dispersion coefficient surprisingly decreases when the roughness coefficient increases.

Fig. 5b shows the dispersion coefficient as a function of the mean velocity of the fluid. The curves obtained for various values of roughness cross together in a single point. That means the order between the curves of Fig. 5a may change regarding to the mean velocity of the fluid. This was clearly not the case when we studied crenelated pores. Another parameter play a key role in the phenomena involved here at an order of magnitude higher than the roughness coefficient: the mean pore size. When we studied crenelated pores in earlier work, we kept constant the distance between the tops of the crenels. Here, in contrast, the mean pore size is kept constant. Fig. 5c proves that the pore opening is the key parameter. Setting the mean pore size constant, it shows the evolution of the dispersion coefficient as a function of the mean velocity of the fluid in different slit pore crenelated geometries having different value of r. We obtain the same behavior as for random roughness.

This means, in terms of materials design, to increase the separation performance of the materials by introducing some roughness, the minimal pore size (or minimal opening diameter) is to be considered as key parameter, rather that the average pore size.

D. Influence of flow on adsorption

In previous work we have shown the fluid flow can create some heterogeneity in the adsorbed density, taking away species from the upstream adsorption sites and accumulating tracers in downstream sites. [15] Here we wanted to see if the same phenomenon occurs at the local scale when disordered roughness is present (in contrast to our previous work on regular grooved pores). Figure 6 plots the relative deviation Φ ads of the adsorbed density:

Φ ads = D ads -D ads D ads ( 14
)
in presence of flow (F ext = 5 10 9 Pa/m) for three different values of the roughness coefficient in a slit pore. We notice on Fig. 6 some disparities in the adsorbed density toward the flow. These heterogeneities are limited locally to a few percents, and we report in table II the value of the heterogeneity coefficient ξ in each case: Unlike for the crenelated slit pore geometry ξ does not increase monotonically with the aggregation condition. The value at A c = 1/32 is lower than the value at A c = 1/16. We already know from previous work the heterogeneity is strongly dependent on the velocity of the flow and here the rise up of the roughness coefficient makes the fluid velocity decrease. We have a competition between the influence of the size of local geometrical cavities (roughness) and the speed of the fluid. This phenomenon, established in prior work on model geometries, is here shown to be generic and applicable to disordered and rough pore surfaces.

ξ = φ 2 ads ( 15
)

III. GEOMETRY COMPARISON

Comparing porous materials for adsorption applications under fluid flow, two main parameters need to be taken into account to judge their efficiency: the total adsorbed quantity (or adsorption capacity), and the permeability (to ensure the lowest pressure drop). We aim here at finding a way to compare materials with the same "chemistry", i.e. locally the same adsorption sites, but with distinct pore geometries, and study the influence of the geometry on the two characteristics of adsorption and transport.

To do so, we choose three totally different geometries, displayed on Figure 7: a honeycomb geometry having Figure 7. Lattice-based models for three different pore geometries, and electron microscopy images from real-life materials with similar geometries. Left: honeycomb geometry; [51] center: inverse opal; [52] right: spinodal decomposition geometry. [53] straight and smooth pores, a sphere replica geometry having spherical interconnected pores, and a very disordered geometry with worm like pores. The honeycomb geometry is typical of materials synthesized using icetemplating methods, [54] the sphere replica is the characteristic pore space of inverse opal materials, [55] and the worm-like porosity is archetypal of materials produced by spinodal decomposition. [56,57] We first describe the procedures we followed to created lattice-based models of these geometries, and then go on to discuss their relative performance for adsorption and fluid transport.

A. Geometries generation on lattice

Honeycomb -The honeycomb geometry is simple due to its translational invariance. To create a latticebased model we use the same technique as the one used usually in Computational Aided Design. In a plane we design the cross section of the geometry: assembled hexagons. To fit with the nodes on the grid the hexagons does not have the same edges length (see Fig. 8). Assuming a is the length of the horizontal edges the basic mesh size is 2(2a -1) horizontally and 2(a -1) vertically. Once the cross section is created we extrude the profile along the direction perpendicular to the plane to have the 3D geometry. To tune the amount of porosity of this geometry we assign some thickness w to the hexagonal profile.

Inverse opal -The inverse opal geometry is also highly symmetric and can be constructed by similarly (see Fig. 8). Starting with a solid block we create a spherical holes of diameter d s . Then we replicate this hole in each direction at a distance d int . This leaves windows of diameter d r between the spherical cavities:

d r = d 2 s -d 2 int ( 16
)
Spinodal decomposition -Materials produced by spinodal decomposition -a phase separation process -feature pore geometries that are very disordered and worm-like channels with no symmetry. There are studies in the existing literature proposing numerical models of spinodal decomposition materials, [58][59][60] but for our purposes taking into account the whole thermodynamic of such process would be to much time consuming. To create lattice-based models of very disordered geometries, like spinodal decomposition materials, we propose here a simple method based on the Ostwald ripening principle. [61] This principle concerns the ability of phases to rearrange themselves to minimize surface energy, with small droplets tending to regroup themselves to form bigger ones, a process which is easy to model on a lattice. Thus, we start by initializing the system as a random distribution of solid/fluid nodes on the cubic lattice, with a fixed ratio (the initial porosity). Then we choose a node randomly, computes the sum of its links with neighbors having the same nature (i.e. the sum of fluid-fluid or solid-solid neighbors), accounting for periodic boundary conditions. The sum of the links, which we note σ, is weighted with the weights of the D3Q19 velocity scheme (see Fig. 1b). We then compare σ to a threshold value, which we set at σ c = 0.4. If σ > σ c , the node stays the same otherwise we switch its nature (liquid node becomes solid, and solid node becomes liquid). This step is repeated a number of N times. After the generation, we use a filter to remove unconnected solid regions and inaccessible porous regions. Figure 9 compare the pore size distribution obtained with our algorithm and some pore size distribution of spinodal decomposition available in the literature. [62][63] The data is computed from a sampling of 10 geometries having the same input parameters. The errors bars represents the standard deviation of the data obtained. The shape of the experimental and the one we have computed are very close. The algorithm we have developed is representative of the geometries obtained experimentally using spinodal decomposition. Moreover, the algorithm is very fast. Its takes less than one minute to create a geometry on a 100 × 100 × 100 mesh. Finally, both the overall porosity (fraction of porous volume) and the amount of tortuosity can be controlled by the two parameters of the algorithm, namely the initial porosity and the length N of the ripening process.

B. Comparing pore geometries

Physical properties such as fluid transport and adsorption in porous materials are strongly linked to the geometry of the material's pore system. Modifying the pore geometry, modifies on one hand the porosity and the specific surface area and on the other hand the behavior of fluid and the motion of species. Transport and adsorption are often inversely linked. Modifying a material to improve adsorption skills generally decrease its transport properties and improve transport properties generally decrease adsorption skills. This makes the comparison between materials tricky. To avoid an enormous generation of data and objectively compare our three geometries in terms of transport and adsorption we can either maintain constant the transport properties and see the influence on adsorption or maintain the adsorption properties and see the influence on transport properties. The easier way,in our case, is to keep constant the adsorption and see the influence on the transport properties. To do so we have to keep constant the adsorbed quantity per mass of material (Q a ). We already know for low concentration of tracers that

F a = 1 + K a S s K d V p -1 (17)
where F a is the fraction of tracers adsorbed (the ratio between the amount of tracers adsorbed and the total amount of tracers), K a and K d are respectively the adsorption and desorption coefficient, S s corresponds to the specific surface area and V p corresponds to the porous volume. Q a can be written as a function of F a :

Q a = C i V p ρ s V s F a (18)
where C i represents the initial concentration of tracers, ρ s is the volumetric mass of the solid part of the material and V s is the volume of solid. Considering Equations 17 and 18:

Q a = C i V p ρ s V s 1 + K a S s K d V p -1 (19)
For our simulations we consider the adsorption sites of the three geometries have the same characteristics (K a = 6.04 m.s -1 and K d = 6.04 10 6 s -1 ), i.e. we study only the influence of geometry at a fixed chemical composition of the porous material's walls. We also consider the solid part of the material as the same nature (ρ s = 4970 kg.m -3 ) and the initial concentration of tracers is constant (C i = 1 g.L -1 ). Considering this, to keep Q a constant we have to keep the ratios V p /V s and S s /V p constant. This is equivalent to keep the porosity Φ and the ratio S s /V p constant.

In order to do so, we have tuned the geometry to have the same porosity Φ and then adjust the ∆x (distance between two nodes) to have the same ratio S s /V p . We adjust the thickness of the wall of the honeycomb to make variation of the porosity. We adjust the distance between the spheres d int to vary the porosity in inverse opal geometry and we vary the number of time step to make variations on the spinodal decomposition geometry. A table with the detailed characteristics of the geometries is available in supplementary information. The adsorption isotherms of the three geometries are also available in supplementary information and show the adsorbed quantity is the same for all the geometry at low and high concentration in solute.

C. Pore geometry influence on adsorption

Having shown earlier that the fluid flow may create some local heterogeneities in the adsorbed density, we want to analyze it here and investigate whether the geometry has any influence on this heterogeneity. Figure 10 presents a 2D cut view of Φ ads for each material. The colored gradient represents the values of the adsorbed density. The flow creates disparities in case of the spinodal decomposition geometry, between the upstream and the downstream part of the internal surface. This effect happens also in the inverse opal geometry with a more visible deviation between the two sides of the geometry. No disparities occur in the honeycomb thanks to its slick surface oriented along the flow. The heterogeneity is one order of magnitude higher in the spinodal decomposition than in the inverse opal.

To compare quantitatively the heterogeneities, Figure 11 shows the probability distribution of Φ ads for each geometry. The distributions are completely different. The distribution is a single peak for the honeycomb, showing absence of any heterogeneity in the adsorbed density. Indeed, the smooth walls oriented perfectly along the flow do not create heterogeneities. The inverse opal distribution has a bimodal shape coming from the two populations of the adsorbate stuck on the upstream and downstream part of the geometry and the spinodal decomposition's one is a Gaussian like function. In this geometry the heterogeneity is "averaged" by the randomness of the geometry.

As a conclusion we see here that both the topology of It's a counter-intuitive effect of keeping the adsorption constant. This means the variations of porosity only influence the solid part of the material (the walls). Increasing the porosity just increases the thickness of the wall and does not modify the void part of the material. This would mean the ratio S s /V p is constant because S s and V p are constant independently but it is not the case here. The phenomenon is more complex than just increasing the thickness of the walls.

The pore size of the honeycomb is almost twice bigger than the others. Straight pores allows to have the largest pore size. As a consequence the value of the honeycomb's permeability (K Φ ) is 2.5 times higher than the spinodal decomposition and 5 times higher than the inverse opal (see Fig. 12b). This demonstrates the real interest of having porous materials with straight pores. At equivalent adsorption skills it gives the highest permeability.

Figure 12c presents the values of the parameter:

ψ = K Φ d 2 p ( 20
)
ratio between the permeability and the squared pore size. Note this ratio is dimensionless. ψ is constant for all the geometry. For this study the permeability only depends on the pore size and the pore geometry has no influence on it. This means to have a material with the best permeability skills for a given adsorption capacity we have to find a way to increase the pore size independently to the porosity and the ratio S s /V p .

IV. CONCLUSION

We have studied here the interplay between adsorption and transport in porous materials under liquid flow, and the impact of the geometry of the pore system on these two properties. By using a Lattice Boltzmann scheme extended to take into account adsorption of tracers in the liquid phase, we showed how adsorption and fluid transport are both affected by global geometric characteristics (pore shape and alignment with the fluid flow) as well as local geometric features (such as roughness of the pore surface). In particular, we showed that roughness of the pore walls effectively modifies the nature of the solid/fluid interface, introducing slippage in a system which would otherwise have a no-slip boundary condition. Moreover, we generated realistic models of complex experimental materials and quantified the impact of geometry on fluid transport and tracer adsorption. This sheds light into the optimization of materials for applications in dynamic separation of species by adsorption under fluid flow. Future work will address the kinetics of fluid adsorption and the dynamics at the scale of the adsorbent sample, to bring the Lattice Boltzmann technique closer to model flow experiments in, e.g., liquid chromatography. Moreover, more work will be necessary to replace the complex model geometries used in this work -as realistic as they may be -with actual 3D images of real-life materials, obtained for example by X-ray tomography.

Figure 1

 1 Figure 1. a. Roughness generation process based on LB weights. b. D3Q19 speed model scheme.

Figure 2 .

 2 Figure 2. Mean value of roughness coefficient as a function of the number of steps of the generation procedure, with different aggregation conditions Ac. Values come from a sampling of 10 geometries having the same input parameters, and error bars correspond to the standard deviation.

Figure 3 .

 3 Figure 3. Sensitivity of adsorption and transport properties on roughness random aggregation. Rr: roughness, Ss: specific surface area, K φ : permeability coefficient, v : average velocity of the fluid, Fa: fraction of tracers adsorbed, dp : mean pore size, Φ: porosity, Vp: porous volume, K: dispersion coefficient, ξ: spatial heterogeneity of adsorbed density.

6 KΦ (m 2 )

 62 2.1 10 -13 2.3 10 -13 1.8 10 -13 1.4 10 -13

Figure 4 .

 4 Figure 4. a: Velocity profile in a slit pore with roughness (Ac = 1/8, 1/16, 1/32) compared with a flat slit pore of equivalent mean pore size dp = 1.6 µm. L corresponds to the distance between the planes of the slit pore. Here Fext is kept constant to 5 10 9 Pa/m. b. Effect of roughness on permeability for Ac = 1/8, 1/16, 1/32.

Figure 5 .

 5 Figure 5. a. Influence of roughness on dispersion coefficient (K) as a function of the fraction adsorbed (Fa) at constant value of Fext = 5 10 9 Pa/m, Ka = 6.0 m.s -1 and K d = [6.0 10 8 ; 3.4 10 8 ; 1.9 10 8 ; 1.1 10 8 ; 6.3 10 7 : 6.0 10 7 ; 2.4 10 7 ; 9.6 10 6 ; 3.8 10 6 ; 1.5 10 6 ] s -1 . b. Influence of roughness on dispersion coefficient (K) as a function of the average velocity of the fluid ( vy ) with Ka = 6.0 m.s -1 , K d = 6.0 10 7 s -1 and Fext = [2 10 8 ; 4 10 8 ; 6 10 8 ; 8 10 8 ; 10 10 8 ] Pa/m. c. Influence of crenelated pores on dispersion coefficient (K) as a function of the average velocity of the fluid ( vy ) with Ka = 6.0 m.s -1 , K d = 6.0 10 6 s -1 and Fext = [2 10 8 ; 4 10 8 ; 6 10 8 ; 8 10 8 ; 10 10 8 ] Pa/m.

Figure 6 .

 6 Figure 6. Relative deviation in the adsorbed density (Φ ads ) in presence of flow with different roughness. a. Adsorbed density, Ac = 1/8. b. Adsorbed density, Ac = 1/16. c. Adsorbed density, Ac = 1/32. Points that appear unconnected to the porosity, in this 2D slice, are actually connected along the direction perpendicular to the plane of the figure.

Figure 8 .

 8 Figure 8. Sketchup of honeycomb and inverse opal geometries and their geometric parameters.

Figure 9 .

 9 Figure 9. Left: pore size distribution of a lattice-based model created through the Oswald ripening procedure. Right: experimental pore size distribution of materials synthesized by spinodal decomposition.[62]

Figure 10 .

 10 Figure 10. Relative adsorbed density (Φ ads ). Grey: solid part, blue: fluid part, color gradient: adsorption sites. a: Spinodal decomposition cut view at y = 1. b: Inverse opal cut view at y = 45. c: Honeycomb cut view at z = 1.

Figure 11 .

 11 Figure 11. Probability distribution function of the relative adsorbed density (Φ ads ), with porosity Φ = 70%.

Figure 12 .

 12 Figure 12. (a) Mean pore size (dp) of the three different geometries at constant adsorbed quantity (Qa). (b) Geometry impact on permeability coefficient (K) at constant adsorbed quantity, Qa. (c) Geometric coefficient (ψ) for the three different geometries at constant adsorbed quantity. The comparison is represented for three different values of porosity Φ (as percentage).

Table II .

 II Heterogeneity of the density adsorbed for aggregation conditions 1/8, 1/16, 1/32. Ka = 6.04 m.s -1 , K d = 6.04 10 6 s -1 .

	Ac	1/8	1/16	1/32
	ξ(×10 -3 ) 4.05 ± 0.04 4.4 ± 0.08 4.11 ± 0.12
	< Vy > (m/s)	1.16	0.892	0.702

ACKNOWLEDGEMENTS

This work was supported by the ANRT through CIFRE sponsorship 1262/2013. We gratefully thank Benoît Coasne for fruitful discussions, as well as Benjamin Rotenberg and Maximilien Levesque for providing us with their Lattice Boltzmann code and many scientific discussions.


Jean-Mathieu Vanson

François-Xavier Coudert
email: fx.coudert@chimie-paristech.fr

Michaela Klotz

Anne Boutin
email: anne.boutin@ens.fr

Kinetic accessibility of porous materials adsorption sites studied through Lattice Boltzmann method

We present here a computational model based on the Lattice Boltzmann scheme to investigate the accessibility of active adsorption sites in hierarchical porous materials to adsorbates in a flowing liquid. By studying the transport and adsorption of tracers after they enter in the pore space of the virtual sample, we characterize their kinetics as they pass through the pore space and adsorb on the solid-liquid interface. The model is validated on simple geometries with known analytical solution. We then use it to investigate the influence of regular grooves or disordered roughness on the walls of a slit pore geometry, looking at the impact on adsorption and transport. In particular, we highlight the importance of adsorption site accessibility, which depends on shape and connectivity of the pore space as well as the fluid flow profile and velocity.

Introduction

Materials with hierarchical porosity, ranging from nanometers to micrometers, are of interest for a wide range of industrial applications, because they can potentially combine high specific surface area and high permeability, i.e. low resistance to fluid flow. Such applications typically rely on the host-guest surface interactions, for example in phase separation, gas mixture separation, or ions exchange and capture. In the liquid phase, practical applications at large scale include, for example, water decontamination and removal of pollutants such as heavy metals or radioactive ions.

Understanding the physical and chemical phenomena involved in fluid transport and adsorption in hierarchical porous materials is crucial for the optimization of existing materials and the design of new materials. In particular, the optimization of the pore space (pore dimensions and pore geometry) for a chosen application is a great challenge. The study of transport and adsorption in hierarchical porous materials raises several open fundamental questions. One of them deals with the kinetics of exchange between interconnected pores of different sizes -between macropores and mesopores, or between mesopores and micropores. The influence of the macroporosity on the accessibility of species to the mesoporosity is both a very fundamental question and yet one directly relevant to applications. It depends in nontrivial ways on the geometry of the macropores and mesopores, as well as the topology and interconnectivity of the two pore networks, each of these factors affecting the overall performance of the material.

Despite its fundamental importance and impact, the pub-lished literature is relatively scarce on investigations of kinetic accessibility of adsorption sites in porous materials with hierarchical pore networks. The problem of fluid flow in complex geometries and chromatographic systems has been widely addressed by computational studies, including many based on the Lattice Boltzmann method. We describe here another approach, building upon a Lattice Boltzmann model recently extended to incorporate an accurate description of adsorption of solutes at the fluidsolid interface. [START_REF] Levesque | Accounting for adsorption and desorption in lattice Boltzmann simulations[END_REF][START_REF] Vanson | Unexpected coupling between flow and adsorption in porous media[END_REF] However, this model cannot in itself answer the question of the accessibility, as it describes the steady state of fluid flow and solute dynamics. In this work, we present an extension of the Lattice Boltzmann model to evaluate the accessibility of active sites by simulating a system out of the steady state. We illustrate it on the kinetics of adsorption on a pore model of a crenelated slit pore, i.e. a larger main pore interconnected with smaller side pockets.

Methods

The method used in this work is based on a well-known numerical simulation scheme for fluid dynamics, namely the Lattice Boltzmann (LB) model. [START_REF] Succi | The Lattice-Boltzmann equation for fluid dynamics and beyond[END_REF] The Lattice Boltzmann model finds its roots in the 1980s and has several advantages compared to other computational fluid dynamic methods, making it suitable to address systems at scales intermediate between atomistic or coarse grained simulations (at smaller scale) and finite element-based computational fluid dynamics (at larger scale), as shown on Figure 1. In particular, the relative simplicity of the computational method make it easy to implement and to parallelize. The LB method works through the propagation of the oneparticle velocity distribution function f (r, c, t) equivalent to the probability of a particle to be at position r with velocity c at time t. For numerical integration, time, space and velocities are all discretized. The distribution function is propagated in time through the following equation:

fi(r + c i ∆t, t + ∆t) = fi(r, t) + (f e i (r, t) -fi(r, t)) τ + F ext i ( 1
)
where fi is the component of f on velocity vector i, i.e. fi(r, t) = f (r, ci, t). f e i corresponds to the local Maxwell-Boltzmann equilibrium distribution, τ is the relaxation time, and F ext i accounts for external forces acting on the liquid and responsible for the fluid flow. This equation is implemented following the method of Ladd and Verberg, 8 relevant for simulations of fluid dynamics in porous materials. To simulate the dynamical properties of solute dispersed in the fluid we use the moment propagation method proposed by Lowe and Frenkel 9,10 and further validated by Merks et al. [START_REF] Merks | The Moment Propagation Method for Advection-Diffusion in the Lattice Boltzmann Method: Validation and P'eclet Number Limits[END_REF] A quantity P (r, t) is defined on the lattice and propagated following:

P (r, t + ∆t) = i P (r -c i ∆t, t)pi(r -c i ∆t, t) + P (r, t) 1 - i pi(r, t) (2)
where pi(r, t) corresponds to the probability of leaving node r with velocity ci:

pi(r, t) = fi(r, t) ρ(r, t) -wi + wiλ with λ = 2D b v 2 T ∆t (3)
Here ρ is the fluid density, wi are constant weights of the velocity model (we chose in this work the D3Q19 model for discretization of velocities), D b is the diffusion coefficient of the tracers in the fluid in bulk phase, and vT is the fluid's speed of sound (v 2 T = 1 3 ∆x 2 /∆t 2 , with ∆x the lattice spacing).

In the published literature, there are relatively few computational lattice-based fluid dynamics models that explicitly take into account the dispersion of solutes and their adsorption at fluid-solid interfaces. [START_REF] Agarwal | A lattice boltzmann model for adsorption breakthrough[END_REF][START_REF] Manjhi | Lattice Boltzmann modelling of unsteady-state 2D concentration profiles in adsorption bed[END_REF][START_REF] Zalzale | Lattice Boltzmann simulations of the permeability and capillary adsorption of cement model microstructures[END_REF][START_REF] Anderl | Numerical simulation of adsorption and bubble interaction in protein foams using a lattice Boltzmann method[END_REF][START_REF] Pham | Transport of nanoparticles and kinetics in packed beds: A numerical approach with lattice Boltzmann simulations and particle tracking[END_REF][START_REF] Hlushkou | How Microscopic Characteristics of the Adsorption Kinetics Impact Macroscale Transport in Chromatographic Beds[END_REF][START_REF] Ren | A Lattice Boltzmann Model for Simulating Gas Flow in Kerogen Pores[END_REF][START_REF] Ning | Numerical modeling of slippage and adsorption effects on gas transport in shale formations using the lattice Boltzmann method[END_REF][START_REF] Long | Modeling adsorption with lattice Boltzmann equation[END_REF][START_REF] Yang | Quantitative determination of molecular propagator distributions for solute transport in homogeneous and heterogeneous porous media using lattice Boltzmann simulations[END_REF][START_REF] Gray | Simulating dispersion in porous media and the influence of segmentation on stagnancy in carbonates[END_REF][START_REF] Gray | Simulation of dissolution in porous media in three dimensions with lattice Boltzmann, finite-volume, and surface-rescaling methods[END_REF] In this work, we account for adsorption at the solid-liquid interface by following the technique developped by Levesque et al. [START_REF] Levesque | Accounting for adsorption and desorption in lattice Boltzmann simulations[END_REF] and refined in later work from our group, 6,24 -accounting for adsorption and desorption in a generic lattice Boltzmann scheme, by modelling these processes through kinetic rates of adsorption and desorption taking place at interfacial lattice nodes. The adsorption of tracers in the fluid occurs on fluid nodes having at least one solid node as a neighbour. It is described using the adsorbed density of tracers D ads , which is defined only on adsorption sites, and the free density of tracers D free , defined everywhere in the fluid. The adsorption kinetic is described by three physical parameters: the adsorption coefficient Ka, the desorption coefficient K d , and the saturation coefficient Dmax (or maximal uptake of a fluid node). The balance between adsorbed and free density is integrated numerically in time following first-order kinetics:

6 D ads (r, t + ∆t) = 1 - D ads (r, t) Dmax D free (r, t)pa + D ads (r, t)(1 -p d ) D free (r, t + ∆t) =D free (r, t) 1 -pa + pa D ads (r, t) Dmax + D ads (r, t)p d ( 4
)
where pa = Ka∆t/∆x and p d = K d ∆t.

In the simulations reported in this work, a slab of porous material (our "sample") is placed within a bigger simulation box, whose dimension is larger along the direction of fluid flow (see Figure 2). This leads to the presence of two fluid reservoirs, one upstream from the porous matrix, and the other one downstream. We consider a no slip boundary condition at the fluid/solid interface, and periodic boundary conditions at four faces of the simulation box (in the xy and yz planes). However, the periodic boundary conditions are removed in the moment propagation on the two faces perpendicular to the flow (in the xz plane). This general approach is relatively simple to implement, and has been used and validated once in earlier literature for lattice-based computational fluid dynamics of flow, dispersion and catalysis in a packed bed of catalytic spherical particles. [START_REF] Zeiser | CFD-calculation of flow, dispersion and reaction in a catalyst filled tube by the lattice Boltzmann method[END_REF] In order to compute the transient kinetics of the tracers in the fluid we proceeded to modify the initialization of the free density at time t = 0. The free density is homogeneously initialized on the upstream face of the simulation box perpendicular to the flow, i.e. the free density is initialized to a constant value on all fluid nodes in the plane y = 0 when the flux is applied along the y direction, while the initial free density on all other fluid nodes is set to zero.

Practical details

The data shown here are obtained in a laminar flow regime. We use periodic boundary conditions on two axes (x and z), and no slip boundary conditions at the liquid/solid interface for the Lattice Boltzmann scheme. We used a convergence criterion of 10 -14 (in relative step-to-step variation) for the average velocity of the fluid along the three directions of space. The Lattice Boltzmann scheme employed here works in reduced units (∆x, ∆t, ∆p). The method to switch between reduced units and SI units is available in We used the open-source Lattice Boltzmann code laboetie, 5 available online at https://github.com/ maxlevesque/laboetie. Simulations reported in this work were performed with a modified version of the code, available at https://github.com/maxlevesque/laboetie/ tree/jmvanson (commit 809d0b7). Input files and additional information can be found at our group's data repository, at https://github.com/fxcoudert/citable-data.

Model validation

This section is devoted to the validation of the modified moment propagation scheme, without periodic boundary condition in the flow direction. The key quantity for this is the dispersion coefficient K, which represents the spreading of the species inside the material, accounting for both advection and diffusion. It characterizes the influence of the material on the motion of species. In a periodic system it can be computed from the velocity auto-correlation function Z as:

K = ∞ 0 [Z(t) -Z(∞)] dt (5)
However, in the specific case of our modified technique with no periodic boundary conditions in one direction, we validate the model by computing the dispersion coefficient using a macroscopic procedure common in the field of chromatography. [START_REF] Wernert | Influence of Molecule Size on Its Transport Properties through a Porous Medium[END_REF] At a given position y = lm in the column (see Fig. 2), we measure the flux of tracers passing through the section as a function of time:

γ(t) = lx i=1 lz k=1 D free (i, lm, k, t) (6)
The evolution of γ with time is the familiar elution curve, an example of which is represented on Figure 3 for a pore with slit geometry, a height of l = 50∆x, where lm has been set to 800δx, with no adsorption (Ka = 0) and Fext = 2E -6∆p/∆x.

Once the simulation has finished and the elution curve is obtained, we fit it to a mathematical model. An impressive number of functions have been proposed in the literature to fit elution curves, and a review on the topic counts up to 90 different mathematical functions used for this purpose. [START_REF] Marco | Mathematical functions for the representation of chromatographic peaks[END_REF] Here, we have chosen simplicity and fit elution curves with a Gaussian function:

f (t) = Ae -(t-t R ) 2 2σ 2 t ( 7
)
The Gaussian is the most widely used model in the literature, and our goal here is not to test in detail the influence of one model or another, but to validate our moment propagation scheme. The results obtained do not depend much on the shape of the fitting function, as we only use the peak position and half-maximum width. From the Gaussian fit, we then calculate the Height Equivalent to a Theoretical Plate (HETP), H, with the following equation: 28

H = lmσ 2 t t 2 R ( 8
)
where σt is the width of the Gaussian fit and tR is the elution time (time t at which γ(t) is maximal). Finally, the dispersion coefficient K is computed using the HETP as: 26,29

K = vy H 2 = vy lmσ 2 t 2t 2 R ( 9
)
where vy corresponds to the average velocity of the tracers in the fluid. We can see in the simplest test case, that of a slit pore with no adsorption (Ka = 0) (Fig. 4), that the value of the dispersion coefficient measured through the above procedure is clearly dependent on the position of the measurement (the value of lm) and converges at high lm. This converged value, obtained far from the sample, increases with external force Fext as is expected. We then compare the values of the dispersion coefficients thus obtained with the analytical solution, which is known for a pore with slit geometry (of width l h ): 5,6

K D b = 1 -Fa + Pe 2 102α 2 + 18α + 1 210(1 + 2α) 3 + D b l 2 h K d 2α (1 + 2α) 3
(10) where Pe = l h v /D b is the Péclet number and α = Ka/K d h.

Figure 5 compares the values obtained from our model with the analytical solution, for cases without adsorption and in presence of adsorption (panels a and b respectively). We can see that the results agree reasonably well with analytical values, especially considering the numerous assumptions made in the indirect "macroscopic" determination (in particular the Gaussian fit). The difference between indirect and analytical results increases with the velocity of the fluid due to the fact that the convergence of K with lm is slower at high value of the fluid velocity. This difference also increases with the ratio Ka/K d because we employed a simplified Van Deemter equation (Eq. 9).

Results and discussion

In this section, we use the extended Lattice Boltzmann model to characterize the evolution with time of the tracers' density inside the fluid and adsorbed at the fluid-solid interface. This yields useful information about the tracer kinetic and the accessibility of the active adsorption sites in two different geometries.

Crenelated pore

We first study the kinetics of solute transport and adsorption in a model geometry of a slit pore with crenels. The geometry is characterized by the aspect ratio of the crenels (or grooves), i.e. the ratio r between the crenel depth and crenel width. Figure 6 shows the time evolution of the density of tracers in the fluid in a slit pore with grooves having an aspect ratio of r = 1 (whose depth and width are equal). Three snapshots of the density are taken at three moments along the simulation, as the "pulse" of tracers, starting from the left side of the picture, moves towards the right side following the fluid flow. We can see that in this case, the tracers' density follows the profile of the fluid flow, with faster flow in the center of the pore. We can also observe the effect of adsorption at the solid-liquid interface on Fig. 6b, where it seems clear that the solute has access to the whole internal area of the pores (i.e., that all adsorption sites are active within the timeframe of the simulation).

We depict the time evolution of tracers in a slit pore with deeper crenels (r = 4, all other parameters being equal) on Figure 7. In contrast with the previous case (aspect ratio r = 1), it is clear that for the deeper crenels the tracers do not have access to the entirety of the pores' internal surface. With the chosen flow, it is clear that advection along the y axis is too fast compared to the diffusion of the tracers along the perpendicular x and z axes, which would give them access to the bottom of the grooves. Thus, a dead volume exits in the material, and not all adsorption sites are active in these given conditions. This result parallels what we have seen in earlier work, [START_REF] Vanson | Unexpected coupling between flow and adsorption in porous media[END_REF] in which we studied the velocity profile for fluid flow in crenelated pore systems and showed the existence of "dead volume", i.e. a non negligible part of the crenel where the velocity of the fluid is close to zero. We observe that this dead volume, which can depend on the velocity of the fluid, becomes significant for r > 1, when the grooves' depth is larger than their width, as the fluid flow cannot enter the channel.

To quantify this effect we have plotted in Figure 8 the evolution of the adsorbed fraction as a function of simulation time, for various values of crenels' aspect ratio r. In all cases, we see that the adsorbed fraction (i.e., the adsorbate uptake) goes through a maximum as a function of time, then decays at long times. The maximum corresponds to the moment when the "cloud" of tracers goes through the sample. For r < 1 the maximum of the adsorbed fraction increase clearly with r: there, all the adsorption sites are accessible to the adsorbate, and the adsorbed fraction is proportional to the specific surface area, which itself grows as a function of r (since we keep the distance between the top of the crennels constant). For r > 1, on the other hand, the adsorbed fraction is lower than for r < 1. Even though, the total specific surface area increases, the adsorbed fraction decreases because not all the sites are accessible for a given fluid velocity. This highlights the fact that, for an adsorbent material to perform best under fluid flow, it is not merely dictated by the absolute specific surface area, but also by its accessibility in specific conditions, as the system is not at thermodynamic equilibrium (as is often considered the case for gas adsorption) or in a steady state.

Random roughness

The crenels studied in the previous section can be considered as an ordered "roughness" or small-scale disorder on the surface of the walls of a slit pore. Here, we extend this study by looking at the impact of a disordered (or "random") roughness of the adsorbent walls on the kinetics of tracers in the fluid. The process by which our rough slit pore models were generated -which relies on aggregation of smaller particles onto pre-existing walls -is described in Supporting Information. The roughness is controlled by the parameter controlling the aggregation, Ac, so that decreasing values of Ac lead to increase in the roughness of the walls.

Figure 9 presents a timelapse of 3D views of the advection and adsorption of a pulse of tracers in a randomly aggregated slit pore. The tracers start at t = 0 on the left side of the sample (upstream), and the fluid flow goes from left to right. As in the simpler geometry, the tracers follows the fluid flow profile, with higher density at the center of the pore, and lower close to the solid-liquid interface. At long times (see t = 9000 on Fig. 9) we observe a direct effect of adsorption on the distribution of tracers, as the tracer density is lower downstream than upstream -in spite of the fact that the adsorption and desorption coefficients are identical for all adsorption sites. A large amount of solute is adsorbed upstream at the solid-liquid interface, thus the density of free tracers decreases close to the interface. There is therefore a lower concentration of available tracers on the adsorption sites located downstream, the adsorbed density is thus lower.

Quantitatively, we show on Figure 10 the maximum value of the adsorbed fraction over time (the peak value of the elution curve, max γ(t)) as a function of the external force Fext. The different values of aggregation criterion Ac represent three different values of the roughness, where Ac = 1/8 corresponds to the lowest roughness coefficient and Ac = 1/32 to the highest roughness coefficient. In all cases, we observe that Fa decreases when Fext increases i.e. the fraction adsorbed is lower for high velocities of the fluid. The tracers close to the interface are adsorbed at the entrance of the material leading to a lower density close to the adsorption sites. Then, the tracers do not have time to diffuse perpendicularly to the flow to reach the active sites. This effect increase with fluid velocity at constant bulk diffusion coefficient and is, in a broader fashion, dependent on the Péclet number of the fluid flow.

We can see that, in the case of random heterogeneous roughness, we do not observe the effect seen in the ordered crenels: as roughness increases, most adsorption sites remain active and thus higher roughness coefficient does not yield to a lower fraction adsorbed (as was the case for crenels). That is because the "side pockets" of the slit pore created by our aggregation procedure are well-connected to the central area of the slit pore, and have an effective aspect ratio (width over depth) that corresponds to the r < 1 case in the crenelated slit pore. We thus see that having disordered multiscale porosity can in this case be more efficient, in terms of fluid transport and solute adsorption, than a well-ordered porous space of similar dimensions.

Conclusion

We have here extended and applied a Lattice Boltzmann scheme for fluid transport in nanoporous materials, taking into account adsorption of solutes on the solid-liquid interface. We use it to study the kinetics of tracer adsorption and the accessibility of adsorption sites in simple models of hierarchical (or multi-scale) porosity. We validated this model on a simple slit pore geometry, then used it to study a slit pore with grooves (or side pockets) of various aspect ratio. We demonstrate the potential limitations due to inacessibility of some adsorption sites in the presence of deep grooves, highlighting the problem of access of the solute to the mesoporosity (here crenels) in the walls of the material in presence of flow. We then studied the effect of a more disordered type of porosity, through computer-generated models of rough slit pores with controlled rugosity. This opens the way to future work on complex, realistic models of hierarchical porous materials, e.g. coming from reconstructed 3D images obtained by X-ray tomography.

Figure 1 .

 1 Figure 1. Experimental and modelisation techniques employed for porous materials characterization.

Figure 2 .

 2 Figure 2. Schematic representation of the simulation box used in this work.

Figure 3 .

 3 Figure 3. Example of an elution curve γ(t) for a pore with slit geometry, with l = 50∆x, lm = 800δx, Ka = 0 and Fext = 2E -6∆p/∆x.

Figure 4 .

 4 Figure 4. Evolution of the dispersion coefficient K as a function of lm for different values of Fext (2 10 -6 , 4 10 -6 , 6 10 -6 , 8 10 -6 , 10 10 -6 ∆p/∆x) in a slit pore geometry, with l = 50∆x, Ka = 0.

Figure 5 .

 5 Figure 5. Comparaison of the dispersion coefficient computed numerically from our model and the known analytical solution. (a) As a function of the average velocity of the fluid on y direction, in the absence of adsorption (Ka = 0). lx = 1∆x, ly = 1450∆x, lz = 52∆x, lm = 1400∆x. (b) In presence of adsorption, as a function of the ratio between the adsorption coefficient Ka and the desorption coefficient K d , with Fext = 5 10 7 Pa.m -1 . Other parameters are: lx = 1∆x, ly = 1450∆x, lz = 52∆x, lm = 1400∆x .

Figure 6 .

 6 Figure 6. Time evolution of the density of a pulse of tracers in a crenelated pore with liquid flow (flow from left to right). Aspect ratio of the crenels r = 1, with h = 5, w = 5, and l = 20. a: t = 1000∆t; b: t = 10000∆t; c: t = 20000∆t. Other parameters of the simulation are: lx = 1∆x, ly = 120∆x, lz = 32∆x, Fext = 10 9 Pa.m -1 , Ka = 6.04 m.s -1 , K d = 6.04 10 6 s -1 , ∆x = 100 nm.

Figure 7 .

 7 Figure 7. Time evolution of the density of a pulse of tracers in a crenelated pore with liquid flow (flow from left to right). Aspect ratio of the crenels r = 4, with h = 20, w = 5, and l = 20 a: t = 1000∆t, b: t = 5000∆t, c: t = 10000∆t. Other parameters of the simulation are the same as in Fig. 6.

Figure 8 .

 8 Figure 8. Effect of the ratio r between the height and the width of the crenels on fraction adsorbed. The simulation parameters are the same as in Fig. 7.

Figure 9 .

 9 Figure 9. 3D view of spatial repartition of Free density in slit pore geometry with random roughness on the walls. The units of time is the reduced ones (∆t). Input parameters for the simulations: lx = 50∆x, ly = 150∆x, lz = 52∆x, Fext = 10 9 Pa.m -1 , Ka = 6.04 m.s -1 , K d = 6.04 10 6 s -1 , ∆x = 100 nm.

Figure 10 .

 10 Figure 10. Maximum values of the fraction adsorbed as a function of Fext for different roughness. The simulation parameters are the same as in Fig. 9.

Acknowledgements

This work was supported by the ANRT through CIFRE sponsorship 1262/2013. The authors gratefully thank Benjamin Rotenberg and Maximilien Levesque for fruitful discussions.

Supporting Information

The Supporting Information is available free of charge on the ACS Publications website. Description of the random roughness aggregation model. Relations between reduced units and SI units.

Graphical TOC Entry