Process tomography of diffusion, using PET, to evaluate anisotropy and heterogeneity

Abstract Anisotropy and compositional and structural heterogeneity in clays are causes of considerable deviations from homogeneous diffusion, in particular in terms of direction-dependent transport rates and preferred transport zones. Conventional diffusion experiments, in which the sample is treated as a homogeneous black box in a concentration gradient, are interminable and insensitive to spatial effects. In contrast, tomographic imaging methods are capable of both reducing the amount of observation time required and revealing space-dependent features of the diffusion process. In the present study, positron-emission-tomography (PET) was applied as the most sensitive quantitative spatiotemporal tomographic modality for direct observation of positron-emitting radiotracers in opaque media at reasonable resolution (1 mm) on a laboratory scale (100 mm). Geoscientific applications of PET, or GeoPET, have revealed anisotropic and heterogeneous effects in diffusion experiments that have been conducted on Opalinus clay samples of different sizes, as well as on other rock types. Applying the Comsol Optimization Module to 2D-image sections of the PET tomograms, effective parameter values were derived, thereby quantifying the anisotropic diffusion.

Molecular diffusion is the dominant transport mechanism of dissolved species in clays, complemented by species interactions with the large internal surface area along the propagation pathways. Due to small pore sizes (in the nm range), advective flow processes play only a minor role. Diffusion in clays is based on Fickian diffusion (Altmann et al., 2012) of some species in a concentration gradient, restricted by the pore-space geometry of the material and the sorption of the diffusing species on the portion of the internal surface area that comes in contact with the species.
Van Loon et al. (2004) designed a conceptual model for mudrocks (which consist of layered clay platelets that are oriented preferentially by overburden pressure). The oriented heterogeneity on the sub-micrometer scale causes anisotropy in many properties of clay on the sub-millimetre scale. Heterogeneities on the larger, millimetre scale such as intercalations of silty or sandy laminations or lenses cause anisotropic behaviour on a larger scale. The diffusion pathways in such laminations are less tortuous and thus establish a quicker route to the less penetrable clay-rich zones. Also, the specific surface area, and thus the density of the sorption sites, is smaller in these zones than in the pure clay. Such heterogeneities are present in Opalinus clay samples of drill-core size and distinguishable in the µCT-image of a comparable sample (Fig. 1). The heterogeneities may also be observed on the field scale as homogeneous anisotropic behaviour with modified diffusion parameters compared to values from the smaller scale.
These multi-scale material characteristics cause complex diffusion pathways that are difficult to predict with simple homogeneous models. In the range below the centimetre scale, which is the size of the largest distinguishable features in the µCT image ( Fig. 1), the propagation is not homogeneous but fluctuates spatially; this is due both to the effects of platelet orientation on the microscale and layering on the millimetre scale. Some of the internal surface area is not affected immediately by the diffusing species because zones with large sorption capacity are frequently bypassed through less tortuous pathways. These multi-scale heterogeneity effects modify the macroscopic diffusion parameters, and thus anisotropic diffusion parameters with strongly increased transport in the direction of the largest axis have to be considered.
The second Fickian diffusion equation in its general form is: where c is the concentration and D a [m 2 /s] is the apparent diffusion tensor, r is the position vector and t is time. The apparent diffusion tensor groups together geometrical effects ( porosity, internal surface area, tortuosity) and chemical parameters (chemical-state variables, surface complexation and exchange). The scalar geometrical parameters must be considered as effective quantities, including only the portion that is involved in the diffusion process. Whether diffusive propagation of species is predictable on the basis of a small number of parameters thus appears questionable.
In the simplest form these are equilibrium distribution coefficient, K d [cm 3 /g] from sorption batch experiments, total porosity -ε, bulk density -ρ b [g/cm 3 ] and the scalar diffusion coefficient in water -D aq [m 2 /s]: where the capacity factor, α considers sorption effects and the formation factor F considers pore geometry. The Archie's equation with an empirical Archie coefficient of m ≈ 2 may be applied as a simple approximation, or better, a more sophisticated tortuosity model (Boving & Grathwohl, 2001). It is proposed here that effective porosity ε eff < ε, thus F eff > F and an effective distribution coefficient K d,eff < K d , taking into account the smaller effective surface area, because some of the bulk material is not involved in the process. Assuming spatial independence of D a , it follows that: Approaches which are more dependable than simple assumptions about the effective parameters are based on empirical observations of the diffusional propagation of a tracer in an undisturbed material under controlled laboratory conditions. Examples of such experiments are in-, out-, and through-diffusion experiments, conducted in diffusion cells with periodic fluid sampling in the fluid reservoirs (Van Loon et al., 2004). The sample is regarded as a black box, and the diffusion coefficient is derived from the input and output tracer concentration. In the case of slow propagation, these experiments can be extended by destructive analysis of the tracer distribution in the sample (Van Loon & Eikenberg, 2005). These experiments yield exact data for D a in one direction; they are elaborate, however, and ignore spatiotemporal effects which could be caused by the natural heterogeneity of the sample or by damage caused during preparation. These inconveniences can be reduced by applying non-destructive tomographic methods for observing tracer propagation within the sample because the impact of withdrawal of in situ conditions and preparation damage can be minimized. For quantitative imaging of tracer concentrations, positron emission tomography (PET) is particularly suitable. The socalled PET is known from medical applications as a method suitable for molecular imaging, which plays a prominent role in cancer diagnosis. A number of geoscientific applications exist which demonstrate the feasibility and potential of the method for in-depth analysis of transport processes, but also demonstrate the problems such as low spatial resolution, adaptability to the material in question and cost. Efforts continue to develop further the GeoPET method, with numerous non-medical applications (Richter et al., 2005;Gründig et al., 2007;Kulenkampff et al., 2008).
The emitted positron travels a slow-down distance of ∼1 mm, which limits the maximum spatial resolution of the method, and then is annihilated when in contact with an electron, sending out two antiparallel photons with energy of 511 keV. This energy is at the minimum of the mass absorption coefficients of most materials, but attenuation and scattering of the photon radiation must still be taken into consideration. From coincident events, which are detected in a cylindrical detector, the 3D distribution of the concentration of the tracer nuclide is reconstructed quantitatively.
The most profitable feature of the method is the sensitivity and selectivity for quantitative detection of the PET tracer. The sensitivity for tracer concentrations is greater by orders of magnitude than that achievable with any other tomographic modality. The detection threshold of some 10 Bq per voxel depends on the background level and the spatial activity distribution. This means that the sensitivity for tracer concentrations is in the order of picomoles per litre, depending on the decay constant of the nuclide.
A high-resolution PET scanner (ClearPET by Raytest) (Sempere Roldan et al., 2007), with a voxel size of 1.15 mm × 1.15 mm × 1.15 mm was applied. This type of scanner has the physical resolution limit of ∼1 mm, in contrast to common clinical PET scanners with a resolution of 3 to 5 mm. The voxel size, which is the integration volume for tracer concentrations, is adequate for observing heterogeneities on the millimetre scale.
Users make the most of the adjustable detector diameter of the scanner being employed, by setting a large diameter for the field-of-view (FOV) of 166 mm. The high resolution in combination with the large FOV is achieved at the cost of imaging artefacts. Monte Carlo simulations of PET observations provided reference data for adjusting the parameters for enhanced image reconstruction methods (Zakhnini et al., 2013) and helped significantly to improve image quality.

