Discrete modelling of contaminant diffusion in porous media with sorption q Microporous and Mesoporous Materials

A meso-scale model for diffusion of foreign species through porous media is proposed. The model consid- ers diffusion as a continuous process operating on a discrete geometrical structure dictated by the pore size distribution. Local diffusivities and hence pore space connectivity are dependent on the size of the diffusing species and the sorption of those species onto the pore walls. The bulk diffusivity of the medium has been analysed to consider the effects of pore structure alone and in combination with sorption. The chosen medium is bentonite, which is being considered for use as a barrier to radionuclide transport in future deep geological repositories for nuclear waste. Results for transient diffusion of U(VI)-complex through bentonite are presented and very good agreement with experiments is demonstrated. Results for diffusion of larger chemical complexes are also presented to illustrate the effect of reduced pore space connectivity on steady-state and transient transport parameters. Diffusion of larger complexes can be used for experimental validation of the model. The proposed methodology can be used for any micro- and meso-porous material with known distribution of pore sizes. It can be extended to other pore space changing mechanisms, in addition to sorption, to derive mechanism-based evolution laws for the trans- port parameters of porous media. Authors.


Introduction
Pore network models (PNM) of mass transport through porous media are an attractive way to analyse the effects of changing pore space structure on macroscopic transport coefficients, such as permeability [1,2] and diffusivity [3,4]. These transport coefficients are required for engineering analyses using, for example, finite element or finite difference methods. The analysis of changes in transport coefficients with changes in the underlying structure can, using pore network models, provide mechanistic as opposed to phenomenological evolution laws for engineering analyses.
Pore network models idealise the medium as a set of pores, some of which are connected by throats or mass conduits, in which transport is allowed only through the system of throats. In ideal circumstances, constructing a PNM requires a wide range of microstructural information, including porosity, size distribution of pores, size distribution of throats, average pore coordination, and the coordination spectrum: the fraction of pores coordinated by different numbers of throats [5]. For macro-porous materials the required parameters are readily obtainable with current experimental techniques, such as computed X-ray tomography [6,7]. For micro-and meso-porous media the distribution of pore sizes can usually be determined but the resolution of the experimental techniques tends not to be sufficient to segment all the throats and calculate their sizes. However, valuable information on the connectivity between resolved pores, the presence of throats in the PNM sense, is available indirectly from macroscopically measured mass transport. A principal objective of this work is to develop a pore network model which can be constructed with the limited microstructural information available for micro-and meso-porous media; in this case the distribution of pore sizes alone. The pore network model is based on a recently proposed network architecture [5] capable of describing rich sets of pore coordination spectra conforming to macro-porous datasets [6,7].
We have chosen to illustrate the development of our pore network model by considering the properties of bentonite; a clay displaying pore sizes ranging between sub-nanometres and several hundred nanometres [8]. Although the micro-and meso-pores (<50 nm) in this material are roughly 20% of all pores, they have a pronounced influence on the mass transport. As a result, bentonite exhibits very low permeability and the principal mass transport mechanism is diffusion [9]. Bentonite is being considered as a backfill material around containers of high-level radioactive waste in most of the current designs for deep geological disposal of nuclear wastes. It is planned to serve as a critical barrier to radionuclide transport between degrading containers and host rock over very long time scales. It is therefore of substantial interest to repository designers and risk assessors to be able to predict diffusivity of radionuclides in bentonite, particularly when the pore space evolves with time as a result of environmental factors.
Sorption of diffusing species changes the pore space in bentonite and may have a significant beneficial impact on the long-term diffusivity [10,11]. In this work we will differentiate between the macroscopic diffusivity of the system with and without sorption [12]. The latter is usually referred to as the effective diffusivity of the medium, D e , and depends on the pore space structure alone. The former is referred to as the apparent diffusivity, D a , and depends on the sorption kinetics in addition to the pore space structure. A second objective of the work is to develop the model to account for sorption effects and demonstrate its use in deriving time-dependent changes in the mass transport. Factors that could change the pore space to increase diffusivity, such as corrosion or mechanical damage/micro-cracking, are also of substantial interest for the long-term repository behaviour. These are not considered here, but are the subject of other investigations.
In view of the role of bentonite as a barrier to radionuclide transport, our focus is on the diffusion of a specific U(VI) complex. Experimental diffusivity data of this complex in bentonite is available from the literature and will provide some of the validation of our model. This complex is smaller than the experimentally measured pore space and percolates easily through the synthetic pore network model. In this case nearly all network pores are fully coordinated. The model establishes a link between the size of the diffusing particles and the apparent network connectivity. To study the effect of reduced pore connectivity, analysis of diffusion of a larger species is also performed. It is shown further how the small-scale numerical model can be used to perform realistic transient diffusion analysis using appropriately derived film coefficients for model boundaries.

