P1Freely available online through the publisher-supported open access option.


Opalinus clay has a high sealing capacity and is therefore considered as a viable candidate for hosting high-level radioactive waste. Assessment of the long-term containment function of clays requires understanding and modelling of mass transport through evolving pore systems. Development of pore network models for diffusion, which can be coupled to models for deformation and micro-cracking, is reported. Effects of clay texture and solid phase constitution are calculated and analysed. The results are in the range of, but slightly over-predicting, experimentally measured coefficients of diffusion in different clay directions. Further model improvements require better knowledge of micro-pore tortuosity, which awaits higher resolution experimental techniques.


Opalinus clay (OPA) is considered as a potential backfill material in the disposal of radioactive nuclear wastes (Boult et al., 1998), and as the main candidate in the decontamination and treatment of heavy metal ions (Aytas et al., 2009). The principal mass transport mechanism in OPA is diffusion. Important for waste repository designers and risk assessors is to operate with reliable data for diffusive-reactive transport of ions in OPA over long time intervals, particularly long-lived radionuclides (Nuclear Decommissioning Authority, 2014).

Pore network models (PNM) provide a suitable description for dealing with mutable pore systems. The pore space is approximated by a set of sites and a set of bonds connecting some sites. The sites positions may form a regular lattice (regular PNM with one independent length scale) or reflect known length variations (irregular PNM). The bonds act as mass conduits. Historically, pores were related to sites and pore throats were related to bonds of regular PNM (Wilkinson and Willemsen, 1983; Meyer and Klobes, 1999; Dillard and Blunt, 2000; Meyers et al., 2001). In these cases, sufficiently detailed 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 (Gao et al., 2012; Jivkov et al., 2013). Laboratory or synchrotron X-ray computed tomography (XCT) is one way to obtain such information (Al-Raoush and Willson, 2005; Dong and Blunt, 2009). The tight pore space in OPA, however, is not fully resolvable by XCT or nano-XCT. Specifically, connectivity data is not available, which requires OPA models based on limited microstructure information.

A regular PNM with pores at sites and throats between each pair of neighbouring pores was firstly proposed (Xiong et al., 2014). The connectivity effect was not considered and diffusivity anisotropy emerging from material texture (preferred mineral orientations) could not be captured. An improved model with variable connectivity was proposed later (Jivkov and Xiong, 2014). This was more realistic in terms of pore-space topology, but the experimental pore-space data was insufficient to rigorously establish a network length scale.

The previous two approaches to tackle incomplete pore space information suffer from lack of additional constraint. This cannot be found from within the pore space information and requires the consideration of the solid phase structure, e.g. the shape and size distribution of mineral grains. This work proposes a methodology for incorporating solid phase characteristics that substantially improves the realism of the constructed PNM, both in terms of geometry and topology. The proposal is conceptually different from the previous works and has the added benefit that the constructed PNM can be paired directly to models of the solid phase developed for analysis of damage evolution via micro-cracking (Jivkov et al., 2012; Zhang and Jivkov, 2014; Wang et al., 2015a,b). This facilitates mechanistic investigations of combined mechanical-thermal-chemical-biological effects on diffusivity, which will be discussed further in the paper.


Experimental data

The transport properties of Opalinus clay are controlled by pore sizes between 1 nm to 100 nm (Marschall et al., 2005). Further, OPA displays anisotropic responses to deformation and transport due to preferred orientation (or texture) of clay minerals (Wenk et al., 2008). Experiments indicated anisotropic diffusion of solute species with slow diffusion perpendicular and fast diffusion parallel to the bedding plane. Our goal is to construct a model truthful to the available structural data, which calculates diffusion in line with macroscopic observations, and can be linked to models for solid phase deformation.

Regarding porosity, we use 3D analysis of OPA reported by Keller et al. (2011, 2013). The pore space, resolved by Focused Ion Beam nanotomography (FIB-nt), consisted of a large number of pores located predominantly within the fine-grained clay mineral matrix. These were elongated in the bedding plane and had approximately cylindrical shapes with radii >10 nm, for the purpose of our work these are called meso-pores. Further, they were largely isolated and did not provide a percolating network through the sample. The meso-pores occupied ~1.8 vol.% of the sample, i.e. their volume density (porosity) was θmes = 0.018. N2 adsorption analysis gave a total porosity of 11.5 vol.%, of which 9.7 vol.% were attributed to pores with sizes <10 nm, called hereafter micro-pores. The volume density of micro-pores is thus θmic = 0.097.

