Measurement and modelling of reactive transport in geological barriers for nuclear waste containment

Abstract

Compacted clays are considered as excellent candidates for barriers to radionuclide transport in future repositories for nuclear waste due to their very low hydraulic permeability. Diffusion is the dominant transport mechanism, controlled by a nano-scale pore system. Assessment of the clays' long-term containment function requires adequate modelling of such pore systems and their evolution. Existing characterisation techniques do not provide complete pore space information for effective modelling, such as pore and throat size distributions and connectivity. Special network models for reactive transport are proposed here using the complimentary character of the pore space and the solid phase. This balances the insufficient characterisation information and provides the means for future mechanical–physical–chemical coupling. The anisotropy and heterogeneity of clays is represented using different length parameters and percentage of pores in different directions. Resulting networks are described as mathematical graphs with efficient discrete calculus formulation of transport. Opalinus Clay (OPA) is chosen as an example. Experimental data for the tritiated water (HTO) and U(VI) diffusion through OPA are presented. Calculated diffusion coefficients of HTO and uranium species are within the ranges of the experimentally determined data in different clay directions. This verifies the proposed pore network model and validates that uranium complexes are diffusing as neutral species in OPA. In the case of U(VI) diffusion the method is extended to account for sorption and convection. Rather than changing pore radii by coarse grained mathematical formula, physical sorption is simulated in each pore, which is more accurate and realistic.

---

1 Introduction

Simulation of reactive transport in both natural (soils, rocks, etc.) and industrial (concrete, bioengineered tissues, etc.) porous media is highly important to a variety of scientific and engineering applications.^1–4 These include environmental contaminant transport,¹ carbon dioxide (CO₂) storage,² bioremediation,³ and stimulation of petroleum reservoirs.⁵ Most applications require further understanding of multiple processes (diffusion, convection, and reaction) and scales (microscopic and macroscopic). Because continuum reactive transport models involve averaging over length scales much larger than typical grain sizes, they are unable to resolve heterogeneity at smaller scales. This heterogeneity has been found important for reactive transport modelling at larger scales.⁶ Progress can be achieved by developing physically based and microstructure-informed models suitable for such analyses. Key initial requirement is that such models are able to predict measurable transport properties at a macro-scale (engineering, geological) from measurable pore space characteristics describing structures at a meso-scale (tens of inter-pore distances).

Pore network models (PNM) are amongst the appealing approaches that provide a suitable description of mutable pore space structures. Such models have been used to describe conservative as well as reactive transport in saturated⁷ or multiphase and unsaturated porous media.^8,9 Moreover, their use has been expanded to complement the application of percolation theory to porous media¹⁰ and to simulate the effect of microbial growth on substrate transport and system permeability.³ The pore space is approximated by a set of sites and a set of bonds connecting some of the sites. The sites positions may form a regular lattice (regular PNM) or reflect known length variations (irregular PNM). Mass transport is allowed through the system of bonds.

Conceptually, PNM are scale indifferent, i.e. they can be applied to any length-scale interval where the structure of the pore space has been experimentally characterised. For example, if a particular experimental technique allows for characterising pore features of sizes between 0.1 μm and 10 μm, then the corresponding PNM will take these into account. The elements of PNM, sites and bonds, are abstract and can be related to the measurable features in different ways depending on the available information.

Historically, pores were related to the sites and pore throats were related to the bonds of regular PNM.^11–13 If such correspondence is to be statistically representative of the material modelled, then sufficiently rich experimental information is required: shape and size distribution of pores and throats, as well as the pore coordination spectrum, i.e. percentages of pores coordinated by different numbers of throats.^13,14 These can be obtained from structures with distinguishable pores and pore throats. One established way is analysing 3D images from synchrotron X-ray tomography.^15,16

It is worth mentioning, that 3D images can be used directly to construct irregular PNM.¹⁷ These allow for validating physical assumptions for flow simulations using 4D imaging of mass transport. However, such irregular PNM are sample-specific and potentially not statistically representative. A regular PNM constructed within a larger volume allows for capturing statistics from a number of imaged samples, i.e. improved statistical correspondence, and for calculating transport at distances closer to the engineering application length scale.

X-ray computed tomography is presently limited to resolving macro-porosity (features larger than 50 nm). However, the resolution of this experimental technique is not sufficient to identify throats for systems dominated by meso-porosity (features between 5 nm and 50 nm) and micro-porosity (features smaller than 5 nm). For pore network construction, this means that pore connectivity data cannot be extracted as throats cannot be identified.

One possibility is to develop a regular PNM with pores from experimental distribution located at sites and notional throats between each pair of neighbouring pores¹ based on a topologically proper lattice proposed by Jivkov et al.¹³ Such a construction allows for the calculation of the lattice length scale from known porosity and pore size distribution, but the pore connectivity is constant. The diffusivity of a notional throat depends on the sizes of the connected pores, the size of the solute molecules, and the interaction of the solute with the pore walls. The approach provided insights into the effects of pore size distribution, solute size, and sorption on diffusivity.¹ It was noted, however, that the spatial randomness of pore sizes yielded variation of calculated diffusion coefficient larger than experimentally measured variability. This indicated that in the real material there was a strong connectivity effect, missed by the model due to lack of detailed information. In addition, diffusivity anisotropy emerging from material texture (preferred mineral orientations) could not be captured.

Another possibility is to assume variable connectivity.¹⁸ While this was more realistic in terms of pore space topology, the observable data were insufficient to establish rigorously a lattice length scale. With variable length scales, improved agreement with the variability in experimentally measured diffusivity was demonstrated, including diffusion anisotropy in clay.

The two approaches to tackle incomplete pore space information, namely predefined connectivity for calculable length scale,¹ and undetermined length scale for realistic connectivity,¹⁸ suffer from the lack of an additional constraint. This cannot be found within the pore space information and requires the consideration of the solid phase structure, e.g., the shape and size distribution of mineral grains. A methodology for incorporating the structure of the solid phase is proposed here. It improves substantially the realism of the constructed PNM, both in terms of geometry and topology. The proposal has added benefit that the constructed PNM can be paired directly to lattice models of the solid phase developed for analysis of damage evolution via micro-cracking.^19–21 This facilitates mechanistic investigations of combined mechanical–thermal–chemical–biological effects on diffusivity, which will be discussed further in the paper.

Typical pore network models idealise the pore space with simple geometries, i.e. throats are assumed to have square or circular cross sections. This makes it difficult to simulate processes like sorption, which may change pore geometries by physical clogging. Recent advances allowed to model pore geometries with irregular cross sectional shapes.⁸ The influence of cross sectional shapes on diffusion and sorption will be analysed in this work.

