Investigation of the 3D pore structure of a natural shale-implications for mass transport

The multiscale pore structure of a natural shale is obtained by three distinct imaging means. First, micro-tomography image data are extended to provide the spatial arrangement of the minerals and pores observable with a voxel size of 700 nm (denoted here as the macroscopic scale). Second, FIB/SEM provides a 3D representation of the porous clay matrix on the so-called mesoscopic scale (10-20 nm); a connected pore network, devoid of cracks, is obtained for two samples out of five, while the pore network is connected through cracks for two other samples out of five. Third, the nanometric pore network is characterized with tomographic STEM. Using these experimental pore structure data, permeability calculations are performed by the Lattice Boltzmann Method on the nanoscale, on the mesoscale, and on the combination of the two. Upscaling is finally done (by a finite volume approach) on the larger macroscopic scale. Calculations show that, in the absence of cracks, the contribution of the pore structure at the nanoscale, on the overall permeability, is similar to that of the mesoscale. The impact of the most recent tomographic STEM measurements on the overall transport properties is discussed.


Introduction
In several industrialized countries, long lived radioactive wastes are currently stored in temporary conditions, awaiting a long term repository solution. Similarly to Switzerland and Belgium, in France, the repository is anticipated in a clayey natural medium, located at 500 m depth in the East of the Paris Basin. Extensive research has been undertaken to determine the safety of the site, and particularly the fluid transport properties of the host rock. The latter is a Callovo-Oxfordian claystone (COx), made of 30-60% clay matrix and of non porous minerals. On the macroscopic scale (i.e. the centimeter scale and above), COx porosity is on average 18%+/-4. However, this parameter alone is not sufficient to characterize the percolating pore network, and hence, the fluid transport ability of the COx. Indirect permeability measurements [1] have proven that it is of the order of 10 -21 m 2 for undisturbed claystone, and several orders of magnitude bigger in the Excavation Damage Zone (EDZ) [2]. In comparison, indirect measurements of the pore network (obtained mainly by combining nitrogen adsorption and mercury intrusion porosimetry) provide two peak pore sizes of 20 nm and 3-4 nm, but no information on the morphology of the pore network (pore path lengths, tortuosity, connectivity, percolating or non percolating parts, distribution, location, spatial organization depending on pore size, etc). Multiscale imaging techniques have the ability to provide 3D pore morphology, but the smaller the pore size (and the voxel size), the smaller the sample size. In order to image pores of sizes ranging between a few nm to tens of nm, samples with a width on the order of a few microns are necessary [3,4]. Generally, one single sample does not provide a representative volume of the whole pore network. Several samples are needed to obtain pore network characteristics in a statistical sense.
For the COx, 3D micro-tomography (micro-CT) has been undertaken by [5] with a voxel size of 700nm, and a relative pore volume of 0.5% is imaged over the investigated volume. This is far below the 18% measured by indirect means. However, this scale provides the spatial distribution of the main components of the COx, namely the clay matrix (60% vol.), non porous minerals (i.e. tectosilicates, carbonates and heavy minerals for 39.5% vol.) and 0.5% porosity. Finer imaging techniques have been undertaken for Opalinus clay, which is the geological medium investigated for underground nuclear waste repository in Switzerland [6,7]. These are FIB/SEM (a Focused Ion Beam coupled to a Scanning Electron Microscope), 2D TEM (Transmission Electron Microscopy) and 3D STEM (tomographic Scanning Transmission Electron Microscopy). FIB/SEM and 2D TEM have also been used for the COx, for five different samples with FIB/SEM and about ten samples with TEM [3,4,8]. FIB/SEM provides voxel sizes on the order of 5-10 nm, enabling to image pore sizes down to 10-20 nm, i.e. the first peak of the pore size distribution (PSD) of the COx. 2D TEM provides the smaller peak pore size of 3-4 nm (using pixel sizes of 0.36 to 1.04 nm). However, in our former research [3,4], only simplistic fluid transport prediction has been undertaken, by using the 1D analytical Katz-Thompson relationship between permeability, porosity, tortuosity and peak pore size [9,10]. The nanometric pores have only been partially characterized and quantified up to now, using 2D TEM and a numerical 3D image reconstruction algorithm [8].
In this context, this particular research aims at providing first results of the 3D characterization of the nanometric pore network of the natural shale, using tomographic STEM. The final aim is to fit the numerical 3D images of the nanometric pore network, and to provide a more reliable prediction of fluid transport through the shale, by scale change techniques, as in [8]. This macroscopic property needs to be fully characterized and understood, in order to analyze the radionuclide retention mechanisms. The latter are critical properties for long term radioactive waste storage.