A single ‘cumulative pore volume fraction – pore radius’ curve given in Fig 1a shows the measured data of the two pore groups (Keller et al., 2011; Jivkov and Xiong, 2014). For model construction, the experimental distribution, Fig. 1a, is re-evaluated as cumulative probabilities for meso- and micro-pores. These are shown in Fig. 1b and Fig. 1c, respectively.

Keller et al. (2013) also reported results from analysis of the solid phase. The non-porous carbonates analysed by FIB-nt occupied ~18 vol.% of the sample with sizes ranging between 100 nm and 300 nm. For construction of the model, the data reported was used to derive the cumulative probability of carbonate grain sizes shown in Fig. 1d. Keller et al. (2013) reported further of the presence of 17 vol.% of non-porous quartz, the size distribution of which was undetermined.

The cumulative probabilities are used to generate pore or particle sizes, belonging to the respective distribution, via random numbers distributed uniformly. More details on the statistics are given by Jivkov and Xiong (2014).

Pore network construction

The available data provide the basis for a rigorous model construction following four steps: (1) selection of the cellular basis for allocation of grains and pores resulting in complementary (dual) lattices for solid and pore systems; (2) allocation of the grains to cell interiors and calculation of the lattice length scale subject to grain volume density and size distribution; (3) allocation of meso-pores along cell boundaries according to their size distribution, relative porosity and in preferred direction(s); and (4) allocation of micro-pores according to their size distribution and relative porosity in domains and directions not occupied by meso-pores.

Statistical analysis of neighbourhoods of arbitrary features, randomly distributed in 3D space, has shown that the average neighbourhood is nearly equivalent to one space-filling polyhedron, the truncated octahedron (Kumar et al., 1992). This has six square faces normal to three perpendicular directions and eight hexagonal faces normal to octahedral directions. The truncated octahedron was the unit cell of a regular space tessellation, proposed for site-bond modelling of solids (Jivkov and Yates, 2012). Considering the available experimental information about the carbonate particles and pores in OPA, and the need to couple the deformation and the transport models in the future, we propose a new network construction based on the regular tessellation of space by truncated octahedrons.

The material is subdivided into cells, truncated octahedrons, representing the neighbourhoods of carbonate particles in OPA. The particles are associated with cell centres (interiors). This is illustrated in Fig. 2a for cells with equal distances between the three pairs of square faces. The 3D graph formed by connecting neighbouring particles, shown in Fig. 2b, was used for deformation and failure analysis of quasi-brittle media (Jivkov et al., 2012; Zhang and Jivkov, 2014). A more general case is considered here to represent texture. The cellular assembly is described by three length parameters, S1, S2, and S3, measuring the distances between square faces in the three orthogonal directions. All other distances, areas and volumes are expressed by S1, S2, S3, e.g. cell volume is Vc = S1S2S3/2.

In an assembly of Nc cells, each cell is assigned a particle with a different radius, ri, generated from Fig. 1d. The volume of all the particles is required to equal the experimentally measured particle volume fraction, ϕ = 0.18, of the assembly volume Vc·Nc. From this requirement  
The calculation of S1, S2, S3 from equation 1 depends on the selection of their ratios, used to represent texture. A non-textured medium has S1 = S2 = S3.

Considering the experimental observation that the meso-pores have elongated cylindrical shapes, we propose a pore network structure complementing the solid phase. Firstly, a skeleton is formed from sites at centres of cell faces and bonds between neighbouring faces. This is illustrated in Fig. 2c on a single cell with S1 = S2 = S3 for clarity. The bonds represent potential diffusion pathways, i.e. they show positions where elongated cylindrical pores are allowed to reside. As pores reside at bonds, sites represent pore junctions – volume-less containers redirecting mass transport.

The transport skeleton, shown in Fig. 2d, is a mathematical graph, dual to the one in Fig. 2b. This facilitates greatly the algorithmic handling of the two graphs, as their matrix representations are related by elementary operations (Grady and Polimeni, 2010). As a consequence, the dual systems, deformation site-bond model and diffusion site-bond model, can be combined efficiently to investigate coupled problems. This is the subject of ongoing work to be reported in a future publication.

The PNM construction maps experimentally measured the porosity to the skeleton. Firstly, meso-pores are assigned to bonds of Fig. 2d in a preferred direction with radii selected from Fig. 1b. Secondly, micro-pores assumed to be cylindrical with radii selected from Fig. 1c are assigned to bonds not already occupied by meso-pores. Assignments terminate when the corresponding pore volume fractions are attained. More details about meso- and micro-pore allocation can be found in Jivkov and Xiong (2014).