Discrete geometrical space for diffusion
The system of pores in a porous medium can be considered as a set of voids of variable size randomly distributed in space. The coordination of each void with neighbouring voids, or the potential pathways for mass transport, can be determined by construction of the Voronoi diagram around all voids [13]. This tessellates the space into polyhedral cells with one void per cell. The neighbours of a particular void are the voids residing in cells having common faces with the cell of the void of interest. Thus a void has as many neighbours, and hence potential mass transport channels, as the number of faces of its cell.
The process of tessellation can be used to construct pore networks from reconstructed 3D-images of a particular sample of a porous medium [14]. However, networks built this way are irregular and quite specific to the imaged sample. We prefer to generalise the problem to avoid this sample-dependency and allow ourselves to scale the results up to sizes relevant to engineers. To do this, the network has to be homogenised in some way to consider the voids to reside in identical cells. Such cells need to fill the 3D space compactly, and be as close to the ''average'' Voronoi polyhedron as possible in a topological sense. Monte Carlo studies of tessellations around random systems of points have provided the statistics of Voronoi polyhedral cells [15]. From there it can be judged that the regular space-filling polyhedron closest to the average cell in an arbitrary spatial distribution of voids is the so called truncated octahedron, or Kelvin solid.
A novel architecture for pore network construction based on the Kelvin solid has been recently proposed and applied to permeability analyses of porous media with known pore and throat size distributions [5]. The discrete regular cellular network based around the pores is adopted in this work and developed further to construct models for diffusion when only knowledge of the pore sizes is available -a situation typical for micro-and meso-porous materials.
The model geometry is illustrated in Fig. 1. A segment of the cellular tessellation of 3D space is shown in Fig. 1a, together with pores of arbitrary size residing at the cell centres. The pores are thought of as spherical containers with volumes dictated by the experimentally obtained pore size distribution, as described below. The entire porosity of the medium is assigned to the pores. Potential pathways for mass transport in this system of pores are depicted in Fig. 1b as cylinders. These represent throats of variable diffusivities, which may be zero, dictated by the sizes of the coordinated pores and the size of the diffusing species, as described later. The throats are volume-less conduits, allowing mass transport via diffusion in the presence of a concentration gradient between the connected pores. Notional throat diameters are considered in the diffusion modelling to link the size of the diffusing species to throat diffusivity, and to account for sorption effects. A single pore in a cell is shown in Fig. 1c to emphasise that the maximum coordination of a pore in the proposed model is 14. Note, that if S is the size of the cell measured between two square boundaries, the cell volume is V c = S 3 /2, and the lengths of the two type of possible throats are: L 1 = S for throats normal to square boundaries; and L 2 = S p 3/2 for throats normal to hexagonal boundaries.
It is important to remember that the coordination number is not a function of spatial coordinates, but a topological feature of each pore: the number of permeable throats adjacent to the pore. All pores are classified with their coordination numbers to give the coordination spectrum; where the spectrum, in mathematical terms, is a chart describing the structure of the pore space.
The proposed network topology allows for higher pore coordination numbers than models based on cubic lattices [3,4]. Further, the pore coordination spectra that can be generated [5] correspond to published experimental data for macro-porous materials [6,7]. This makes the proposed topology versatile and capable of representing large classes of porous media, with the exception of materials that possess noticeably large fractions of pores coordinated by more than 14 pores. Since there is no experimental evidence for the existence of such materials the model promises to be sufficiently realistic for the purposes of our work.