Imaging techniques and analysis
The Slice and View procedure [3] is used to provide FIB/SEM stacks of between 200 and 310 bidimensional grayscale images (16 bits, in a xy-plane) and regularly spaced along the z-axis by a constant distance of 10 or 20 nm. This corresponds to volumes ranging between 27.8 and 146.7 μm 3 ( Table 1). 3D STEM images are recorded from a finely polished claystone sample (down to about 30 microns), followed by ion beam milling performed under a liquid nitrogen environment, in order to avoid sample damage due to the ion beam heat. The selection of the observation zone is done purposely inside the clay matrix (neither in the non porous grains nor in their interfaces with the clay matrix). As detailed in [12], each bidimensional grayscale image is recorded in STEM mode (Scanning Transmission Electron Microscopy), by using a HAADF (High-Angle Annular Dark Field) detector. The intensity in the HAADF image depends on the density, thickness and atomic number Z [13]. The probe semi-convergence angle is taken at 16 mrad and the collection angles for the HAADF detector are between 50 and 100 mrad. The beam size is on the order of 500 pm. The 3D information is retrieved from a tilted series of 2D HAADF images acquired at different angles of the sample (between -60° and 60°), at a constant step of 2°, in continuous mode acquisition. All projection images are aligned for drifts in x and y with a 1-4 pixels precision with Digital Micrograph (Gatan software) and ImageJ, to produce a faithful representation of the original object. The reconstructions have been performed using the simultaneous iterative reconstruction technique (SIRT) of the tomoJ plug-in [14]. The resulting reconstruction consists of grayscale voxel intensities exhibiting the 3D distribution of the sample density, with the darkest parts corresponding to pores ( Fig. 1 top left).
For both image types, the 17 segmentation algorithms available in the free image analysis software ImageJ are tested and the most realistic one is selected (by visual inspection and as belonging to a group of at least two other algorithms providing close porosity for the stack, by less than 1%). The Yen algorithm [15,16] is retained for all FIB/SEM images. For 3D STEM images, results are presented herein for all segmentation algorithms, with a peculiar focus on the RenyiEntropy algorithm [17] (leading to the most frequent, i.e. median, porosity value) and on the IsoData algorithm [18] (leading to a significantly greater porosity value, but a still plausible one). Fig. 1 (top right) shows an example of segmentation result for a STEM image, using the RenyiEntropy algorithm. Prior to numerical calculations, a 3D pore structure is reconstructed and may be quantified, e.g. for its relative volume (Fig. 1 bottom). For FIB/SEM, Fig. 2 shows the 3D rendering of the pore volume (with the Volume Viewer tool in Fiji) for a typical sample (out of five).

Permeability derived from the Stokes equation
Whenever the pore structure is available, the corresponding permeability can be derived by solving the Stokes equation. This procedure is applied to the nanoscale and mesoscale structures, and to their superposition. This equation is usually expressed as −∇ + ∇ = 0 , ∇ • = 0 in Vp (1) where p and v denote pressure and velocity, respectively; is the fluid viscosity, and Vp the pore volume. This system should be supplemented by boundary conditions. Assuming that the solid matrix is impermeable, the fluid velocity should vanish on the solid surface Sp, which limits the pore space v=0 on Sp (2) In addition, spatially periodic boundary conditions are applied at the lateral boundaries of the sample. The driving force is a macroscopic pressure gradient ∇ imposed on the sample.
The seepage velocity is the average of the local velocity over the pore volume and it is proportional to ∇ The components of the permeability tensor K are denoted by the names of the corresponding axes, for instance the diagonal component along the x-axis is noted Kxx. The numerical resolution of the Stokes equation can be done in different ways, namely a classical Computational Fluid Dynamic (CFD) technique applied to cubic voxels [19], a Lattice Boltzmann Model [20] and the finite volume technique applied to tetrahedra [21].
In this contribution, the second technique is applied since the corresponding code is parallel, and thus significantly faster for the sample volumes considered (in tens of million voxels). According to the standard terminology, it is a D3Q19 code working in three dimensions with nineteen velocities. For a better precision, the model is TRT, i.e., with Two Relaxation Times. The classical bounce-back condition is used at the solid interface. More details can be found in [20]. The LBM code has been compared in the past to the two other techniques just mentioned.
Since the algorithm assumes a spatially periodic structure, each sample is completed by its mirror image along the direction under consideration. This detail will be made clear in each calculation.

Permeability on the macroscale
On the macroscopic scale (macroscale), the medium is not any more composed of solid and pores. Rather, it is a continuum with a locally variable permeability KM(x) yielding the local macroscopic permeability. Then, the local seepage velocity is related to the local pressure gradient by Darcy's law Since for an incompressible fluid, the seepage velocity satisfies the continuity equation, one has for a constant viscosity ∇ • ∇ = 0 (5) Usually, the total block of volume V is submitted to a macroscopic pressure gradient ∇ . The resulting macroscopic seepage velocity is defined as follows, and it is related to ∇ by the macroscopic permeability KMA The elliptic equation (6) is discretized by the so-called box integration method [22] and the resulting linear system is solved by a conjugate gradient method, as detailed in [23].