For safety assessment of repositories for high-level nuclear waste, in particular the migration behaviour of radioisotopes with long half-lives (e.g., ²³⁸U, τ_1/2 = 4.468 × 10⁹ a; main fraction of spent nuclear fuel) is in the focus of research. Opalinus Clay (OPA) is discussed as potential host rock for nuclear waste repositories, where it is considered as one of the barriers to radionuclide transport, because it exhibits very low hydraulic conductivity.²² This makes molecular diffusion the main transport mechanism of radionuclides through the barrier. In the present work, OPA is selected to demonstrate the development of the new pore network model. The pore space information for model construction is adopted from Keller et al.^23,24 and Nagra.²² Validation data on macroscopic diffusivities come from recent literatures.^25–28 Although OPA has a very low permeability, the effects of velocities on mass transport are also analysed to demonstrate how the methodology will be applied to other porous media. In the diffusion experiments so far the neutral Ca₂UO₂(CO₃)₃(aq) complex was assumed to be the main diffusing species through OPA,²⁵ because it was spectroscopically verified to be the dominant U(VI) species in the experimental source reservoir solution. The model described in the present study shall clarify whether U(VI) diffuses as cationic, anionic, or neutral species through OPA.

While the focus of this work is on the reactive transport of radionuclides through Opalinus Clay, the modelling approach proposed is applicable to analysis of transport phenomena in other porous media, where the pore space cannot be fully characterised by current experimental techniques. Examples include reactive transport of species in micro-porous and meso-porous materials which may not be easily studied by classical in-diffusion or through-diffusion methods, biomass growth in bioengineering materials and dissolution–precipitation processes in carbon capture and storage.

2 Geometric modelling and data

2.1. Cellular basis for pore network construction

Statistical analysis of the neighbourhoods of arbitrary features, randomly distributed in 3D space, has shown that the truncated octahedron represents the average neighbourhood.²⁹ The truncated octahedron was the unit cell of a regular space tessellation, representing the neighbourhood of one type of measurable feature, a solid particle²⁰ or a pore,¹³ which has been used to analyse deformation and failure of quasi-brittle media^19,30,31 and transport in porous media,^1,18 respectively. Considering the need to couple the deformation and the transport models in the future, we choose the truncated octahedron as cellular basis for network construction. This has six square faces normal to three perpendicular directions, say (1, 0, 0), (0, 1, 0), and (0, 0, 1), and eight hexagonal faces normal to the corresponding octahedral directions, (±1, ± 1, ± 1).

The material is reduced to a mathematical graph by placing particles or grains at cell centres (interiors). This is illustrated in Fig. 1(a) for cells with equal distances between the three pairs of square faces, a setup used in the previous works. In order to represent texture anisotropy, three different length parameters, S₁ [m], S₂ [m], and S₃ [m], are used to measure the distances between the square faces in directions (1, 0, 0), (0, 1, 0), and (0, 0, 1), respectively. All other distances are then dependent on the three length parameters, e.g., the volume of each cell is V_c = S₁S₂S₃/2 [m³].

Fig. 1 Pore space representation: (a) unit cell illustrating potential diffusion paths (bonds, yellow and red) in the neighbourhood of central particle (green); these join neighbouring cell faces and show where elongated pores may be assigned to the experimental pore system information; (b) a fragment of a lattice model for deformation and failure analysis, illustrating sites on particles and bonds between neighbouring particles.

In an assembly of N_c [—] cells, each cell is assigned to a particle with different radius, r_i [m]. The assignment of particles terminates when their accumulated volume fraction equals to the experimentally measured value. From this requirement the volume of a cell assembly is calculated by:

where ϕ [—] is the measured volume percentage of particles. The calculation of the three length parameters from eqn (1) depends on the selection of their ratios used here to represent texture. A non-textured medium has S₁ = S₂ = S₃.

Each particle (site, cell centre) is connected to its 14 neighbouring particles by bonds representing relative deformations between cells, illustrated in Fig. 1(b). This graph can be analysed very efficiently using discrete exterior calculus to investigate the response of the material to loading, including micro-crack generation, competition, and coalescence.³² The presence of pores in cell neighbourhoods is reflected in pore-size-dependent reduction of bond strengths. Further, the assigned particle sizes may also contribute to spatial variation of bond properties within the model. However, for the construction of the PNM the particular particle sizes do not play a role; only the length parameter calculated by eqn (1) is of essence.

The pore network construction is as follows. Firstly, a new skeleton is formed using sites at the centres of cell faces and bonds between neighbouring faces. This is illustrated in Fig. 1(a) on a single cell with S₁ = S₂ = S₃ for clarity. The bonds represent potential transport pathways, i.e. they show the positions where elongated pores are allowed to reside. Note that there are no bonds cutting through cells, hence the pore system can follow paths around particles. The elongated pores can be assigned with different cross sectional shapes, such as circular, equilateral, and square. Since pores reside on bonds, the sites at face centres represent pore junctions – volume-less containers redirecting mass transport.

Notably, this skeleton allows for pores with maximum coordination of 8 or 12 depending on whether the bond connects square to hexagonal or hexagonal to hexagonal boundary. Importantly, the length parameters of the skeleton, determining the lengths of possible pores, are now dictated by the solid phase structure.

2.2. Pore shape

A key characteristic of real porous media is the angular form of pores. It has been demonstrated that having pores with a circular cross section, and thus single-phase occupancy, causes insufficient connectivity of the wetting phase and as a result poor representation of experimental data.³³ Angular cross sections retain the wetting fluid in their corner and allow two or more fluids to flow simultaneously through the same pore. Pores which are angular in cross section are thus a more realistic model than the commonly employed cylindrical shape. The shape of an angular pore cross section is prescribed in terms of a dimensionless shape factor, G, which is defined aswhere A [m²] and P [m] are the area and the perimeter of the cross section,³⁴ respectively. The shape factor replaces the irregular and complicated shape of a pore by an equivalent circular, triangle, or rectangle shape. The value of shape factors ranges from a minimum of zero corresponding to a slit to a maximum of 0.08 corresponding to a circular cross section. Values of shape factors for triangular cross sections vary from 0 to 0.048 (with maximum for equilateral triangle), and rectangular cross section from 0 to 0.062 (with maximum for square shape).

2.3. Assignment of porosity

Pores with the desired characteristics, e.g., porosity and pore size distribution, will be generated in a particle or grain structure that best represents the measured material's size and orientation. Usually, the pores in porous material can be divided into several groups, such as macro-pores, meso-pores, and micro-pores. For example, there exist meso-pores and micro-pores in the material. Firstly, meso-pores are assigned to available bonds of Fig. 2(a) according to the experimental data. An available bond is a bond with no previously assigned pore. The volume of all allocated meso-pores is required to equal experimentally determined meso-pore volume, i.e. whenwhere V_i [m³] is the pore volume, θ_mes [—] is the volume fraction of meso-pores. The assignment of meso-pores can be in different directions or constrained to a selected direction, to simulate the anisotropy and heterogeneity. In the selected direction, the pores are randomly assigned to bonds. All bonds in the chosen direction are available for meso-pore assignment.