Experimental pore size distribution
It is not currently possible to observe the pore connectivity of bentonite directly, like many micro-and meso-porous media, but the distribution of pore sizes can be determined experimentally. We use recently reported data [8], where four different sample preparation techniques have been tested and analysed with a number of experimental techniques to resolve pore size distribution. Of particular interest are the results obtained with high pressure frozen samples which realistically reflect the clay pore space under hydro-geological conditions. The data used in this work is the 3D pore size distribution obtained for hydrated bentonite with focused ion beam nano-tomography. These data are presented as cumulative pore volume versus pore radius in [8], revealing total porosity of the sample / = 0.68. For the construction of the pore network model, the cumulative pore volume against pore radius is converted into the cumulative probability function of pore radii, F(r), using a standard statistical technique, see for example [16].
The numerically derived F(r) is shown in Fig. 2 with data points shown as circular symbols. A generator of uniformly distributed random numbers 0 6 p < 1 can be used to assign pore sizes. For a given random number p, the pore size will be r = F À1 (p), which could be determined numerically by linear interpolation between experimental points. For the particular experimental data considered, it was found that F(r) can be fitted accurately with the function: where r 0 is the centre and n the shape of the distribution. The fit is also shown in Fig. 2. Hence, for a given uniformly distributed random number p, the resolved pore size is: The process defined by Eq. (2) is used to distribute pore sizes to all cells in the model. This ensures that the distributed pores follow the same size distribution as the experimentally measured population. From the volume of the distributed pores and the prescribed porosity of the system, /, the length-scale of a discrete network with N c pores is determined by the formula

Local and global diffusivities
The driving force for local diffusion is the concentration gradient between connected pores with resistance provided by the throats. Although the throats are assumed volume-less, they are allocated a geometrical property to assist the description of local diffusion and sorption effects. In our work, the throats are considered to be cylindrical conduits, with notional radii R i,j , where i and j denote pore labels. The mass flow, J i,j , through a throat is propor-  Cumulative probability distribution of pore sizes in bentonite. Result obtained by converting cumulative volume versus pore radius reported in [8]. Continuous fit to data shown is used to distribute pore sizes to model cells.
tional to its diffusivity, D i,j , to its cross-section area, A i,j = p R 2 i;j , to the concentration difference between connected pores, DC i,j = C i -C j , and inversely proportional to its length, L i,j , so that Note that L i,j are the characteristic lengths between pores i and j, and can be either L 1 = S or L 2 = S p 3/2, depending on the throat orientation in the network structure. Since the network consists of discrete pores, the corresponding model is also discrete. The material balance for pore i with volume V i can be written as where the sum is taken over all throats coordinating pore i. Thus the mass transport through the network is formulated as a system of equations for dynamic mass conservation at pores. For given boundary conditions, these are solved with a standard matrix solver to obtain pore and throat concentrations and fluxes.
In order to account for the steric hindrance and frictional resistance within the throat for diffusing species of radius r A , the effective pore diffusivity is defined similarly to [4,17,19] by where D 0 is the free molecular diffusion coefficient of the species. This is calculated with the Einstein-Stokes equation [18] where k B is the Boltzmann constant, T is the absolute temperature, and v d is the dynamic viscosity of the liquid. The use of Eq. (7) is justified by the fact that our model represents the pore space in a discrete manner, with the assumption that it is fully saturated with water and the foreign species diffuse in the water along the available pathways or channels. From Eq. (6), as the solute size approaches the throat radius, the steric hindrance effect will increase, and when r A P R i,j , the solute molecules cannot pass through the throat. This means that for a given pore and throat size distribution, the pore space connectivity will be controlled by the solute size: from low connectivity for large solvent molecules to high connectivity for small solvent molecules.
To illustrate the effect of solute size, two diffusing species are considered in this work. One is a neutral U(VI) complex, Ca 2 UO 2 (CO 3 ) 3 . The molecular size of this complex is r A = 0.524 nm [11]. The other one is a hypothetical colloidal solute with r A = 20 nm. The bentonite is considered fully saturated with water to reflect repository conditions, and the solutes diffuse in the water with molecular diffusion coefficients calculated with Eq. (7) at room temperature: D 0 = 4.66 Â 10 À10 m 2 /s for U(VI); and D 0 = 1.22 Â 10 À11 m 2 /s for the colloid.
The model network occupies the region (0 6 X 6 10S, 0 6 Y 6 10S, 0 6 Z 6 20S) with respect to a coordinate system (X, Y, Z) oriented along the principal directions of the unit cell. This network contains 3539 pores and 22,006 potential throats. The size is sufficiently large to reduce the effect of the random spatial distribution of pore sizes on calculated steady-state diffusivity of the system to less than one order of magnitude. This has been confirmed by a series of random spatial distributions for which the analyses provided nearly identical results. The boundary conditions reflect a particular experimental setup, where solute is supplied with a given concentration at one boundary and flux is measured at the opposite boundary, with the remaining four boundaries not permitting flux. The applied boundary conditions are: fixed concentration C 0 = 5 Â 10 À4 g/m 3 on the boundary Z = 0; prescribed zero flux on the boundaries X = 0, X = 10S, Y = 0 and Y = 10S; prescribed concentration C 1 or film coefficient h on the boundary Z = 20S. Prescribed C 1 are used with steady-state analyses to calculate macroscopic diffusivity using Fick's law: where J m is the total flux through pores on Z = 20S, L m = 20S is the diffusion path, and A m = 10S Â 10S is the cross-section area for diffusion. Selected C 1 are used to derive film coefficients for the transient analyses, as described below. Prescribed h are used with transient solutions to allow for mass outflux, mimicking repetitive extension of the network in the Z direction.
It is clear from the above setup that the mass transport of species of given size through a throat depends crucially on the decision for the hypothetical throat radius. This needs to be determined from the available information on the radii of connected pores, and throat length calculated from the prescribed porosity using Eq. (3). A benchmark for such a decision is the well-known scaling law between macroscopic diffusivity, D m , and porosity, /, expressed by D m /D 0 $ / b , where D 0 is given by Eq. (7) and b can be determined experimentally [18]. We propose the following expression for the calculation of the hypothetical throat radius, R i,j , as a function of the radii of the connected pores, r i and r j , and the throat length, L i,j where a is a calibration parameter to fit a particular scaling law. To illustrate the suitability of this decision, simulations have been performed with the discrete model populated using the experimental pore size distribution in Fig. 2, but with variable porosity values ranging from 0.3 to 0.8, including the value 0.68 measured experimentally. Note that the variation of porosity changes the distance between pores but not pore sizes distributed in the network. The The results in Fig. 3 are important to demonstrate that a wide range of scaling laws can be reproduced accurately by our proposal for hypothetical throat sizes and the methodology for calculation of macroscopic diffusivity. The most widely applied scaling law uses b = 4/3 [19], which is reproduced with a = 1 in our model. This value is adopted for the rest of the work. Note that the D m obtained with these simulations is in fact the effective diffusion coefficient of the system, D e . It incorporates the effects of pore sizes and connectivity, but excludes non-linear effects, such as sorption, which may change the pore space with time.