Diffusion in bonds (pores) is driven by the concentration difference between sites (junctions). For a pore with radius Re, connecting sites n and m, the mass flow of diffusing species, Je [kg/s], is described by Fick's first law  
where De [m2/s] is effective pore diffusivity, Ae = πRe2 [m2] is pore cross-section area, Le [m] is pore length, and ΔCnm = CnCm [kg/m3] is concentration difference. The effective pore diffusivity of species of radius r0 is given by (Bryntesson, 2002):  
where D0 is the free molecular diffusion coefficient of the species in the medium filling the pore system.

Specifically for this work, OPA is considered fully saturated with water to reflect repository conditions. The diffusing species is Hydrogen Tritium Oxygen (HTO), a neutral solute, selected for comparison with published experimental data on its diffusivity in OPA. The molecular diffusion coefficient of HTO is D0 = 2.24 × 10−9 m2/s (Wersin, 2010). Pore diffusivities enter the graph representation of PNM detailed by Jivkov and Xiong (2014).

Results and discussion

With respect to coordinate system (X1, X2, X3), the skeleton within the boxed region (0 ≤ X1G·S1, 0 ≤ X2G·S2, 0 ≤ X3G·S3) was used. Here, S1, S2, S3, are the cell sizes and G is the number of cells in each direction. The boundary conditions reflect experimental setup: concentrations are prescribed on two opposite boundaries, while the remaining four boundaries do not permit flux. Effective (system) diffusivity is found from the fluxes at the boundaries with concentrations; details in Jivkov and Xiong (2014). Analyses with 10 PNM realizations on the skeleton have shown that for G = 20 the variation of calculated macroscopic diffusivity is reduced to <10%. Results reported below are obtained with G = 20.

Texture effect on diffusivity

Constructed PNM exhibit macroscopic tortuosity (path lengthening), introduced by the selection of diffusion pathways along interfaces between solid phase regions. This depends on the material texture, represented here by the cell length parameter ratios. Experimental results show that tortuosity is smaller in the bedding direction (Van Loon et al., 2004). Therefore, we investigated the effect of larger cell lengths in the bedding direction (S1) and smaller cell lengths in directions perpendicular to bedding (S2, S3). The out-of-bedding directions are not differentiated, i.e. S2 = S3. Three ratios of the cell length parameters are considered: S1/S2 = 2, 1.5, 1. Results are reported in Fig. 3; each point is the average of 10 random PNM realizations on the corresponding skeleton.

The effective diffusivity in bedding direction, D1, is shown to increase as the ratio S1/S2 increases. Because the cell volume is fixed by the solid phase characteristics, particle size distribution and density, the volume of meso-pores is fixed by their volume density. As the meso-pores are entirely in the bedding, increasing S1 reduces the number of meso-pores in PNM. Sparser meso-pore sets could be anticipated to reduce the effective diffusivity D1. However, there are two factors that counteract such an expectation. The meso-pores in sparser sets are longer as S1 increases. This alone reduces the macroscopic tortuosity. Further, the significance of the micro-pores in the bedding direction is increased, because the smaller number of meso-pores leaves space for allocation of larger numbers of micro-pores. The combined effect is a reduction of macroscopic tortuosity and correspondingly, albeit not large, an increase of diffusivity. Clearly this outcome is affected by the assumption of straight micro-pores. Local tortuosity of micro-pores could lead to changes in diffusivity values and trend observed in Fig. 3.

The effective diffusion coefficient perpendicular to bedding direction, D2, is shown to decrease as the ratio S1/S2 increases. This is in agreement with expectations for two reasons. Firstly, the decrease of S2 alone causes an increase of macroscopic tortuosity in out-of-bedding directions. Secondly, the allocation of a larger number of micro-pores in the bedding direction, as the meso-pore set becomes sparser with increasing S1, leads to a smaller number of micro-pores in out-of-bedding directions. In addition, the out-of-bedding micro-pores become shorter, which reduces the effect of approximating them with straight segments. Overall, the texture effect on both diffusion coefficients predicted by the model is smaller than one order of magnitude.