Image analysis results
Results of FIB/SEM image analysis are presented in [4]. For 3D STEM, the analysis is performed on the first data series available (Table 1). After segmentation by each of the 17 standard thresholding algorithms available in ImageJ, porosity is calculated and plotted (Fig. 3). The median porosity value is of 14.8%. It varies by +/-1.4% when only the four closest algorithms to this median are considered (these are MaxEntropy, RenyiEntropy, Yen and Triangle). The RenyiEntropy algorithm provides the closest value to this median, with a porosity of 14.7%. Porosity obtained with the IsoData algorithm, which corresponds to the biggest porosity value remaining in a plausible range, is of 28.9%. Earlier 2D TEM experiments had provided consistent porosity values of 10.2-25.2% [3]. However, these values are smaller than the expected porosity of the clay matrix at this scale, which is around 35-40% (using that clay porosity is macroscopic porosity divided by clay content) [5,24]. This means that a non negligible part of porosity is not assessed when focusing on the clay platelet arrangement (as is done with 3D STEM). The pore space between clay platelets and other minerals (non porous tectosilicates, carbonates and heavy metals) is also significant. In the following, it is assumed that fluid transport is not affected significantly by this remaining porosity. For fluid transport, they are considered dead spaces (i.e. fluid accumulation spaces), because they are usually bigger than the pore spaces between clay platelets (not shown here). The latter are quantified by 3D STEM.

Numerical results
On the pore scale, the Stokes equation is solved with spatially periodic boundary conditions. If numerical samples can be easily generated with such conditions, this is not true for measured samples. When such conditions are directly used for experimental samples, this results in a systematic underestimation of permeability, since the pores on two opposite faces of the experimental samples do not generally coincide. In order to avoid this bias, the experimental samples are completed by their mirror images with respect to the plane which is perpendicular to the axis used to determine the permeability. Calculations are performed on the completed samples. This procedure is illustrated in the following Figs.4 and 5. The original segmented sample displayed in Fig.4a is discretized by 453x468x311 cubic voxels of side 3.91 nm. The total porosity is equal to 29%, and the open porosity to 28% (thresholding with the IsoData algorithm). The diagonal components of the permeability tensor are determined by solving the Stokes equation and by using the LBM algorithm. Their values are: Kxx= 6.5 10 -19 , Kyy= 7.7 10 -19 , Kzz= 2.3 10 -17 m 2 (7) Therefore, one can conclude that the x-and y-directions are equivalent, while the permeability along the zdirection is about two orders of magnitude larger. The sample presents a very strong anisotropy.

Effect of thresholding algorithm on permeability
The very same sample was re-discretized using a different thresholding algorithm (RenyiEntropy instead of IsoData). The size of the cubic voxels is identical (3.91 nm). However, total porosity is now equal to 15%, and open porosity to 13%. The resulting samples are displayed The same anisotropy as before is observed on these data. The x-and y-axes are roughly equivalent, and the permeability along the z-axis is about two orders of magnitude larger.

Discussion
Numerical simulation results show that using the experimental 3D STEM series 1 allow to numerically reproduce a significant anisotropy to fluid transport, with a permeability smaller by about two orders of magnitude along the x-and y-axes, compared to the z axis. This may be attributed to experimental artifacts, because the sample thickness along z is significantly smaller than along x and y.
The choice of the thresholding algorithm, to compute the pore structure from the original grayscale images, also has a significant effect on permeability calculation at the nanometric scale. When the two discretizations are compared (either from using the IsoData or the RenyiEntropy algorithm), permeabilities are roughly two orders of magnitude smaller if the open pore relative volume decreases from 28 to 13%, i.e. if open porosity is divided by a factor of 2.15.
Compared to former research using 2D TEM data numerically reconstructed to provide a 3D pore system [8], the obtained permeability values using actual 3D data for the pore structure are only one order of magnitude bigger than when reconstructing the 3D pore system from 2D TEM images at the nanoscale. In the latter case, permeability along the x and y axis is of 10 -21 to 10 -20 m 2 (porosity of 14%), and it is 10 -20 -10 -18 m 2 with 3D STEM, for very close relative pore volume (of 13%). It is inferred that when upscaling at the mesoscale of tens to hundreds of nm, each individual 3D STEM pore system, obtained at the nanoscale, has a different spatial orientation than the other, providing an apparent tortuosity to the fluid flow.

Conclusion
This contribution relates to multiscale fluid transport through a natural shale, from the nm to hundreds of nm. It has focused on the determination of fluid permeability, in the Darcy's sense, from actual 3D STEM images (voxel size of 3.91 nm). Results show that, compared to former nanoscale pore system, this realistic pore system (given by 3D STEM) is more permeable, by about one to two orders of magnitude. Former nanoscale pore system, combined to mesoscale FIB/SEM and micro-CT data, allowed to reproduce the macroscopic permeability of the shale. These new results mean that each individual sample analysed by 3D STEM provides a realistic pore system, with greater permeability than reconstructed data from 2D TEM imaging. In order to recover the expected macroscopic permeability, this result mean that very probably, each 3D STEM sample has a different spatial orientation from one another, providing an additional apparent tortuosity for fluid transport. Owing to supplementary 3D STEM datasets, this will be reproduced in further research.