E X P E R I M E N TA L M E T H O D
A horizontal drill core (diameter = 100 mm, length = 80 mm) of Opalinus clay from Mont Terri was cast in epoxy resin (Fig. 1). A blind hole (diameter = 5 mm, length = 50 mm) was drilled into the core perpendicular to the bedding, to allow filling with synthetic Opalinus (OPA) pore water (Pearson et al., 2003). After equilibration with synthetic OPA-water over 4 months the hole was refilled with synthetic OPAwater that was labelled with the PET-tracer, 22 Na. Subsequently, it was closed with a screw, establishing a physically sealed source. The sample was then stored at 20°C. Beginning daily, with increasing time lag, a sequence of 20 PET images was produced over a period of 150 days until the tracer was roughly equally distributed throughout the core. In the following, an initial period of 1 month was considered, before the spreading tracer reached the sample circumference.
The raw data ( projections) were corrected for attenuation and scattering of the photon pair. The tomographic image was reconstructed using the iterative OSMAPOSL algorithm supplied by Opensource STIR-software (Thielemans et al., 2012). The remaining concentric-ring artefacts, caused by large gaps between the detector modules, were reduced using image-processing methods. The images were decay-corrected, and the total activity was calibrated with respect to the administered activity, thus yielding quantitatively activity concentrations per voxel.
The accuracy of the method was evaluated with the help of Monte Carlo PET simulations. The initial tomogram (image taken at 3 h after the start of the experiment) was reproduced successfully by the Monte-Carlo simulation (Zakhnini et al., 2013). However, imaging artefacts were still present in the tomograms, caused by multiple factors including tracer distribution, spatial distribution of the signal-to-noise ratio, and the precision of applied corrections. Improvements of imaging quality are more elaborate in dense geomaterials than in in light biomaterials of medical PET applications. Further success at eliminating artefacts is being achieved on an ongoing basis.
The sample had to be considered as damaged by stress and fluid-pressure release during coring, fluidloss during transport and storage before the experiments, and by preparation. In principle, it is possible to minimize these effects by conducting the experiments in a pressure vessel and by reducing the preparation impact. However, the aim of the present study was to optimize the measurement procedure and to develop evaluation methods.