The calculated diffusion coefficients are in the following ranges: D1 = 7.7 × 10−11 ~9.2 × 10−11 m2/s; D2 = 2.9 × 10−11 ~ 5.7 × 10−11 m2/s, see Fig. 3. Values obtained experimentally and reported for HTO diffusion in Opalinus clay are D1 = 5.5 × 10−11 m2/s (Wersin, 2010); D1 = 3.1 × 10−11 ~ 8.0 × 10−11 m2/s (Van Loon et al., 2004; Wang et al., 2005); D2 = 1.0 × 10−11 ~ 2.3 × 10−11 (Van Loon et al., 2003; Wang et al., 2005; Joseph et al., 2013). These ranges are given in Fig. 3. Direct comparison shows that the simulation and experimental results are in reasonable agreement, with the model somewhat over-predicting the diffusivity of HTO in OPA.

Clearly, 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. (2011) and Keller et al. (2013) and used for model construction, and the clays used for experimental measurement of effective diffusion coefficients (Van Loon et al., 2003; Wang et al., 2005; Joseph et al., 2013). The analyses of the trends in Fig. 3, however, suggest significant micro-pore system control on diffusivity. It is reasonable to assume that the OPA systems used for imaging and for diffusion experiments did not differ significantly in the micro-pore range. This leaves the local tortuosity of micro-pores, which is not presently accounted for in the model, as the major reason for over-prediction of diffusivity. One additional possibility is that the carbonates statistics alone are not sufficient for determining cell lengths. This is investigated below.

Cell volume effect on diffusivity

The OPA contains 17 vol.% quartz, which together with the 18 vol.% carbonates, can be considered as 35 vol.% solid particles; see ‘data’ section. Due to the lack of quartz-specific experimental data, quartz is assumed to follow the size distribution of carbonates. The error with such an approximation is probably acceptable for the diffusion model, as the non-porous phases have a cumulative effect on the cell volume. Models with 35 vol.% particles following the procedure described in the sub-section ‘Pore network construction’ were therefore generated. The increased particle volume density reduced the length parameter S1 to 398, 520 and 631 nm, respectively for the textures given by the ratios, S1/S2 = 1, 1.5, 2. The calculated effective diffusivities are shown in Fig. 4 together with the experimentally obtained ranges.

The trends in the macroscopic diffusivity coefficients due to texture in Fig. 4 are similar to those shown in Fig. 3. This is expected, because the only difference between the two types of systems is the cell volume and the arguments provided for these trends in the previous sub-section remain valid. The effect of the reduced cell volume is seen in the smaller effective diffusion coefficients. Specifically, the following ranges are calculated: D1 = 5.8 × 10−11 ~ 8.4 × 10−11 m2/s; D2 = 2.0 × 10−11 ~ 3.7 × 10−11 m2/s, see Fig. 4. The results are closer to the experimentally measured diffusivities, demonstrating that a more realistic account for the solid phase constitution improved the model performance. This understanding will help the methodology for construction and dimensioning of the coupled deformation and fracture model.

The methodology is particularly suitable for modelling and simulation of micro- and mesoporous materials, such as bentonite and shale rocks. On the one hand, the existing experimental techniques do not allow for complete quantitative analysis of their pore systems. On the other hand, the effect of damage (micro-cracking) on mass transport in such materials is fundamental to the understanding of their retention or release functions.

The effect of the clay's texture on the transport properties has been analysed using diffusion of HTO. Calculated diffusion coefficients are within the ranges of reported experimental data, particularly when clay heterogeneity is simulated via different length parameters in bedding and perpendicular directions. The effect of a solid phase constitution on model length parameters and hence diffusivity has also been analysed. It has been shown, that the additional volume of non-porous inclusions, albeit not resolved as particles, generates models with improved predictions for diffusivity. Models mimicking tortuosity are possible and are the subject of ongoing work. These are expected to provide a reverse-engineering technique for determining an average micro-pore tortuosity.


This work presented (1) a new method for pore network construction, where solid phase characteristics of porous media are considered in the model to partly address the incomplete knowledge of geometrical and topological pore system characteristics; and (2) a framework for coupling the proposed model with models for the deformation and fracture of porous media, where the pore system and the solid phase, represented as dual graphs, will interact effectively.


Jivkov acknowledges the support from BNFL for the Research Centre for Radwaste & Decommissioning and from EdF R&D for the Modelling & Simulation Centre.

The publication of this research has been funded by the European Union's European Atomic Energy Community's (Euratom) Seventh Framework programme FP7 (2007–2013) under grant agreements n°249396, SecIGD, and n°323260, SecIGD2.