Fig. 2 Illustration of pore assignment to diffusion site-bond model: (a) exclusive assignment of meso-pores along (1, 0, 0) with different pore cross sectional shapes taken as the clay bedding direction; (b) pore network model after assignment of meso-pores (red) and micro-pores (yellow) – micro-pores allowed in all directions. Examples are given on lattice with equal length parameters for bedding and out-of-bedding directions. (c) Schematic depiction of meso-pores and micro-pores assigned to bonds.

Secondly, micro-pores are assigned to bonds not already occupied by meso-pores. The micro-pores are also assumed to be of different shape and each new pore assignment increments the pore volume according to selected cross section area and bond length. The process terminates when the total pore volume fraction, from meso- and micro-pores, attains the experimental porosity. Notably, micro-pores are allowed at bonds in all directions, including bonds in the bedding direction not occupied by meso-pores, and more than one micro-pore can be assigned to a bond if required to reach prescribed porosity, illustration in Fig. 2(c). A set of micro-pores assigned to a bond can be considered as providing an “effective” transport between two pore junctions, which is necessary to create transport paths in the otherwise non-percolating meso-pore system. Importantly, they act as sieves for transport of species with different molecular sizes.

An example of a PNM after assigning meso-pores and micro-pores with different pore shape sets is shown in Fig. 2(b). The micro-pore sets are presented as bundles with scaled-up cross-sections and the network is constructed with equal length parameters in the bedding and out-of bedding directions for improved illustration. In the simulations different length parameters in bedding and out-of bedding directions were used to represent texture.

2.4. Experimental characterisation data

OPA has anisotropic transport properties due to preferred orientation (or texture) of clay minerals attained during sedimentation and compaction.³⁵ Specifically, solute species diffuse slowly perpendicular and fast parallel to the bedding plane. Our goal is to construct a regular PNM, which calculates diffusion in line with macroscopic observations, and can be linked to models for deformation and failure of the solid phase, based on available structural data.

Regarding the pore space, we make use of the 3D analysis of OPA reported by Keller et al.^23,24 and Nagra.²² The pore space was divided into meso-pores and micro-pores based on the resolution of different experimental techniques in our work. The meso-pores were resolved by Focused Ion Beam nano-tomography (FIB-nt) with elongated pores in the bedding plane. These pores had approximately cylindrical shapes with diameters >10 nm and were largely isolated (i.e. solute species cannot percolate through the sample). The meso-pores occupied approximately 1.8 vol% of the sample, i.e. their volume density (porosity) was θ_mes = 0.018. The micro-pores were resolved by N₂ adsorption with pore sizes <10 nm, but were unresolved by FIB-nt. The analysis of N₂ adsorption gave a total physical porosity of around 11.5 vol%, of which 9.7 vol% were micro-pores. The volume density of micro-pores is thus θ_mic = 0.097.²³ However, micro-pores narrower than 2 nm cannot be probed by N₂ adsorption. For OPA from Mont Terri, most estimates of the total porosity are on the order of 16 vol%.^25,36,37 In this work, we assume that micro-pores smaller than 2 nm occupy 4.5 vol% and follow the same distribution as pores occupying 9.7 vol%. These measurements show that the fine-grained clay matrix contains an extensive micro-pore network. The meso-pore system can be seen as part of the effective “flow relevant” porosity that surrounds the clay mineral particles and is interconnected by micro-pores.

The two measurements were combined into a single “cumulative pore volume fraction – pore radius” curve given in Fig. 3(a).^18,23 This shows that the geometrically resolved meso-pores (symbols) and the adsorption-based curve coincide in their common length interval, with meso-porosity θ_mes = 0.018 and total porosity θ = 0.16. The radii are derived with the assumption of cylindrical pore shapes truthful to experimental evidence for meso-pores²³ and approximate for micro-pores. For the purpose of model construction the experimental distribution from Fig. 3(a) is re-evaluated as cumulative probability separately for meso- and micro-pores. These are shown in Fig. 3(b) and (c), respectively. The use of the cumulative probability is to produce sizes from the measured distribution in Fig. 3(a) using a generator for uniformly distributed numbers 0 ≤ p < 1. This is illustrated by the red arrows on Fig. 3(b) and (c); more details on statistics are given by Jivkov and Xiong.¹⁸ For pores with other shapes such as square, the radii are converted to length of side with equal pore volume.

Fig. 3 Microstructure characteristics of OPA: (a) cumulative pore volume fraction versus pore size determined by FIB-nt and N₂ adsorption analyses;²³ (b) cumulative distribution of meso-pore sizes; (c) cumulative distribution of micro-pore sizes with assumed probability density from ref. 18; (d) cumulative distribution of carbonate grain sizes.²⁴

Regarding the structure of the solid phase, Keller et al.²⁴ reported 18 vol% of non-porous carbonates and 17 vol% of non-porous quartz in the sample. The sizes of carbonates ranged from 100 to 300 nm, while the size distribution of quartz was not measured. The reported data of carbonate were converted to the cumulative probability of its grain sizes shown in Fig. 3(d) for model construction. The same cumulative probability is assumed for the non-porous quartz, due to the lack of more accurate data, and used to generate particle sizes of both types in a way similar to pore size generation.

3 Transport modelling and data

3.1. Diffusion equation

The driving force for local diffusion is the concentration gradient between connected nodes with resistance provided by the pores. In reality, pores have different shapes. However, the experimental data of cumulative pore fraction from Keller et al.²³ substituted the irregular shape of pores by equivalent circular shape. To reflect realistic situation, circularity dependent on shape factor is introduced to consider the effect of pore shapes when all pores are assumed to be circular.³⁸ In our work, the pores are considered to be elongated conduits with different pore shapes. The diffusive flow, F_ij [kg s⁻¹], through a pore is described by the Fick's first lawwhere D_ij [m² s⁻¹] is the pore diffusivity, which is calculated in the same way as in Xiong et al.,^1,39A_ij [m²] is the pore cross-section area, l_ij [m] is the pore length, c_i − c_j [kg m⁻³] is the concentration difference between nodes, 4πG_ij [—] is the circularity.

To verify the pore network model, two diffusion species are considered in this work. One is HTO, a neutral solute with molecular size of r_A = 0.1 nm.¹⁸ The other is a neutral U(VI) complex, Ca₂UO₂(CO₃)₃, which dominates the U(VI) speciation in OPA pore water at room temperature to nearly 100%.⁴⁰ The molecular size of this complex is r_A = 0.524 nm.⁴¹ These two species are selected for comparison with experimental data on their diffusivity in OPA. The free molecular diffusion coefficients calculated with the Einstein–Stokes equation at room temperature, are D₀ = 2.24 × 10⁻⁹ m² s⁻¹ for HTO⁴² and D₀ = 4.66 × 10⁻¹⁰ m² s⁻¹ for U(VI),⁴³ respectively.

3.2. Flow equation