Sorption effects and apparent diffusivity
Compacted bentonite offers a very high sorption capacity for radionuclides, suggesting that long-term diffusivity of such species can be significantly affected by sorption [12]. The molecules of diffusing species M diffuse in the pore liquid and are adsorbed onto the walls of the throat. The fractional saturation of the active sites, h, is given by an adsorption isotherm, h = h (C M ), where C M is the local concentration of species M in solution and h is the concentration in the solid. Different expressions are commonly used for correlation and prediction of equilibrium isotherms [19]. In most laboratory-scale or in situ experiments, linear isotherms or a constant distribution coefficient are used to describe the adsorption process [20].
However, under constant chemical conditions the adsorption isotherms of actinide elements and transition metals are often nonlinear and best characterised by the Freundlich isotherm [21]. For example, experiments with U(VI) have shown that the adsorption isotherms at different pH values are nonlinear and each can be described by a Freundlich isotherm, h = K f C 1/n [9]. In the current work a recently reported isotherm for U(VI) in compacted bentonite [11] is adopted. This was obtained with experiments conducted at solid/solution ratios of 94-100 g/L, pH of 8.0 ± 0.1, Pco 2 = 10 À3.5atm, 22.5°C. After 15 days, the estimated Freundlich isotherm was h = 18.6 C 0.81 . The conditions of compaction and alkalinity are close to those expected in the bentonite when used as a barrier in deep geological repositories for nuclear wastes. Since we also wish to study diffusion of larger colloidal particles, for which no sorption isotherm is available, the U(VI) sorption isotherm is used for the simulations of diffusion with sorption for the colloid.
An estimate for the effect of adsorbate obstruction on diffusivity can be obtained by reducing the radii of the throats as a function of the concentration of the diffusing species. In principle, systems with dilute adsorbate region and systems with high adsorbate concentration can be considered separately [3]. As compacted bentonite offers a very high sorption capacity for radionuclides, an estimate for the obstruction effect at high adsorbate concentration is used here. This can be obtained by considering the adsorbate to be uniformly smeared into a layer of thickness, t. The smeared layer assumption becomes increasingly more realistic as h increases. The thickness of the adsorbed layer, t, can be calculated from the equation [3]: where R is the throat radius, r M is the radius of the diffusing species, and h is the adsorption isotherm. Thus, the radius of a throat after adsorption at given concentration, C A , and hence h, becomes R ⁄ = R À t. The concentration C A in the throat is taken to be the average between the concentrations of species A in the connected pores. This is a reasonable assumption considering the linear approximation of the concentration gradient within the throat and the fact that local morphology effects are outside the scope of the proposed model. The flux of species A through a throat connecting pores i and j, with sorption taken into account is therefore given by: This expression has been used to simulate diffusion through the system with sorption, providing an apparent diffusion coefficient, D a . For comparison with available experimental data [10], the system with a = 1 has been solved for porosities / = 0.62-0.64. The results of our simulations, obtained with 10 random distributions of pore sizes within the system, were in the range D a = 2.83 Â 10 À12 to 2.85 Â 10 À11 m 2 /s. Experimentally obtained values were D a = 1.16 Â 10 À12 m 2 /s to 2.30 Â 10 À12 m 2 /s [10]. Clearly, with some of the random distributions of pore sizes the model predictions are in very close agreement with experimental values. With other distributions the outcomes are in poorer agreement. This could be attributed to potential differences of pore size distribution between the bentonite used in [10] and the imaged bentonite that formed the basis of the model. Further, the sorption isotherm applied in the model was determined from bulk measurements [11] and could be inaccurate when describing solute sorption on throat walls. It is possible that the sorption rates within confined and curved spaces are different from those measured macroscopically. Further studies to clarify whether such differences exist could be performed with appropriate atomic scale models. The comparison suggests stronger sorption effect, assuming that the pore size distribution used in the model and the experiments were identical. Finally, the size of the model was smaller than the experimental test specimen and the model was solved for steady-state conditions with fixed inlet and outlet concentrations. In the absence of data for microstructure and apparent diffusivity for one and the same system, further analysis of the factors creating the differences is not possible. However, the potential size effect is addressed by the introduction of film coefficient at the outlet boundary, intended to represent model extension in the direction of diffusion. This is discussed in the next sub-section. In the remaining of the work sorption is modelled with the isotherm h = 18.6 C 0.81 .

