Quantitative experimental monitoring of molecular diffusion in clay with positron emission tomography

Clay plays a prominent role as barrier material in the geosphere. The small particle sizes cause extremely small pore sizes and induce low permeability and high sorption capacity. Transport of dissolved species by molecular diffusion, driven only by a concentration gradient, is less sensitive to the pore size. Heterogeneous structures on the centimetre scale could cause heterogeneous effects, like preferential transport zones, which are difficult to assess. Laboratory measurements with diffusion cells yield limited information on heterogeneity, and pore space imaging methods have to consider scale effects. We established positron emission tomography (PET), applying a high-resolution PET scanner as a spatially resolved quantitative method for direct laboratory observation of the molecular diffusion process of a PET tracer on the prominent scale of 1–100 mm. Although PET is rather insensitive to bulk effects, quantification required significant improvements of the image reconstruction procedure with respect to Compton scatter and attenuation. The experiments were conducted with 22Na and 124I over periods of 100 and 25 days, respectively. From the images we derived trustable anisotropic diffusion coefficients and, in addition, we identified indications of preferential transport zones. We thus demonstrated the unique potential of the PET imaging modality for geoscientific process monitoring under conditions where other methods fail, taking advantage of the extremely high detection sensitivity that is typical of radiotracer applications.


Introduction
Natural clay typically has a heterogeneous composition and a spatially variant anisotropic structure.Due to grain sizes of the clay fraction in the micrometre range, pore sizes are extremely small.This effectively inhibits advective flow and causes high internal surface area, forming natural geological barriers.However, natural sediments have to be considered as heterogeneous material.In particular, in Opalinus Clay (OPA) we not only observe heterogeneity on the micrometre scale (Keller et al., 2011), but also laminations and concretions on the millimetre scale (Fig. 1) which potentially provide either preferential diffusion pathways or specific sorption sites, respectively.Therefore, standard millimetre-sized samples in diffusion cell experiments cannot be considered as homogeneous.As smallest possible size of the representative elementary volume (REV), we consider the size of standard drill cores because their size is just above the largest observed heterogeneities.It still is an unanswered question as to whether a REV exists beyond which the material may be considered as homogeneous, or whether analogical heterogeneities exist on larger scales.
This question is relevant for a number of safety cases because clay layers are generally considered as a geological barrier for the safe enclosure of hazardous substances.Preferential zones with higher diffusion rates generally are less tortuous, which also means a smaller effective internal surface area.Therefore, these zones delimit the barrier function both by increased diffusion rate and decreased retention, which have to be considered as a worst case scenario.
Determination of these heterogeneous parameters with common diffusion-cell experiments is one possible approach, but it is laborious and requires a large number of studies in order to match these delimiting zones (Van Loon et al., 2004).Another approach is based on the detailed characterization of the pore space using multiscale imaging, applying X-ray µ-CT (e.g.Kaufhold et al., this issue), and imaging modalities with higher resolution (Hemes et al., 2015).These methods directly reveal structural heterogeneities.However, derivation of diffusion parameters from structural imaging relies on involved modelling.It has to consider hierarchical structures over many orders, from surface characteristics on the molecular scale to fine layering on the centimetre scale.Benefits of this method are readiness because it does not rely on long-lasting experiments, and comprehensiveness, in particular by ab initio process simulation.However, the complexity of the models requires verification, preferably by direct spatio-temporal observation of the process.
We introduced a method for quantitative visualization of the diffusing species (Kulenkampff et al., 2015), which can serve as an experimental complement and for the verification of structural models.We applied positron emission tomography (PET) as a totally selective and the most sensitive imaging method for positron-emitting radionuclides.PET responds to the spatial distribution of positron decay events, which are detected as coinciding antiparallel photons with energy of 511 keV.The mass attenuation is at a minimum in this energy range, which allows the reconstruction of the spatial distribution of positron decay events in rock samples up to diameters of 100 mm.
Since the introduction of PET as clinical imaging modality, a number of PET studies have been published on fluid flow in rocks, including dispersion phenomena (e.g.Kulenkampff et al, this issue).These studies have been conducted generally with short-living PET tracers and clinical PET scanners with rather low spatial resolution.Our application of non-standard long-living PET tracers and a high-resolution scanner, which was originally designed for biomedical research, is a means for long-term studies at higher spatial resolution.This is a prerequisite for monitor-ing slow molecular diffusion processes without other driving forces than the concentration gradient of the tracer.
However, strong blurring effects of Compton scattering and attenuation of the 511-keV-radiation have to be considered, which could delimitate the significance of the method.Earlier, we identified these effects as an issue that must be considered with regard to quantification of the tracer concentration (Zakhnini et al., 2013); the resulting images rather had to be considered as qualitative.Scattering and some deficiencies of the reconstruction procedure had caused blurring and imaging artefacts, and the fit of a conceptual diffusion model to these data had produced unexplainable high-diffusion coefficients (Kulenkampff et al., 2012).These issues could be solved to a large extent with a proper scatter correction procedure that is outlined in this article, and we are now able to make quantitative use of it.
2 Materials and methods