We assume saturated flow through the network. The flow field in the network is described by imposing two different pressures on two opposing boundaries of the network. The remaining boundaries of the network do not permit flow. For a pore with area A_ij [m²] and shape factor G_ij, connecting nodes i and j, the volume flow rate of diffusing species, q_ij [m³ s⁻¹], is described by the Hagen–Poiseuille equationwhere m [—] is a pore-shape coefficient (for a circular, an equilateral and a square tube, m is 0.5, 0.6, and 0.5623, respectively⁴⁴), μ [Pa s] is the dynamic viscosity, p_i and p_j [Pa] are pressures at pore i and j, respectively. Eqn (5) is considered to be appropriate for describing flow in pores⁴⁵ and is valid for laminar flow in a wide range of Reynolds numbers. For incompressible, steady state flow, the sum of pore discharges at node must be zerowhere z_i [—] is the number of pores connected to node i. Then the average pore water velocity through the pore network, [small nu, Greek, macron] [m s⁻¹], can be determined aswhere Q [m³ s⁻¹] is the total discharge through the network; L [m] is the network length in the flow direction; V [m³] is the volume of the network.

3.3. Mass transport equation

The transient dynamics of a solute subject to sorption and transport are described by the following equations:where t [s] is time, c [kg m⁻³] is the concentration of the species in the system, u [m s⁻¹] is the fluid velocity vector, D [m² s⁻¹] is the diffusion coefficient, s is the adsorbed mass onto the solid phase [kg m⁻²], and L is the throat length [m].

A temporal and spatial discretization of above equations is used to describe mass transport through a pore

Substituting eqn (4) and (5) into eqn (12), the mass flow of species in a throat iswhere q_ij [m³ s⁻¹] is the volumetric flow within the throat, c_ij [kg m⁻³] is the average concentration in the pore, s_ij [kg m⁻²] is mass adsorbed per unit area of the pore wall, S_ij [m²] is the surface area, V_ij [m³] is the pore volume.⁴⁶

Linear equilibrium adsorption is described by:

s_ij = k_a,ijc_ij, (14)

where k_a,ij [m] is the distribution coefficient normalized to the surface area (k_d/S_ij).⁴⁷

Then eqn (13) becomes

where β = S_ij/V_ij [m⁻¹].

For nonlinear equilibrium adsorption the Freundlich isotherm⁴⁸ is used:

s_ij = ac^b_ij, (16)

leading towhere a [m^3b−2 (kg^b−1)⁻¹] and b [—] (b > 1) are coefficients for sorption.

3.4. Network transport analysis

The system illustrated in Fig. 2(b) is a 3D mathematical graph – sites are graph nodes, bonds are graph edges. The graph topological structure (connectivity) is described by the so-called incidence matrix of dimensions E × N, where E is the number of edges (pores) and N is the number of nodes (pore junctions). An element of this matrix, a_en, can be either −1, +1, or 0, if node n is the first, second or not incident of edge e.

The incidence matrix describes both the topological structure of the system, which can be used for connectivity analysis, and the derivative of a discrete function defined on the nodes.⁴⁹ Specifically, a concentration field on the graph nodes, vector C̲ of dimension N, has a gradient given by the matrix product , which is a discrete field on the edges, vector of dimension E. As importantly, the transpose of the incidence matrix describes the derivative of a discrete function defined on the edges. Specifically, a mass flow field on the graph edges, a vector J̲ of dimension E, has a gradient given by , which is a discrete field on the nodes, a vector of dimension N.

The components of the discrete field of edge flows, J̲, are calculated by eqn (15) or (17), i.e. each component is the edge gradient scaled by an edge weight given by W_e = (q_ij + D_ijA_ij/L_ij)/(1 + k_aβ) or W_e = (q_ij + D_ijA_ij/L_ij)/(1 + ac^b−1_ijβ) for linear and nonlinear sorption at pore scale, respectively. Clearly, the edge weight depends on the assigned pore radius, velocity in the pore, selected length scale, radius of diffusing species, and sorption coefficient. For convenience the edge weights are arranged in a diagonal matrix, , of dimensions E × E, where the element in row and column e is the weight of edge e. Thus, the discrete flow field is given by .

The mass transport through the graph is governed by mass conservation at nodes, i.e. the rate of change of diffusing species mass at the nodes equals the resultant mass flows through the nodes. The latter is the algebraic sum of the flows through adjacent edges, which is the divergence of the flow field given by . Thus, the process of mass transport through a weighted graph is formulated mathematically by

where V_n is the node volume, which can be taken as cell volume, t is the time. This must be supplemented by initial and boundary conditions, similarly to continuum problems, but now defined on a discrete set of nodes. It is sufficient to point out that the matrix system of eqn (18), supplemented by initial concentrations at all nodes and prescribed concentrations or fluxes at boundary nodes, is solved numerically to find nodal concentrations and boundary fluxes at any future time. The simulations terminate when steady-state is reached, i.e. when concentrations at all nodes are constant.

The formulation of flow via discrete analysis on graphs is highly beneficial for incorporating system evolution compared to traditional finite element or finite difference formulations. Conceptually, it shows the intrinsic relation between the pore space topological structure and transport, both described via the incidence matrix. Practically, it differentiates clearly between effects of connectivity, captured by the incidence matrix , and geometry and physics, captured by the edge weight matrix . This facilitates the incorporation of potential pore space changing mechanisms in a computationally effective way. For example, local pore geometrical changes arising from corrosion, sorption, bacterial film growth, deformation, etc., will affect corresponding coefficients of only but leave constant. Material failures, such as micro-cracks, during evolution present a more complex case. These change the connectivity, i.e. they alter , and introduce geometry and physics for the new edges, i.e. modify .

Immutable pore systems have been analysed for the calculation of macroscopic diffusivity. Comparing the calculated results with experimental diffusion coefficients reported in the literature, the realism of the model construction in terms of initial pore space geometry and topology will be verified. In this work, the effects of pore shapes, sorption coefficient, and velocities are analysed based on previous work.

3.5. Diffusion experiments

OPA bore core samples from the Mont Terri rock laboratory in Switzerland were used in all diffusion experiments applying the experimental set-up described in Van Loon et al.⁵⁰ and shown in Fig. 4. The samples (thickness: 11 mm, diameter: 25.5 mm) were placed in stainless steel diffusion cells between two stainless steel filter plates (316L, pore diameter: 0.01 mm; MOTT industrial division, Farmington, USA). Synthetic OPA pore water⁵¹ was used as background electrolyte in all diffusion experiments. Each diffusion cell was coupled with a peristaltic pump (mod. Ecoline, Ismatec, IDEX Health & Science, Glattbrugg, Switzerland) and a source and receiving reservoir filled with 200 mL and 20 mL synthetic OPA pore water, respectively. The solutions were circulating through the end plates of the cells in order to saturate the samples. The saturation time amounted to three weeks. Subsequently, the solutions were replaced by fresh ones, whereby the source reservoir contained the tracer and thus, tracer diffusion perpendicular to the clay bedding was started.