Film coefficients
In many lab-scale and field-scale situations, fixed concentrations are usually prescribed on two opposite sides of the system, where the concentration is constant on one side and the concentration on the other side is zero. However, analysis of a model as large as the experimental one is impractical. In this study, we use a much smaller model but large enough to represent the pore space of the material. In order to simulate practical situations where the far-field concentration, C 1 , approaches zero, we introduce a film coefficient, h, on the boundary Z = 20S. It should be noted that the far-field concentration is the down-stream concentration far away from the outlet surface. The film coefficient introduced in this study is used to allow for appropriate out-flux of species through the boundary Z = 20S, so that the concentration on the outlet boundary can change with time. The film coefficient relates the flux J m and the concentration of solute C 1 to the far-field concentration, C 1 , via J m = h (C 1 -C 1 ). Noting that for a given C 1 , the steady-state solution of the system obeys J m = K (C 0 -C 1 ), where K = D m A m /L m , the film coefficient can be written as a function of C 1 : Note that for the linear problem of diffusion without sorption, D m and hence K is independent of the prescribed concentrations and depends solely on the pore network structure and size of diffusing species. Hence, for diffusion without sorption the film coefficient dependence on C 1 can be determined with a single steadystate solution to calculate K. For the non-linear problem of diffusion with sorption, K is a function of (C 1 /C 0 ), which can be derived numerically by performing a number of simulations with a fixed inlet concentration, C 0 , and a variable outlet concentration, C 1 .
Such simulations were performed for both smaller and larger diffusing species in one and the same network system: constant pore size distribution and initial porosity. The results for the film coefficient, h, as a function of (C 1 /C 0 ) are presented in Fig. 4. For these plots it is assumed that C 1 = 0.
The results in Fig. 4 show how the film coefficients decrease with increase in the solute concentration at the outlet boundary. Considering transient diffusion through the pore system, the film allows for large out-flux for low outlet concentrations which decreases with the increase of the outlet concentration. Since the solutions are obtained with identical pore distributions, the application of the film coefficient as a boundary condition represents a repetition of the same structure in the Z-direction enough times that the far-field concentration is maintained at zero. Clearly, the sorption effect is in the reduction of the film coefficient from the case without sorption. This is due to the obstructive effect of adsorbed substance, which changes local diffusivities as well as connectivity as demonstrated in the next section. It can be observed that the obstructive effect of sorption of larger diffusing species yields substantially larger reduction of the film coefficient. Whilst the general trend is expected, the actual film coefficient values cannot be obtained without knowledge of the pore system structure and the size of the diffusing species. The film coefficient dependencies on the concentration at the outlet boundary are used in the corresponding transient analyses presented in the next section. The film coefficients are used to apply generalised flux boundary conditions via J m = h (C 1 -C 1 ), thus incorporating their effect in the initial boundary value problem formulation. The problem becomes non-linear, even in the absence of physical non-linearity, such as the sorption, and is analysed with incremental non-linear solver.