Measurements
An OPA horizontal drill core (diameter 100 mm, length 80 mm) from Mont Terri was cast in epoxy resin (Fig. 1).The sample is not considered to be undisturbed because fracturing due to the stress release and drying during storage is probable.This situation is representative rather of the excavation damage zone than of undisturbed host rock in the far field.
Perpendicular to the bedding, an axillar blind hole (diameter 5 mm, length 50 mm) was drilled into the core and filled with synthetic OPA porewater (Pearson et al., 2003).After an equilibration period, which included observations using 124 I-labelled OPA water, the hole was filled with synthetic OPA water, now labelled with 22 Na, with an initial activity of 18 MBq.Then it was closed with a screw, establishing a physically sealed source.The sample was stored at 20 • C. Beginning daily, with increasing time lag, a sequence of 20 PET images was recorded over a period of 143 days until the tracer was roughly equally distributed throughout the core.The length of these frames was 40 min.
We applied a ClearPET scanner manufactured by Raytest (Sempere Roldan et al., 2007).This is a highresolution PET scanner, designated for biomedical research on small animals, that was originally developed by the Crystal Clear Collaboration (http://crystalclear.web.cern.ch/crystalclear/Default.html).In contrast to clinical PET scanners with a spatial resolution of 3-5 mm, this class of scanners reaches the fundamental physical resolution limit of about 1 mm, depending on material density; the voxel dimensions are 1.15 mm, accordingly.This high resolution is achieved at the expense of technically unavoidable gaps between the detector crystals, which cause an inhomogeneous sensitivity matrix.This loss of homogeneity is mitigated by rotating the gantry.

Image reconstruction
The raw PET data are projections of the spatial tracer distribution.These projections correspond to lines of response (LORs) connecting the detection points of the antiparallel photon pair that was transmitted by the annihilating positronelectron pair.These photons undergo attenuation and Compton scattering, which are both controlled by the mass attenuation coefficient, which depends directly on material density.Because the density of geomaterials is considerably higher than of biological tissue, we have to consider corrections for these effects more thoroughly than in common biomedical or clinical PET applications.Details of the method are reported in Kulenkampff et al. (2016).Here, we focus on the implementation and calibration of the scatter correction procedure, which showed significant impact on the experimental results.
In Zakhnini et al. (2013), we applied a Monte Carlo (MC) simulation procedure for determining the impact of scattered events and for developing a tentative correction procedure.This procedure is extremely demanding with respect to computing resources and takes at least several hours for each image.The MC simulations were conducted with the open-source simulation platform GATE (Jan et al., 2004).Input parameters are scanner geometry, sample geometry ("phantom"), material parameters, and source parameters.The source parameters include geometry, activity, and the significant contributions of the decay spectrum of the particular nuclides.In addition to the experimentally applied nuclides 18 F, 22 Na, and 124 I, we also considered 48 V, 58 Co, 64 Cu, 86 Y, and 132 La as possible tracers for studying geochemical effects.These isotopes are either commercially available or can be produced with our in-house cyclotron, Cyclone ® 18/9 (IBA Molecular, Belgium) (Mansel, 2015;Mansel et al., 2014).
In contrast to the real physical world, MC provides a means to study the history of all recorded events, from the positron decay to the detection.This includes a comprehensive parametrization of scattering, taking into account the magnitude and origins of scatter, and the distributions of energy, the order of scatter, and scattering angles.We also studied the effects of different potential PET nuclides because different initial positron energies and additional gamma radiation have different impacts on the recorded data.We have to distinguish -N tot : total number of coincidences; -N rand : number of random coincidences (coincidences that are not from the identical positron decay event); -N sc : number of scattered coincidences; -N trus : number of true (non-random) unscattered coincidences; -N fals : false coincidences, coinciding spurious gamma radiation.
The global MC results with respect to the nuclides are depicted in Fig. 2  by material effects (N trus /N tot ) is considerably higher in water (52.3 %) than in clay (31.8 %).This is also reflected by the scatter fraction The high scatter fraction, in combination with attenuation up to 80 %, is the reason for the annoying blurring effect.The deviations of the LORs from the straight connection between source and detector positions are on the order of centimetres.This contradicts the presupposition of the reconstruction algorithm and causes considerable circular imaging artefacts.
The mean order of scattering (number of Compton interactions per registered coincidence) is more than 3; therefore the premises for single-scatter modelling as a basis for a scatter correction procedure are not satisfied, and the applicability has to be verified.However, the mean scatter angle is in the range of small-angle scattering (below 10 • ), with maximum values around 20 • .
In contrast to the preceding paper (Zakhnini et al., 2013), we abandoned the version 1.4 of STIR, which was the originally supplied image reconstruction library, and switched to the recent version 3.0 (Thielemans et al., 2012).As before, we applied corrections for random coincidences and attenuation.The normalization process was improved according to Weber et al. (2006) in order to better consider the effect of void bins in the projections caused by the gaps between the detectors.
The most recent version of the reconstruction software, since STIR 2.1 (Tsoumpas et al., 2004), supplies a simpler and faster analytical scatter modelling method than MC that is based on Watson et al. (1996).It is an approximation of the deviations of the coincidences along each LOR according to the Klein-Nishina equation for Compton scatter.Scatter correction of multiple time frames becomes practicable with an optimized procedure that typically requires 30 min on a standard CPU for each frame.
This single scatter simulation algorithm (SSS) does not account for multiple scattered events and scatter from sources and scatter points outside the field of view (FOV).It requires the distribution of mass attenuation coefficients, as a measure of the scatter cross section, and an estimate of the source, usually an uncorrected image, as input data.Multiple scattering is approximately considered by scaling of the scatter estimate and by iterative application of the procedure (Polycarpou et al., 2011).The customary calibration of the SSS correction is then conducted by fitting the simulated scatter events to the measured events outside the region of the source distribution.Because of the large number of void bins in this marginal region of the FOV of the ClearPET scanner, this method is unreliable.We therefore take advantage of the results of the MC simulations, with which we are able to determine the scatter fraction SF quantitatively, and vary the scaling factor of the SSS iteratively until the scatter fraction from the MC simulation is matched (Fig. 3).In order to avoid overcompensation of the scatter, we assume a target value which is 10 % lower than the simulated scatter fraction because the simulation includes additional sources of scatter (e.g. in the detector crystals and in material outside the FOV).
In order to calibrate the scatter correction model, we thus either need an estimate of the scatter fraction, which is the result of one single MC simulation of a rough model, or we consider a standard calibration factor of 5 for our scanner as a conservative estimation.
The outcome of the SSS is in general accordance with the result of the MC simulation.We reconstructed images from the MC data sets:  22 Na diffusion in an OPA drill core (diameter: 100 mm, length 80 mm).Each frame is depicted as isosurface of the one-tenth maximum value (grey), three horizontal slices through the source region, and the axial (vertical) maximum projection at the bottom.The colour ranges are scaled frame-wise; the maximum value of the colour range is half maximum total amplitude (for the corresponding amplitude distribution, see Fig. 6).
all true coincidences tot (no scatter correction) unscattered unsc (reference image) only scattered coincidences sc -SSS-corrected with scaling factor 5 and 10.
In Fig. 4 we computed the circular mean values of the axial projections of these images, yielding the intensity vs. the radial distance from the central line source in clay.The impact of scattering becomes clear at distances above 5 mm, and without scatter correction, the intensity becomes dominated by scatter (tot and sc).The SSS correction scaled with a factor 10, which is the intersection point in Fig. 3, yields data below the line of the reference image, which therefore  are overcorrected.With a factor 5 they are undercorrected, and scatter effects dominate in regions far from the source.We assume that this is the region of high scattering angles and multiple scatter that is not considered by the SSS algo- rithm.It should be noted that the ordinate of Fig. 4 is scaled logarithmically, and that these deviations are more than 2 orders below the source intensity.We recommend slight undercorrection to avoid the creation of additional void bins.

Results
As example, we report the measurements on one sample (BLT 137/3): Beginning daily, with increasing time lag, we produced a sequence of 20 PET images over a period of 150 days until the tracer was roughly equally distributed over the core (Fig. 5) (t = 0, 3, 6, 10, 13, 16, 20, 22, 27, 31, 35, 41, 48, 55, 69, 93, 112, 127, 143, 161 days).Because of the large amplitude range, the images were scaled frame-wise.Figure 6 shows the evolution of the maximum value, which is the scaling reference, and the 90 % and 99.9 % quantile.The first frame of Fig. 5 shows the clearly distinguishable initial source distribution, i.e. the 5 mm blind hole, filled with the labelled synthetic OPA porewater.In the subsequent frames, the tracer concentration in this source zone decreases and the vertical slices show the roughly elliptical spreading of the tracer distribution into the clay.The major axes of these ellipses are identified as the anisotropy axis of the diffusion process.A gas bubble is visible in the first frame, which at the second frame had moved from its initial position at the bottom of the hole.Figure 7 shows vertical slices through the three-dimensional data set in the direction of these major anisotropy axes of seven selected frames (t = 0, 3, 6, 13, 20, 35, 143 days) with a uniform colour scale, in order to display the decrease of concentration with time.The maximum projection in vertical direction is plotted as height map at the Figure 9. 15 PET frames of 123 I diffusion in an OPA drill core (diameter: 100 mm, length 80 mm).Each frame is depicted as isosurface of the one-tenth maximum value (grey), three horizontal slices through the source region, and the axial (vertical) maximum projection at the bottom.The colour ranges are scaled frame-wise; the maximum value of the colour range is half maximum total amplitude (for the corresponding amplitude distribution, see Fig. 10).
bottom of each frame.It indicates that the diffusing tracer reaches the sample surface after 35 days at frame 10.Then, the tracer distribution homogenizes over the sample; therefore, the quantiles stabilize (Fig. 6) and the anisotropy of the distribution decreases until frame 18 after 143 days, which shows almost homogeneous tracer distribution.
Figure 8 is an enlargement of the maximum projection of frame 9, after 1 month and just before the spreading reaches the circumference of the sample.The shape of this distribution appears as a roughly 10 mm thick rectangular block, rather than as a smooth anisotropy ellipsoid, notwithstanding remaining indications of circular artefacts.This type of shape could be interpreted as an indication of diffusion along fine layering, rather than homogenous transversal anisotropic behaviour.However, these findings should be confirmed by additional investigations (e.g.µ-CT, radiography and analysis of thin sections, as well as model simulations) and after further improvement of image quality.
Prior to the investigations with 22 Na, we tested 124 I as a PET tracer (Fig. 9).In principle, its decay time (4.176 days) allows observation periods up to 1 month, sufficient for parameter derivations according to the procedure outlined below.However, the images are noisy because of other "parasitic" γ radiation (false coincidences), and the decreasing count rate causes an increasing error with time.Therefore, application of 124 I merely is a qualitative test method, rather than being suited for quantitative estimations.

Discussion
In the past we have experienced problems with the correction methods that had been supplied with the scanner because these had been designed for material with low attenuation and scatter.Attenuation model and scatter correction had been derived from measured data, a procedure that increased the noise level unacceptably.Instead, we implemented a method which is based upon attenuation and scatter modelling.Their establishment is simple because of the plain geometry of the samples.The expedience of the STIR-SSS for scatter modelling has been proven with MC simulations, and calibration factors have been found.
With these corrections we could reconstruct clear and sharpened images which retained a minimum scatter component.Features became visible that had been covered by blurring effects, and new levels of detail were identified.A gas bubble became clearly distinguishable.The extension of the tracer cloud, which had been apparently widespread due to scatter effects, became localized and bounded.An optimum fit of an updated finite-element model now yields anisotropic diffusion coefficients in the range of literature results (D xx = D zz = 10.2e−11 m 2 s −1 and D yy = 2.2e −11 m 2 s −1 ) (Lippmann-Pipke et al., 2016).We therefore conclude that the only transport process is molecular diffusion.This is in contrast to our first hypothesis, which was misled by inadequate consideration of the scatter effect, where we presumed that suction into partially saturated zones of the sample caused accelerated spreading of the tracer.
We fitted anisotropy ellipses to the axial projections of the data in order to determine direct experimental anisotropy parameters and to estimate the propagation velocity of the tracer front.The position of this front was defined as the lo- cation of the FWHM line (full width half maximum).From Fig. 11 we read a propagation velocity of 0.5 and 2 mm d −1 for the smallest and the largest axis, respectively.This is in accordance with results that were calculated from literature values by Lippmann-Pipke et al. (2016).
We also computed the ratio b/a of these major axes of the FWHM for all frames, as shown in Fig. 12. Apparently, this anisotropy ratio is not constant, but increases steadily with time until the tracer propagation front reaches the sample surface after 31 days.This initial increase of anisotropy could be caused by preferential diffusion continuing along a rather linear pathway, in contrast to homogeneous anisotropic behaviour.The subsequent decrease of b/a is caused by increasing equilibration of the tracer concentration over the complete sample.
Figure 8 may explain this situation.Here, the amplitude distribution appears as a roughly rectangular region with a thickness of about 2 cm, rather than as a lineament.Therefore the anisotropic behaviour could be caused by a thin layer with elevated diffusion properties.
It is not yet possible to compare these findings with structural µ-CT images before the tracer activity had decayed.To date, we have therefore not been able to correlate the observed diffusion pattern with alternating sandy and clay layers that are present in Opalinus Clay.This will be the topic of a subsequent study.
As heterogeneous effects are controlled by the topology of the pathways, they are only observable on sufficiently large samples.Nevertheless, these effects should be taken into account for transport simulations on the field scale because preferential diffusion zones could accelerate diffusional transport, not only by increased directional transport.Zones which are not effectively taking part in the diffusional transport process are inefficient with respect to sorption as well.Therefore, the retention parameters decrease simultaneously with the increasing diffusional transport through preferred zones.This could further amplify the acceleration of the diffusional propagation.We therefore recommend conducting diffusion experiments on the decimetre scale, where these heterogeneous effects emerge, preferably with PET.

Conclusions
In the past it was shown that PET is applicable as a nondestructive method for quantitative determination of the spatio-temporal tracer distribution in order to observe flow phenomena in complete drill cores.We have now established PET as a tomographic method for the imaging of molecular diffusion.We take advantage of the extreme sensitive detection of decay radiation, which is a beneficial characteristic of radiometric methods in general.In particular, we applied long-living PET tracers, a high-resolution PET scanner, and specially adapted correction methods of the raw data, in order to gain suitable quantitative tomograms of the temporal evolution of the tracer distribution.From these data, integral anisotropic diffusion coefficients can be computed that are comparable to the results of a large set of single onedimensional diffusion-cell experiments on small specimens.Beyond that, it is presently the only quantitative method which enables the effects of heterogeneities to be studied on the centimetre scale, in particular preferential diffusion zones.Such zones could have a significant impact on the diffusional propagation of chemical species on the field scale.

Data availability
The voluminous spatio-temporal image data were archived at the HZDR data repository which is currently prepared for open access.The HZDR has registered at datacite.org to provide open access data.Using the short name of the HZDR (TIB.HZDR) -data will be accessible via the datacite search (http://search.datacite.org).The official HZDR DOI-Prefix is doi:10.14278.

Figure 1 .
Figure 1.Left: OPA drill core before it is embedded in resin ("plug size" refers to the dimensions of the µ-CT image).Right: µ-CT image of a mini-plug from the same formation (recorded at the Federal Institute for Material Research and Testing BAM, Berlin). Figure taken from Kulenkampff et al. (2015).

Figure 3 .
Figure 3. Relative number of scatter corrected coincidences N corr (c)/N tot vs. the scatter calibration factor c. The horizontal reference line was derived from the MC simulation (Fig. 2, Na-22).

Figure 4 .
Figure 4. Intensity I vs. radial distance r from the central line source of images, in dependence of the type of scatter correction (MC_tot: total true coincidences tot, MC_unsc: only unscattered coincidences unsc, MC_sc: only scattered coincidences sc, SSS_Fakt_10/SSS_Fakt_5: all true coincidences (tot) with SSS correction and scaling factors 10 and 5).

Figure 5 .
Figure5.20 PET frames of22 Na diffusion in an OPA drill core (diameter: 100 mm, length 80 mm).Each frame is depicted as isosurface of the one-tenth maximum value (grey), three horizontal slices through the source region, and the axial (vertical) maximum projection at the bottom.The colour ranges are scaled frame-wise; the maximum value of the colour range is half maximum total amplitude (for the corresponding amplitude distribution, see Fig.6).

Figure 6 .
Figure 6.Maximum, 99.9, and 90 % quantiles of the amplitude vs. frame time of Fig. 5.The colour scale in Fig. 5 refers to the maximum value.

Figure 7 .
Figure 7. Vertical slices and axial projection of selected frames of Fig. 5, with uniform colour scale.At t = 35 days, the tracer arrived at the cylinder surface.

Figure 8 .
Figure 8. Axial maximum projection of frame 9 (27 days), when the tracer distribution reached the circumference.

Figure 10 .
Figure 10.Maximum, 99.9, and 90 % quantiles of the amplitude of Fig. 9 vs. frame time of Fig. 9.The colour scale in Fig. 9 refers to the maximum value.

Figure 11 .
Figure 11.Lengths of the major axes of the FWHM vs. frame time (red: fast axis, green: slow axis) of the 22 Na experiment.

Figure 12 .
Figure12.Left: axial projection of the activity distribution with elliptical fit of the FWHM curve of the ninth frame after 27 days (see Fig.8).Right: temporal evolution of the ratio of the major anisotropy axis.