Fig. 4 Schematic depiction of the experimental set-up used in the diffusion experiments (based on Van Loon et al.⁵⁰).

At first, HTO through- and out-diffusion experiments were performed according to the procedure described by Van Loon et al.⁵⁰ in order to determine values for the total porosity (θ) of the clay samples (c₀(HTO) = 1000 Bq mL⁻¹). After that, the U(VI) in-diffusion was studied (c₀(²³³U(VI)) = 1 × 10⁻⁶ mol L⁻¹). After 3 months, the diffusion experiments were stopped and the clay samples were removed from the cells. With the help of the abrasive peeling technique,⁵² U(VI) diffusion profiles were determined. The peeled layers were extracted for U(VI) content by 1 mol L⁻¹ HNO₃ (p.a., Roth, Karlsruhe, Germany). The tracer concentrations in the extracts were determined by liquid scintillation counting (mod. TriCarb 3100 TR, Perkin Elmer, Freiburg, Germany). All experimental results were evaluated using the commercial software COMSOL Multiphysics: finite-element software package 3.5a (2008) and 4.2a (2011). A detailed description of all diffusion experiments can be found in Joseph et al.^25,27

4 Results and discussion

4.1. Boundary value problems on graphs

With respect to a coordinate system (X₁, X₂, X₃) normal to the square faces of the unit cell (see Fig. 1), a pore network skeleton within the boxed region (0 ≤ X₁ ≤ N·S₁, 0 ≤ X₂ ≤ N·S₂, 0 ≤ X₃ ≤ N·S₃) was used. Here, S₁, S₂, S₃, are the cell sizes in the three coordinate directions and N is the number of cells in each coordinate direction. Pore networks were built on skeletons with increasing N using the process from Sections 2.1–2.4. Statistically representative distribution of pore sizes was obtained with relatively small number of cells, e.g., it was sufficient that N > 5 for the distributed pores to follow the experimental pore size distribution closely. Further, the effect of the random spatial allocation of pores on results was analysed using 10 different realisations for each N. It was found that for N = 20, the variation of results from the average was reduced to under 10%. This was judged to be an acceptable accuracy and the diffusivity results reported hereafter were obtained as the average values from 10 realisations on a skeleton with N = 20. This skeleton contains nodes = 110 860 nodes, i.e. potential pore junctions, and edges = 533 520 bonds, i.e. potential pore locations.

Zero concentration of diffusing species is taken as initial conditions in all nodes. The boundary conditions reflect a particular experimental setup, where concentrations c₀ and c₁, are prescribed on two opposite boundaries, while the remaining four boundaries do not permit flux. The selection of the two boundaries depends on the macroscopic transport properties being analysed. Specifically, the boundary conditions used to calculate the macroscopic transport properties parallel to bedding direction (1, 0, 0), D₁, are: prescribed concentration c₀, in all nodes on plane X₁ = 0; prescribed concentration c₁ in all nodes on plane X₁ = 20S₁; zero flux through all nodes on planes X₂ = 0, X₂ = 20S₂, X₃ = 0, X₃ = 20S₃. For calculating the macroscopic transport properties perpendicular to the bedding direction, say (0, 1, 0), D₂, the boundary conditions are: prescribed concentration c₀ in all nodes on plane X₂ = 0; prescribed concentration c₁ in all nodes on plane X₂ = 20S₂; zero flux through all nodes on planes X₁ = 0, X₁ = 20S₁, X₃ = 0, X₃ = 20S₃.

One example of diffusion in the bedding direction for a system with unequal length parameters and circular cross sectional pores is shown in Fig. 5 by the solute concentration profiles at select times. The prescribed concentration c₀ on plane X₁ is 1 mol m⁻³ and the prescribed concentration c₁ on plane X₀ is 0. These are snapshots taken from the whole history obtained by numerical simulation based on eqn (4). The profiles are shown on the cellular assembly, rather than on the graph representing the pore network, for clarity.

Fig. 5 Concentration profiles of HTO [mol m⁻³] at time t = 0.1 s (a), t = 0.5 s (b), and t = 1 s (c). The pores are assumed with circular cross section.

4.2. Effective diffusivity

The constructed pore networks exhibit macroscopic tortuosity (path lengthening), introduced by the selection of transport pathways along the interfaces between solid phases regions. This tortuosity depends on the material texture, represented here by the ratios of the cell length parameters in three perpendicular directions. Experimental results show that tortuosity is smaller in the bedding direction.²⁶ Therefore, we investigate the effect of larger cell length in the bedding direction and smaller cell lengths in the directions perpendicular to bedding. The results of different cell length ratios have been discussed in a previous paper.³⁹ It was shown that the calculated effective diffusivity considering both carbonates and quartz as particles is closest to the experimental value when the ratio of the cell length parameters is: S₁/S₂ = 2 and the out-of-bedding directions are not differentiated, i.e. S₂ = S₃ is assumed. In this work, the same ratio of the cell length parameters and volume percentage of particles are adopted. For cellular assembly of this given ratio, 10 realisations of pore spatial distributions were analysed to obtain the transport in the bedding, S₁, and out-of-bedding, S₂, directions. The results reported are the averaged values of these analyses.

The diffusion of tracer molecules in clay has shown that mainly neutral and cationic species diffuse in the interlayer water, diffuse layer water and the free pore water.⁵³ Experiments reported that the accessible porosity for HTO^25,26,37,54 agrees very well with total porosity in OPA,²² which implies that almost all the pores in OPA are available. The total porosity has been used to simulate the effective diffusion coefficients of HTO in this work. In Fig. 6, the effective diffusivities of HTO and U(VI) parallel and perpendicular to the bedding, D₁ and D₂, are shown to increase as the average shape factor G increases. The calculated effective diffusivity of HTO is in the following ranges: D₁ = 5.41 × 10⁻¹¹ ∼ 17.10 × 10⁻¹¹ m² s⁻¹; D₂ = 1.58 × 10⁻¹¹ ∼ 5.03 × 10⁻¹¹ m² s⁻¹ (see Fig. 6(a)). Reported experimentally obtained values for HTO diffusion in OPA shown in Fig. 6(a) are in the range of D₁ = (5.4 ± 0.4) × 10⁻¹¹ m² s^−1 26 and D₂ = (1.48 ± 0.03) × 10⁻¹¹ m² s⁻¹,²⁵ (1.26 ± 0.09) × 10⁻¹¹ m² s⁻¹,²⁵ (1.23 ± 0.06) × 10⁻¹¹ m² s⁻¹,²⁶ (1.67 ± 0.12) × 10⁻¹¹ m² s⁻¹,²⁷ (1.33 ± 0.11) × 10⁻¹¹ m² s^−1 27 determined for different bore core samples. In OPA, the accessible porosity is different for the various species resulting in different effective diffusivities. The increased diffusion of cations is explained either by interlayer or by diffusion in the electrical double layer^37,55 that surrounds the negatively charged clay surface and contains an excess of cations.^54,56–61 The diminished diffusion of anions is due to their limited access to the interlayer space. This is caused by the overlay of double layers. As a result the electric potential becomes considerably large.⁶²