Results and discussion
The results presented hereafter were obtained with the experimentally determined porosity of bentonite, / = 0.68, and one random distribution of pore sizes obeying Eq. (2) with r 0 = 86 nm and n = 2.646, see Fig. 2. This determined a model length scale S = 520 nm from Eq. (3). Throat radii were calculated with Eq. (8) with a = 1. Local fluxes were defined with Eq. (10), where R ⁄ is either the initially assigned radius when sorption is not considered or the current radius when sorption is accounted for. The transient solutions were performed with fixed concentration C 0 = 5 Â 10 À4 g/ m 3 on the inlet boundary Z = 0 and C 1 -dependent film coefficient h on the outlet boundary Z = 20S, as derived in Eq. (11) and shown in Fig. 4.
Structural and physical parameters were monitored during the simulations. The structural parameters were the pore space connectivity, expressed by the pore coordination spectrum, and the average pore coordination number. The idea behind the introduction of structural parameters is to investigate the link between the pore space structure and the macroscopic transport behaviour of the medium. The transport behaviour of the pore space is characterised by the model diffusivity. The initial values of the structural parameters for the two diffusing species studied are shown in Fig. 5. The pore coordination spectra represent fractions of pores, f c , coordinated by given number of throats permitting mass transport, c, where c varies between zero for isolated pores and 14 for fully coordinated pores. The average pore coordination number, z, is the sum of all products (c Â f c ) for c = 0. . .14. This is depicted in the boxes of the two plots in Fig. 5. For the smaller diffusing species, U(VI)-complex, the pore space initially permits transport through practically all possible throats. The spectrum for the  colloidal species is more complex, with a substantial fraction of isolated pores and average coordination number close to those obtained with resolved throats in macroporous geological materials [6,7].
For diffusion without sorption the pore coordination spectra and average numbers shown in Fig. 5 do not change. When sorption is accounted for, the spectra evolve and the informative way to illustrate this is to give the time dependence of the average coordination number z. This is shown in Fig. 6 for the two diffusing species studied. The results provide an averaged account for the reduction of pore space connectivity with continuous solute adsorption at throat walls, which leads to progressive blockage of throats. The process appears to be evolving the structure towards a lower limit of constant pore connectivity. The outcome suggests that in order to be comparable to experimentally obtained values, the apparent diffusivity of the model needs to be determined within the tail of the connectivity evolution.
To illustrate the effect of sorption, the change of the systems diffusivities with time is shown in Fig. 7. It is clear that the reduction of connectivity in both systems yields reductions of the longterm apparent diffusivity. In the case of U(VI) diffusion, Fig. 7a, the reduction of diffusivity is about two fold. If applied to the system with 62% porosity measured experimentally in [10] the model predicts diffusivity values closely matching the experimental values. Notably, the initial apparent diffusivity of the colloidal species, Fig. 7b, is three-orders of magnitude lower than the U(VI), whilst the size of this solute is about 40 times larger. This is an illustration of the strong non-linear effect of the system connectivity. Further, the sorption reduces the diffusivity of the colloid by nearly one order of magnitude. This reflects the faster reduction rate of the coordination number observed in Fig. 6.
Apart from the ability of the model to predict macroscopic diffusivity with good accuracy, it can provide information about the invasion of the solute in the pore system with time. In principle, such invasion can be studied experimentally using 4D computed tomography of a porous specimen with a diffusing tracer element. The technique is still under development and the outcomes will be used for enhanced model validation. The experiments will provide an image of the domain where the concentration of the tracer, C, is above certain fraction of the inlet concentration, C 0 , at each time point. The imaged domain is discrete containing the pores and the channels, or throats, forming invasion pathways. The parameters that can be used to characterise the discrete domain evolution and be compared to the model predictions are: the invasion depth; the total length of the invaded pathways; and the total pore space volume containing tracer with concentration greater than C. Clearly, these are still macroscopic parameters with respect to the pores and throats. However, the model is intended to describe and simulate the macroscopic behaviour by averaging the localised pore space features. Hence comparisons of diffusion behaviour at pore/throat length scale between model and experiment are out of the scope of the validation process.
In preparation for the validation process, model results are presented for a select tracer solute concentration C = 50% C 0 . The results are from simulations with and without sorption as different tracers may be adsorbed or not. Hence, comparison between model simulations and 4D imaging could also provide indirect information about sorption kinetics. The evolution of solute penetration depth is shown in Fig. 8 for the two types of diffusing species. The penetration depth is the maximum distance from the inlet where the concentration of the diffusing species is greater than the prescribed C. This is normalised by the total model length, the distance between the inlet and outlet boundaries, and plotted as relative penetration. A value of one is reached when the  concentration at the outlet boundary becomes greater than the prescribed C. In the prospective experiments this will mean that the detectable tracer concentration will be measured somewhere at the outlet boundary. Note that the time for full penetration of the colloidal species is less than one order of magnitude larger than the time for U(VI) species. This suggests that percolating pathways exist for both species and the difference in the penetration times can be attributed predominantly to the different molecular diffusivities of the two species. One can also note that the rate of penetration changes with time, from higher during the first 1-2 s to moderate rate controlled by the pore space topology, to zero at full penetration. This illustrates the importance effects of the structure on the properties, which cannot be analysed using classical continuum description of the porous medium.
The evolution of the invaded path length is shown in Fig. 9 for the two types of diffusing species with and without sorption. The invaded path length is the sum of all lengths of throats with concentration larger than the prescribed C. This length is normalised by the total length of throats, or potential diffusion pathways, and plotted as relative length. Because of the almost complete connectivity of the system with diffusing U(VI), the whole system of pathways is invaded almost as quickly as the specimen is penetrated in the absence of sorption, compare Figs. 8a and 9a. The sorption effect on invaded path by U(VI) is predominantly manifest by blocking some throats, i.e., reducing the length of diffusion pathways, whilst the time for pathways invasion is similar to the time for penetration with sorption. The result for the colloidal species is very different as the time to invade all throats permitting diffusion is about two orders of magnitude larger than the time for full penetration, compare Figs. 8b and 9b. The sorption effect is similar to the case of U(VI) diffusion.
The evolution of the invaded pore volume is shown in Fig. 10 for the two types of diffusing species with and without sorption. The invaded volume is the sum of all volumes of pores with concentration larger than the prescribed C. This volume is normalised by the total accessible pore space volume and plotted as relative volume. The results show an interesting behaviour, different from the penetration depth and invaded length. The jump-like changes at certain times indicate that the smooth invasion along pathways observed in Fig. 9 leads to sudden breakthroughs in the invaded pore space. The breakthroughs reflect the complexity of the pore space. Whilst the invaded length is changing smoothly as it depends on the lengths of invaded throats, the invaded pore volume experiences jumps due to the inclusion of large pores in the penetrated paths. This again shows the effect of the pore space structure on the meso-scale transport behaviour of the medium. The time for invading the entire pore space is approximately the same as the time for penetration for U(VI), but is about two orders of magnitude larger for the colloidal species, similar to the time for invading the accessible pathways as expected.
The results presented here capture the transport behaviour at a ''meso'' length-scale, which is the length of several to tens of interpore distances. They illustrate that at this length scale the effect of the pore space geometry and topology is strong, leading to different evolutions of penetration depth, invasion length and volume. Such behaviours cannot be observed with standard bulk diffusion experiments and our plan is to validate the results with 4D imaging.
The comparison between the smaller and the larger diffusing species is of practical importance for planning future experiments. The penetration time and the times for volume and length invasion . Penetration depth is the maximum depth at which solute with concentration C = 0.5 C 0 is found. Depth is normalised with total model length, i.e., distance between inlet and outlet specimen boundaries. of the smaller species suggests that the use of a solute of comparable size as a tracer would not be practical for experiments, as the current experimental techniques may not allow for sufficient number of 3D images to be obtained over a short enough time interval to validate the model. The use of a tracer with a molecular size comparable to the larger species studied here will allow a longer time interval over which 3D images can be taken and smoother evolution analysed. The difference in the time-scales for penetration and path length or volume invasion provides an excellent opportunity for experimental validation of our model. Simultaneous matching of the three parameters between experiments and simulations at different times will strongly support the model.
In this work, pores of different sizes are assigned to all cells of a regular space tessellation. The sizes of the assigned pores follow experimentally determined pore size distribution, but their positions are fixed by the predefined cellular topology. The results could therefore be thought as mean-field solutions, albeit the averaging is done in advance by selecting a physically realistic topology. The length scale of the model, in our case the cell size, represents the average inter-pore distance of the real pore space, when calculated with specific porosity. Hence, the results presented here are applicable to materials with the same pore size distribution and average inter-pore distance.