R E S U LT S A N D D I S C U S S I O N
The diffusive propagation of the tracer is shown in four time steps as orthogonal slices of spatial activity concentration, from the initial image to that taken after 27 days, when the tracer reached the circumference (Fig. 2). At the bottom of each image, the respective projections in the axial direction of the mean value are shown. Isosurfaces of activity concentrations of the tomogram taken after a period of 10 days are shown in Fig. 3. At this time step, a transversely isotropic finiteelement model (FEM) was fitted to the data collected, applying the Comsol Optimization module (Fig. 4). This method (Schikora, 2012;Gerasch, 2015) yielded the apparent diffusion coefficients in the principal directions, D a,xx = D a,zz = 1.1E-10 m 2 /s and D a,yy = 3.4E-11 m 2 /s, which is in accordance with results Distant from the source, in the range of low activity, the axial projections are smooth and elliptically shaped, in accordance with the FEM result. These low-activity values are more affected by background noise than the higher-activity values adjacent to the source. In the zones with higher-activity values, local structures become visible which deviate from the smooth elliptical shape. In particular, the perimeter of the isosurfaces appears oblate, in the shape of rounded planes, which also are slightly eccentric from the bore (Fig. 5). These deviations from simple anisotropic propagation are indications of heterogeneity, probably because of a fine silty layer adjacent to the bore. The diffusional transport appears to be controlled by laminations, lenses and other structural heterogeneities that are present in the OPA clay, causing preferential diffusion pathways at the cm scale.
C O N C L U S I O N S Compared to input-output methods for deriving diffusion parameters, process tomography of diffusion in clays has certain advantages: continuous sampling of fluid samples containing radioisotopes is not required, the complete spatio-temporal evolution of the tracer distribution pattern is recorded and samples with representative sizes with respect to major structural features can be investigated non-destructively at simulated temperature and pressure conditions. Meaningful anisotropic diffusion coefficients were derived. In addition, indications of heterogeneity effects were found, i.e. alignment of zones with greater concentrations of tracer in contrast to blank zones with small concentrations of tracer. The nature of these zones will be clarified after decay of the tracer, using µCT and destructive methods.
This spatial variance of diffusion at the cm scale will have implications for our understanding of the fundamental process and for diffusional transport modelling. Zones of the material which are dominated by high tortuosity and a larger surface-to-volume ratio are less effective for diffusional transport than zones with lower tortuosity and a smaller specific surface area. Thus, total transport is amplified and retention is reduced with respect to homogeneous conditions. FIG. 5. Accentuated display of heterogeneity of the tracer distribution after 20 days: iso-surface and axial projection (right) and cylinder slice at r = 12 mm (left).