Fig. 6 Calculated effective diffusivity of HTO (a) and U(VI) (b) in OPA as a function of shape factor and accounting for carbonates and quartz as non-porous particles with combined volume fraction of 0.35. Experimentally determined ranges of OPA effective diffusion coefficients are shown for comparison.

If in the model the diffusing U(VI) species is neutral, then the same initial accessible porosity as for HTO can be assumed. The simulated effective diffusivity of U(VI) in the bedding direction, D₁, (cf.Fig. 6(b)) is in the range of 8.04 × 10⁻¹² ∼ 2.56 × 10⁻¹¹ m² s⁻¹, and in the perpendicular direction, D₂, is from 1.98 × 10⁻¹² to 7.26 × 10⁻¹² m² s⁻¹. Experimental data for U(VI) in the perpendicular direction are 1.5 × 10⁻¹² ∼ 2.3 × 10⁻¹² m² s⁻¹.^25,27 The simulated results for HTO and U(VI) are in good agreement with experimental investigations when average shape factor, G, is equal to 0.025 (cf.Fig. 6). In OPA, three types of water accessible for diffusion can be distinguished. The solute, for instance a U(VI) species, enters freely the larger pores, since it resides in an electrically neutral solution, the so-called free pore-water. On the net negatively charged clay surface, the water turns into a diffuse double layer. The charge deficiency can be balanced by exchangeable cations, anions are repulsed from the surface. The third type is interlayer water in the smectite fraction of OPA. It is considered to contain exchangeable cations only, which behave chemically and physically comparable to the surface complexed cations. The types of water or pore space accessible for diffusion depends on the diffusing species.^22,37,63,64 Specifically, diffusion of cations is enhanced, and that of anions diminished, relative to a neutral tracer such as HTO.³⁷ The here discussed effective diffusivity excludes sorption effects of the U(VI) species, Ca₂UO₂(CO₃)₃(aq), on OPA during their migration, only the pure transport process is considered. Since the same approach was used for modelling HTO and U(VI) diffusion except that the radius of the diffusing species was changed, the observed difference of one order of magnitude in D₂ for HTO and U(VI) can be attributed solely to size exclusion effects and not to charge effects of attracted or repulsed charged U(VI) species. It has been demonstrated in previous works,^54,65 that the amount of electrical double layer water is approximately equal to half of the porosity. This means that anions cannot access almost half of the porosity, manifested in anion exclusion effect on diffusivity. Hence, the effective diffusivity of UO₂(CO₃)₃⁴⁻ and U(VI) species should be different despite of their similar molecular sizes. Therefore, it can be concluded that the neutral Ca₂UO₂(CO₃)₃(aq) complex is not only the dominant U(VI) species in the OPA pore water solution in the source reservoir but also the main diffusing U(VI) species through OPA.

Overall, the model in this work can predict diffusion in OPA very well. The difference between the computational and the experimental results could be partially due to a difference between the microstructure characteristics of OPA obtained by Keller et al.^23,24 used for model construction, and the clays used for experimental measurement of effective diffusion coefficients by Van Loon et al.²⁶ and Joseph et al.^25,27

4.3. Apparent diffusivity

To investigate the effects of the shape factor on sorption, apparent diffusivities of U(VI) with different average shape factors are analysed and presented in Fig. 7(a). It is obvious that the increase of the average shape factor in the system yields an increase of apparent diffusivities of U(VI) in both directions investigated. Because the shape factor is the ratio between cross section area and perimeter, smaller shape factors lead to smaller ratios between pore volume and surface area. The larger the surface area the more species adsorb onto the pore wall. The variation of U(VI) apparent diffusivities is bigger than the growth of U(VI) effective diffusivities when the average shape factor increases from 0.025 to 0.08. In addition, the change of apparent diffusivities in the bedding direction is larger than apparent diffusivities perpendicular to the bedding plane. This is due to the fact that elongated pores are distributed in the bedding direction. The shape factor only affects cross sectional shape.

Fig. 7 Simulated apparent diffusivity of U(VI) in OPA. (a) The dependence on shape factor when distribution coefficient, k_a, equals 200 nm; (b) as a function of k_a, when G = 0.025, experimentally determined ranges of D_a values are shown for comparison.

Because effective diffusivities of U(VI) agree well with experimental values when the average shape factor, G, equals to 0.025, this value of G is adopted to study sorption effects of U(VI) in OPA. We assumed that linear sorption occurs at pore-scale according to experimental investigation.³⁵ Apparent diffusivities of U(VI) species are shown in Fig. 7(b) as a function of distribution coefficients, k_a. It is clear that the increase of k_a yields reductions of the apparent diffusivity, D_1a and D_2a. Because the adsorbed mass is dependent on the value of k_a, more species will be adsorbed when k_a increases. As a consequence, migration of the U(VI) species retards because of their stronger interaction with the OPA pore walls. In addition, since more species are adsorbed onto the pore-wall, the cross section of pores allowed for diffusion becomes smaller and thus the accessible pore space will be reduced. Notably, the decrease of the apparent diffusivity becomes slowly as the adsorption coefficient increases. This is an illustration of the non-linear effect of the sorption in the system.

The calculated results of apparent diffusivities are in the following ranges: D_1a = 8.35 × 10⁻¹⁴ ∼ 2.05 × 10⁻¹³ m² s⁻¹; D_2a = 1.98 × 10⁻¹⁴ ∼ 4.87 × 10⁻¹⁴ m² s⁻¹, see Fig. 7(b), when k_a decreases from 250 nm to 100 nm. Reported experimentally obtained values for U(VI) complexes diffusion in OPA are D_1a = 9.5 × 10⁻¹⁴ m² s⁻¹,²⁸D_2a = 2.8 × 10⁻¹⁴ ∼ 3.4 × 10⁻¹⁴ m² s⁻¹,^25,27 2.8 × 10⁻¹⁴ m² s⁻¹.²⁸ These ranges are given in Fig. 7(b). Direct comparison shows that the simulation and the experimental results are in reasonable agreement when k_a is about 200 nm. The experimental distribution coefficients k_d measured by experiments were parallel to the bedding: 0.0142 m³ kg⁻¹,²⁸ perpendicular to the bedding: 0.025 ± 0.003 m³ kg⁻¹,²⁵ 0.038 ± 0.005 m³ kg⁻¹,²⁷ 0.0138 m³ kg⁻¹.²⁸ Since k_a is the distribution coefficient, k_d, normalised by surface area, and the surface area obtained by ethylene glycol monoethyl ether for OPA varies from 80 to 107 m² g⁻¹ (the average value is 94 m² g⁻¹),²² the experimental k_a values range from 128 nm to 350 nm. This is very close to the simulated value of k_a. The sorption of U(VI) species is therefore linear.