Conclusions
A discrete model for diffusive transport through porous media has been proposed with the ambition to provide a bridge between the short length scale of individual pores and the long length scale of interest to engineering applications. The methodology has been developed subject to two constraints: limited experimental information about the pore space; and limited size of the modelled pore space. The former has been tackled by selection of local diffusivities as functions of pore sizes and total porosity alone. The need for connectivity data which can be unobtainable for micro-and meso-porous materials is thus avoided. The latter has been tackled by the introduction of boundary film coefficients that allow us to prescribe a far-field concentration by hypothetical repetition of the modelled pore space. In addition the idea has been extended to account for the adsorption of the diffusing species and study its effect on the long-term diffusivity of the system.
The model has been applied to bentonite using a reported experimental pore size distribution. With a combination of film boundaries and sorption kinetics the model has been shown to give macroscopic diffusivity in very close agreement with reported experimental values for radioactive contaminant with small molecular radius. Whilst the agreement for this macroscopic parameter is promising, the performance of the model requires validation at a ''meso-scale'', with measurements involving diffusion through several tens to hundreds of inter-pore distances. From the results presented for the invasion of the small contaminant with time, it has been concluded that diffusion experiments with tracer of comparable size would not be able to provide validation data. Therefore results for larger diffusing species, hypothetical colloidal molecules, have also been presented. Through these, the strong effect of the pore space connectivity on the system diffusivity has been demonstrated. Further, it has been concluded that diffusion experiments with a tracer of similar size could be used for the required advanced model validation at the meso-scale.
Although this work has used a particular material to illustrate the methodology and its application, the model is not limited to bentonite. The proposed approach can be applied to any other micro-and/or meso-porous material after analysis of the pore size distribution and total porosity. Further, the model offers clear separation between effective diffusivity, dependent only on the pore space structure, and apparent diffusivity, resulting from local non-linear effects. This can be used in conjunction with experiments to investigate internal sorption by appropriate selection of non-adsorbing and adsorbing tracers. The approach is not limited to sorption effects and could be used for a range of other pore space changing mechanisms. For example, our immediate plans are to study the effect of micro-cracking due to thermo-mechanical loads on the diffusivity of the medium with application to long-term behaviour of bentonite in deep repository for nuclear wastes. Such a coupled problem is of significant interest to the oil and gas industries, particularly when hydraulic fracturing is used to enhance hydrocarbon recovery. In this application, the evolution of permeability rather than diffusivity of the medium is required, but the methodology developed in this work would not change. The most promising area for the application of our approach to such problems is the shale gas extraction. Further, the model can be extended to investigate possible use of different porous media for uranium separation from waste as well as for analysis of multicomponent systems. Analyses with colloidal-porous media systems could be of interest to the pharmaceutical industry.
An important element of future model development is the better understanding and the more realistic description of sorption strength. It is in the plans of the authors to use atomic scale modelling techniques for the analysis of sorption in order to be able to derive efficiently sorption isotherms and improve the predictive capabilities of the model.