The pore-scale Freundlich isotherm sorption, s = ac^b (cf.eqn (16)), is also studied in this work. Since this sorption equation is dependent on the pore-local concentration variation of solutes, the concentration applied at the boundaries must be clarified. This is different for the analysis of effective diffusion and linear sorption where the results are independent of the prescribed concentration difference on the boundaries. We apply C₀ = 0.05 mol m⁻³ at X₀ = 0, C₁ = 0 at X₁ = 20S₁, and no flux at other boundaries. The effects of parameters a and b on the apparent diffusivities are shown in Fig. 8. In Fig. 8(a), the increase of a yields a reduction of the apparent diffusivities of U(VI) in both directions when b = 1.5. The changes are similar with the situation of linear sorption, s = k_ac. In Fig. 8(b), the apparent diffusivities of U(VI) in both directions increase as b increase. The increase of apparent diffusivities becomes slowly with increasing b. This is due to the fact that the adsorbed mass becomes less as the parameter of b increases. The maximum values of apparent diffusivities will be very close to their effective diffusivities when b increases.

Fig. 8 Simulated apparent diffusivity of U(VI) in OPA by assuming a nonlinear sorption. Data are presented as a function of the Freundlich coefficient a (a) and the Freundlich exponent b (b). The shape factor equals to 0.025.

4.4. Velocity effect on transport property

Although hydraulic conductivity is very slow in OPA, it may be crucial for mass transport in other porous media. To illustrate the mechanism of diffusion–sorption–convection in porous media, we still use OPA as an example. The effect of velocity on transport of uranium species is considered. The shape factor is 0.025 and the distribution coefficient k_a is 200 nm. The evolution of uranium concentration at the outlet surface with time for different velocities is described in Fig. 9. Obviously, convection accelerates transport of species through the pore network model. The larger the velocity through the pore network model, the faster the species arrive the outlet surface. There is not much variation in the time evolution of the outlet concentration when velocity increases from 0 to 0.186 mm s⁻¹. This is in agreement with the fact that not convection but molecular diffusion is the main transport mechanism for radionuclides in OPA.

Fig. 9 Time history of uranium concentration at the outlet surface for different velocities, the distribution coefficient k_a is 200 nm, the pore shape factor is 0.025.

5 Conclusions

A discrete model for mass transport through porous media has been proposed to balance incomplete knowledge of geometrical and topological pore system characteristics.

The effect of clay's texture on the transport properties has been investigated by analysing the transport of HTO and U(VI) through OPA. The anisotropy and heterogeneity of the clay can be simulated via different length parameters and percentage of pores in different directions. Different pore shapes are also taken into account in this model, which is very important for multi-phase transport. Calculated diffusion coefficients are within the ranges of experimental data. The results suggest that U(VI) diffuses as neutral species, Ca₂UO₂(CO₃)₃(aq), through OPA. Because only species in the source reservoir solution are known in experiments, the methodology proposed in this work could clarify the kind of U(VI) species (anion, cation or neutral species) diffusing in OPA. In addition, the idea has been extended to account for sorption and convection. Physical sorption is simulated in each pore, which is more accurate and realistic compared to previous coarse grained mathematical formula by changing pore radii.¹

Although this work has used a particular material to illustrate the methodology and its application, the model is suitable for modelling and simulation of other micro- and meso-porous materials. The methodology proposed in this work can be used to study sorption of species that may not be easily studied by classical in-diffusion or through-diffusion methods, for example, when the species is adsorbed onto laboratory materials (e.g., cells, vessels, and tubing).

The proposed model can be potentially used for coupling with deformation and fracture of porous media, where the pore system and the solid phase, represented as dual graphs, will interact effectively. Models mimicking the coupling are possible and are subject of ongoing work. For example, our immediate plans are to study the effect of micro-cracking due to thermo-mechanical loads on the diffusion with application to long-term behaviour of OPA in a nuclear wastes repository. The goal with the concept development is to have a set of discrete models, i.e. PNM accounting for transport mechanism and lattice models of the solid phase developed for analysis of damage evolution via micro-cracking,¹⁹ that evolve concurrently at their own rates and inform each other about changes in their structures. Such a coupled problem could be of great interest in many disciplines, e.g., shale gas extraction, dams, and underground cavities in hydropower engineering. In these applications, the evolution of permeability rather than diffusivity of the medium is required, but the methodology developed in this work would be the same.

Disclaimer

This document was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes (LLNL-JRNL-676844).

Acknowledgements

Q. Xiong acknowledges gratefully the Doctoral Scholarship Award from the President of The University of Manchester. A. Jivkov acknowledges the support of EPSRC via grant EP/K016946/1, “Graphene-based membranes” and from EdF R&D for the Modelling & Simulation Centre. C. Joseph and K. Schmeide thank the German Federal Ministry for Economic Affairs and Energy (BMWi) for financial support (no. 02 E 10971). This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

References

1. Q. Xiong, A. P. Jivkov and J. R. Yates, Microporous Mesoporous Mater., 2014, 185, 51–60.
2. M. J. Blunt, B. Bijeljic, H. Dong, O. Gharbi, S. Iglauer, P. Mostaghimi, A. Paluszny and C. Pentland, Adv. Water Resour., 2013, 51, 197–216.
3. C. Ezeuko, A. Sen, A. Grigoryan and I. Gates, Biotechnol. Bioeng., 2011, 108, 2413–2423.
4. N. C. Marty, C. Tournassat, A. Burnol, E. Giffaut and E. C. Gaucher, J. Hydrol., 2009, 364, 58–72.
5. S. Whitaker, AIChE J., 1967, 13, 420–427.
6. A. F. White and S. L. Brantley, Chem. Geol., 2003, 202, 479–506.
7. J. M. Köhne, S. Schlüter and H.-J. Vogel, Vadose Zone J., 2011, 10, 1082–1096.
8. M. J. Blunt, Curr. Opin. Colloid Interface Sci., 2001, 6, 197–207.
9. M. J. Blunt, M. D. Jackson, M. Piri and P. H. Valvatne, Adv. Water Resour., 2002, 25, 1069–1089.
10. B. Berkowitz and R. P. Ewing, Surv. Geophys., 1998, 19, 23–72.
11. D. Wilkinson and J. F. Willemsen, J. Phys. A: Math. Gen., 1983, 16, 3365.
12. J. Meyers, O. Crosser and A. I. Liapis, J. Chromatogr. A, 2001, 908, 35–47.
13. A. P. Jivkov, C. Hollis, F. Etiese, S. A. McDonald and P. J. Withers, J. Hydrol., 2013, 486, 246–258.
14. S. Gao, J. N. Meegoda and L. Hu, Int. J. Numer. Anal. Methods Geomech., 2012, 36, 1954–1970.
15. R. Al-Raoush and C. Willson, J. Hydrol., 2005, 300, 44–64.
16. H. Dong and M. J. Blunt, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 036307.
17. M. Piri and M. J. Blunt, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 026302.
18. A. P. Jivkov and Q. Xiong, Transp. Porous Media, 2014, 105, 83–104.
19. A. P. Jivkov, M. Gunther and K. P. Travis, Mineral. Mag., 2012, 76, 2969–2974.
20. A. P. Jivkov and J. R. Yates, Int. J. Solids Struct., 2012, 49, 3089–3099.
21. A. P. Jivkov, D. L. Engelberg, R. Stein and M. Petkovski, Eng. Fract. Mech., 2013, 110, 378–395.
22. NAGRA, Technischer Bericht, 2002, 02–03, 250–252.
23. L. M. Keller, L. Holzer, R. Wepf and P. Gasser, Appl. Clay Sci., 2011, 52, 85–95.
24. L. M. Keller, P. Schuetz, R. Erni, M. D. Rossell, F. Lucas, P. Gasser and L. Holzer, Microporous Mesoporous Mater., 2013, 170, 83–94.
25. C. Joseph, L. R. Van Loon, A. Jakob, R. Steudtner, K. Schmeide, S. Sachs and G. Bernhard, Geochim. Cosmochim. Acta, 2013, 109, 74–89.
26. L. R. Van Loon, J. M. Soler, W. Müller and M. H. Bradbury, Environ. Sci. Technol., 2004, 38, 5721–5728.
27. C. Joseph, K. Fritsch, R. Steudtner and K. Schmeide, Workshop des Verbundprojekts Rückhaltung endlagerrelevanter Radionuklide im natürlichen Tongestein und in salinaren Systemen, München, Germany, 2012, vol. 10, pp. 29–30.
28. M. García-Gutiérrez, U. Alonso, T. Missana, J. L. Cormenzana, M. Mingarro, J. Morejón and P. Gil, in Informes Técnicos Ciemat 1168, ed. CIEMAT, Departamento de Medio Ambiente, Madrid, 2009.
29. S. Kumar, S. K. Kurtz, J. R. Banavar and M. Sharma, J. Stat. Phys., 1992, 67, 523–551.
30. X. Wang, Z. Yang, J. R. Yates, A. P. Jivkov and C. Zhang, Constr. Build. Mater., 2015, 75, 35–45.
31. X. Wang, Z. Yang and A. P. Jivkov, Constr. Build. Mater., 2015, 80, 262–272.
32. M. Zhang and A. Jivkov, Constr. Build. Mater., 2014, 66, 731–742.
33. D. Zhou, L. A. Dillard and M. J. Blunt, Transp. Porous Media, 2000, 39, 227–255.
34. G. Mason and N. R. Morrow, J. Colloid Interface Sci., 1991, 141, 262–274.
35. H. R. Wenk, M. Voltolini, M. Mazurek, L. R. Van Loon and A. Vinsot, Clays Clay Miner., 2008, 56, 285–306.
36. M. Mazurek, A. Gautschi, P. Marschall, G. Vigneron, P. Lebon and J. Delay, Physics and Chemistry of the Earth, Parts A/B/C, 2008, 33(suppl. 1), S95–S105.
37. C. A. J. Appelo, L. R. Van Loon and P. Wersin, Geochim. Cosmochim. Acta, 2010, 74, 1201–1219.
38. M. E. Houben, G. Desbois and J. L. Urai, Appl. Clay Sci., 2013, 71, 82–97.
39. Q. Xiong and A. P. Jivkov, Mineral. Mag., 2015 , in press.
40. C. Joseph, K. Schmeide, S. Sachs, V. Brendler, G. Geipel and G. Bernhard, Chem. Geol., 2011, 284, 240–250.
41. J. Bai, C. Liu and W. P. Ball, Environ. Sci. Technol., 2009, 43, 7706–7711.
42. MontTerriProject, Report TR, 2010, 2009–2004.
43. M. Cappelezzo, C. Capellari, S. Pezzin and L. Coelho, J. Chem. Phys., 2007, 126, 224516.
44. T. Patzek and D. Silin, J. Colloid Interface Sci., 2001, 236, 295–304.
45. B. Jacob, Soil Sci., 1972.
46. M. Gharasoo, F. Centler, P. Regnier, H. Harms and M. Thullner, Environ. Model. Softw., 2012, 30, 102–114.
47. T. E. Payne, V. Brendler, M. J. Comarmond and C. Nebelung, J. Environ. Radioact., 2011, 102, 888–895.
48. H. Freundlich, J. Phys. Chem., 1906, 57, 385.
49. L. J. Grady and J. R. Polimeni, Discrete calculus: applied analysis on graphs for computational science, Springer, 2010.
50. L. R. Van Loon, J. Soler and M. Bradbury, J. Contam. Hydrol., 2003, 61, 73–83.
51. F. Pearson, Opalinus clay experimental water: A1 type, Version 980318. Technical Report TM-44-98-07, Paul Scherrer Institute, Villigen Switzerland, 1998.
52. L. R. Van Loon and J. Eikenberg, Appl. Radiat. Isot., 2005, 63, 11–21.
53. I. C. Bourg, A. C. M. Bourg and G. Sposito, J. Contam. Hydrol., 2003, 61, 293–302.
54. C. A. J. Appelo and P. Wersin, Environ. Sci. Technol., 2007, 41, 5002–5007.
55. M. A. Glaus, M. Aertsens, C. A. J. Appelo, T. Kupcik, N. Maes, L. Van Laer and L. R. Van Loon, Geochim. Cosmochim. Acta, 2015, 165, 376–388.
56. J. Lehikoinent, A. Muurinen and M. Valkiainen, MRS Proc., 1999, 556, 663,  DOI:10.1557/PROC-556-663.
57. M. Ochs, B. Lothenbach, H. Wanner, H. Sato and M. Yui, J. Contam. Hydrol., 2001, 47, 283–296.
58. J. Van Schaik and W. Kemper, Soil Sci. Soc. Am. J., 1966, 30, 22–25.
59. P. Leroy and A. Revil, J. Colloid Interface Sci., 2004, 270, 371–380.
60. P. Leroy, A. Revil and D. Coelho, J. Colloid Interface Sci., 2006, 296, 248–255.
61. D. Jougnot, A. Revil and P. Leroy, Geochim. Cosmochim. Acta, 2009, 73, 2712–2726.
62. P. Wersin, L. R. Van Loon, J. Soler, A. Yllera, J. Eikenberg, T. Gimmi, P. Hernan and J.-Y. Boisson, Appl. Clay Sci., 2004, 26, 123–135.
63. M. H. Bradbury and B. Baeyens, J. Contam. Hydrol., 2003, 61, 329–338.
64. C. A. J. Appelo, A. Vinsot, S. Mettler and S. Wechner, J. Contam. Hydrol., 2008, 101, 67–76.
65. Q. Xiong, A. P. Jivkov and S. M. Ahmad, Appl. Clay Sci., 2015 , in press.

---