Geoscientific process monitoring with positron emission tomography ( GeoPET )

Transport processes in geomaterials can be observed with input-output experiments, which yield no direct information on the impact of heterogeneities, or they can be assessed by model simulations based on structural imaging with μCT. Positron emission tomography (PET) provides an alternative experimental observation method which directly and quantitatively yields the spatiotemporal distribution of tracer concentration. Process observation with PET benefits from its extremely high sensitivity together with a resolution that is acceptable in relation to standard drill core sizes. We strongly 10 recommend applying high-resolution PET scanners in order to achieve a resolution in the order of 1 mm. We discuss the particularities of PET applications in geoscientific experiments (GeoPET), which essentially are due to high material density. Although PET is rather insensitive to matrix effects, mass attenuation and Compton scattering have to be corrected thoroughly in order to derive quantitative values. Examples of process monitoring of advection and diffusion processes with GeoPET are illustrating the procedure and the 15 experimental conditions, as well as the benefits and limits of the method.


Introduction 1.Aims
It is common practice to observe and quantify transport in geological materials with input-output experiments, which treat the material as a black box.However, these methods miss heterogeneous effects, like preferential transport, and are also disadvantageous when the process is slow, or when an applied tracer is retained by internal reactions.One approach for unravelling transport processes in geological materials is a very detailed description of the structure and composition of the matrix and the pore space.These provide boundary conditions for model simulations whose final results are compared and validated with input-output experiments.Such model simulations require specific geochemical parameters for the tracer species and the particular material, which are derived from other experiments.
Direct observation of the propagating species during the course of the process is a complementary approach which yields more intrinsic information on the process.Such process-tomographic methods are particularly suitable when complexity or fineness impede the compilation of an applicable structural model.
Positron emission tomography (PET) was originally a nuclear medical method that takes advantage of the simple but most sensitive detectability of radiotracers.We have striven, successfully, to broaden the scope of PET applications to geoscientific and technical use for more than 15 years.The potential of the method was proven and communicated as an important experimental method for enhanced understanding of geochemical processes in soils and rock formations, for modelling and upscaling.

Process tomography modalities
Geological materials are generally heterogeneous with respect to structure and composition.Moreover, the internal structure is highly complex, with contributions from hierarchical scaling (fractal behaviour) and non-scaling contributions on a particular scale.On larger scales, this behaviour is considered by anticipating a representative elementary volume that is supposed to combine all heterogeneities in a ho-mogeneous parameter set.Disregarding tomographic methods, the characteristic size of this volume had to be derived from laboratory experiments on a large number of differently sized samples.Structural tomographic methods, like µ-CT, now allow for the characterization of the structure comprehensively and over a large scale.Then, model simulations can be conducted on these structures, taking the parameters of the particular process into consideration.However, the huge size of hierarchical tomographic data currently precludes composition and processing of pore-scale simulation models in detailed structures.Currently, typical model simulations are conducted on structural models that are based on the limited resolution range of a particular tomographic method.
This type of bottom-up approach could miss significant features of transport processes because the chemical species are progressing in a highly heterogeneous material with disorder from the molecular to the centimetre scale.Often, preferential transport and reactions at distinct sites of the internal surface have to be considered, with major contributions from structures below the resolution of the tomographic modality.These limitations with respect to transport processes could be overcome by direct tomographic observation of the transport processes.
The number of applications of industrial process tomography is rapidly growing.Industrial process tomography means observation of a large variety of processes by spatio-temporal monitoring of process parameters for process understanding and control (Wang, 2015;Williams and Beck, 1995).Applications of these methods are also growing in geosciences; some of these methods were even originally developed for geoscientific purposes, like electrical resistivity tomography.
The tomographic observation of substances in disjoint phases is possible with a large number of imaging modalities that respond to specific properties of the phases, like CT (density) and ultrasonic imaging (elastic velocity), but they rely on significant contrast.CT data are commonly processed in order to produce structural images, segmenting the grey-level images into two or multiple domains with different density.Therefore, the most significant characteristics in facilitating the segmentation are spatial resolution and homogeneity of the mapping.
A major class of processes in the geosphere are transport phenomena, including advection, dispersion, and molecular diffusion of dissolved substances or particles.Frequently, the most meaningful observable is the concentration of dissolved or dispersed species, which is an intensive thermodynamic parameter that is most relevant for aligning observations with simulation results from models of the processes.Concentration is a continuous value, referring to a finite test volume which is defined by the spatial resolution of the measuring method.Here, spatial resolution is merely the response function which controls the smoothing of the image.
Ionic concentrations can be monitored with different types of electrical impedance tomography (EIT).These are frequently applied to investigate transport processes of conduc-tive solutions in soils and rocks.However, electrical resistivity, as one of these parameters, depends on a variety of material parameters and rock-fluid interactions; therefore the result is not univocal.Other imaging methods apply tracers for labelling, which facilitates detection and quantification of the species (Hevesy and Paneth, 1913;Garrett, 1963).Spatio-temporal visualization of labelled substances requires tomographic methods with selective response to the tracer.In the field of biomedical research, these techniques are called molecular imaging methods (Weissleder et al., 2010).Out of these methods, magnetic resonance imaging (MRI), single photon emission computed tomography (SPECT), and positron emission tomography (PET) are suited for opaque media.In contrast to EIT, these methods selectively respond to the concentration of protons or specific tracer isotopes, but they are mere laboratory applications and are barely suitable for field applications.The methods differ with regard to the matrix effect: the applicability of MRI is generally restricted by distortions of the magnetic field which are caused by paramagnetic compounds of the matrix.Compared to PET, SPECT offers a wider choice of applicable radiotracers and is not restricted by a fundamental limit of resolution.However, the due consideration of radiation attenuation and scattering in the material is currently constricting the applicability, spatial resolution, and quantification of SPECT.PET is least affected by material interferences.Boutchko et al. (2012) presented applications of both tomographic nuclear imaging modalities.
PET is perfectly selective and by far the most sensitive method because the response is directly and only related to the number of decaying positrons, and thus also the number of tracer atoms.Depending on the source distribution, it only requires a small number of events (on the order of 100 per voxel) to reconstruct the tomogram.A typical voxel size around 1 mm means a cube of water containing about 10 19 atoms.This means that the sensitivity is better than picomoles per microlitre, depending on the decay characteristics of the PET nuclide and the noise level.With this ultimate sensitivity, PET can serve as a "gold standard" for tomographic imaging of tracer concentrations.Its spatial resolution of 1 µL is appropriate for typical sample sizes of 1 L.
Currently, structural imaging using µ-CT brings a great increase in the understanding of transport processes in geomaterials.The geometry of the pore space from segmented µ-CT images provides realistic boundary conditions for model simulations.Commonly, the results of such simulations are verified with rather simple input-output experiments on macroscopic samples that only consider the integral response to an input signal.Instead, we promote the application of process observation with PET for directly quantifying mass transport inside the sample.
PET is a particularly suitable method when mass transport through argillaceous rocks is dominated by the small sizes of the pores, which complicate the creation of a realistic model.Methods with very high spatial resolution, like FIB-SEM or nano-CT, are able to visualize the pore structure in the nanometre range, but they are only applicable to millimetre-sized samples.Hence, these are not suited to consider structures in the millimetre range, like fine layering and other heterogeneities.Furthermore, N 2 adsorption porosimetry (e.g.Nagra, 2002) suggests that a significant portion of the pores falls below the resolution of even the highest resolving imaging methods.This means that connecting pores are not segmented and therefore excluded from the structural model.Here, PET is a dependable method for monitoring the propagation of the tracer through the fine-grained material, although its spatial resolution is far beyond the largest pore structures.
Most of these geoscientific PET studies have been conducted with clinical scanners of cooperating medical departments.Typically, this limits the possibilities for hardware and software tuning, constrains the spatial resolution to 3 to 5 mm, impedes long-term studies (longer than some hours), and application of proprietary software might distract from pushing forward the development of algorithms for enhanced image quality.
Our work group systematically developed and applied PET scanners for geoscientific research in the course of over 1 decade (Richter and Gründig, 2000;Richter et al., 2005;Wolf et al., 2010;Gründig et al., 2002Gründig et al., , 2007;;Kulenkampff et al., 2008Kulenkampff et al., , 2010Kulenkampff et al., , 2013Kulenkampff et al., , 2015;;Barth et al., 2014a, b).The first published experiments were conducted with a selfconstructed PET scanner for soil columns and rock cores, designed for slow transport studies (Richter and Gründig, 2000).We have also applied clinical PET and PET/CT scanners, and since 2007 we have used a high-resolution PET scanner (ClearPET by Raytest, Germany), designed for biomedical research on animals (Sempere Roldan et al., 2007).

Principle
PET responds to the emission of radiation caused by positron decay.A positron-emitting isotope is applied as a radiotracer for labelling the substance which propagates through the sample during the observed process.
In order to compute the concentration c of a chemical species in one voxel, we have to consider the decay law Eq.( 1) with time t, decay constant of the tracer λ, and the number of decaying atoms N (Lieser, 2001): (1) With the activity a, we can calculate the concentration of the decaying nuclide (with the voxel volume V and Avogadro's number N A ) from the measured activity: Often, the tracer isotope is used for labelling a carrier solution with the natural isotopic composition, in order to avoid specific interaction effects acting exclusively on the tracer.Then the concentration of the labelled carrier solution is determined by the fraction of tracer atoms in the carrier solution.
Positron-emitting isotopes of virtually all natural elements exist.However, only a small fraction of theses nuclides can be produced with acceptable expenditure, and only a few nuclides are practically useful because decay time, positron energy, and other types of decay radiation have to be considered (Hawkesworth and Parker, 1995;Richter et al., 2005).We applied a number of GeoPET isotopes, which are listed in Table 1.
PET makes use of the coincidence detection of the photon pair produced by the annihilation of the positron which occurs when a slowed down positron interacts with an electron.The initial kinetic energy of the positron determines its free path length in matter.Therefore, the positron energy determines the fundamental limit of resolution, which is on the order of 1 mm in water and decreases with increasing electron density (Levin and Hoffman, 1999).
The annihilation radiation of the positron has the remaining energy of an electron, 511 keV.Because of the poor energy resolution of PET scanners, background radiation or portions of the γ spectrum of the nuclide may interfere with the coincidence detection.Particularly, γ energies in excess of 1022 keV may produce positron-electron pairs anywhere in the material, with subsequent annihilation of the positron without relation to the tracer position ("false" events).These false events are correlated and remain unfiltered by the coincidence detection.This is in contrast to background radiation, single γ quanta, which produces a background noise of random coincidences.
The radiation is usually detected in detector crystals that are mounted on multichannel photomultipliers.The detection probability depends on the crystal parameters, is rather small, and is also strongly reduced by inevitable gaps between the crystals.Two events which are detected in opposite detector positions, within an energy window around 511 keV, e.g.350-700 keV, and within a small coincidence window (typically 20 ns), are accepted as coincidences.This type of coincidence filtering is highly efficient, and only about 3 % of the registered events are recorded as coincidences."Delayed" coincidences, coinciding between the reference window and a delayed window with a time shift greater than the coincidence time, provide an estimate of the random coincidences.
The geometrical parameters of the connecting line between the detection positions (line of response: LOR) are stored as projection files of "bins" from which the threedimensional tracer distribution is reconstructed.The total number of required coincidences for image reconstruction is typically on the order of 10 7 , depending on the source distribution.Another algorithm allows fast particle tracking (PEPT) (Parker et al., 1993); when only a few distinct particles are labelled, their position can be reconstructed from a small number, at least two linearly independent LORs.
Currently, the detector crystals of clinical PET scanners are comparatively large, in order to optimize detection efficiency, which also accelerates the examination time and reduces the dose for the patient.As a consequence, these scanners do not accomplish the maximum spatial resolution (typical resolution of clinical PET scanners: 3-5 mm).High resolution is the domain of biomedical PET scanners with a smaller field of view (FOV) that are designed to achieve the maximum resolution (around 1 mm) (Cherry and Chatziioannou, 2004).Generally, the number of coincidences, and thus the frame length, determines image quality.Typically, the maximum frame rate is on the order of minutes and depends on the quantity and distribution of activity in the FOV.
The sample material plays a minor role and eventually may be neglected when the density ρ is low (ρ ≈1 g cm −3 ).In denser material, however, both attenuation and scattering should be considered, although interaction effects between 511 keV photons and matter are low compared to the energy range of common µ-CT sources (Fig. 1).
We have to assume a loss of photon radiation intensity I , which is described by the linear attenuation coefficient µ : In the energy range of annihilation photons the most probable interaction with matter is Compton scattering.The linear Compton attenuation coefficient µ c is related to the Compton Solid Earth, 7, 1217-1231, 2016 www.solid-earth.net/7/1217/2016/Table 2. Attenuation properties at 511 keV of relevant materials (µ/ρ: mass attenuation coefficient, ρ: density, µ: linear attenuation coefficient, q: intensity through material thickness of 5 cm).Data assembled from the XCOM database (Berger et al., 2010).
with the electron density ρ e , which depends on bulk density ρ, mean atomic number Z, mean atomic mass A of the material, and Avogadros's number N A .The mass attenuation coefficients µ/ρ are material parameters in the range of 0.081-0.096cm 2 g −1 (at 511 keV) and can be looked up in the XCOM database at NIST (Berger et al., 2010).Therefore, it is obvious that the main factor controlling attenuation is density ρ.Attenuation coefficients of typical materials derived from this database are tabulated in Table 2. Therefore, in geological materials with a density above 2 g cm −3 we have to consider an attenuation of 80 % through a sample with diameter of 10 cm.Compton scattering means that photons are deflected from the straight LOR, they simultaneously lose a portion of their energy, and they eventually leave the energy window -which is considered as attenuation.The solid deflection angle θ is implicitly described by the Klein-Nishina formula: with the classical electron radius r e (2.818 × 10 −15 m) and α e = E 0 m e ×c 2 ≤ 1 (E 0 : initial photon energy, m e c 2 : energy equivalent of the electron).
There is no simple way to apply this formula for scatter correction.A number of scatter correction algorithms are available (Zaidi, 2001;Basu et al., 2007) which generally were developed for medical PET-applications, where the material density is low.We found that the application of simple algorithms based on some type of filtering of the projections or the resulting images is not appropriate in highly scattering material.More recent algorithms based on approximations of a scatter model are more suited.In order to evaluate such algorithms for dense material, we conducted Monte Carlo (MC) simulations, with which all relevant reactions are considered, and which allow the history of all events to be traced back and scattered and unscattered events to be differentiated and quantified (Zakhnini et al., 2013).MC is a method for computing the PET records corresponding to any source distribution, considering all relevant material properties and the characteristics of the scanner.It thus allows the results of image reconstruction algorithms to be verified and principally can serve as forward-modelling operator for tomographic inversion.

Specifications and requirements
A radioisotope laboratory is required that allows safe handling of activities and high local dose rates.Short lived or unconventional PET tracers have to be produced on-site with a cyclotron, like IBA Cyclone 18/9 in our laboratory (Mansel et al., 2014;Mansel and Franke, 2015); longer living nuclides, like 22 Na, can be purchased.
Our ClearPET scanner was optimized for geoscientific questions.It features a rotating gantry with 20 detector cassettes, each equipped with four 8×8 multichannel detectors and double-layered miniaturized detector crystals (LYSO/LuYAP, both 2 × 2 ×10 mm 3 ), in order to achieve a spatial resolution at the physical limit.The gantry diameter can be adjusted to 13.6 or 23.7 cm, allowing samples to be investigated with the typical dimensions of soil columns and drill cores.Its adjustable FOV has a maximum diameter of 160 mm and a length of 110 mm.Compared to medical PET scanners, we take advantage of higher resolution and sensitivity.In the case of ClearPET, this increase of sensitivity and resolution is achieved at the cost of less homogeneity of the field of view.This adverse effect for image quality is caused by inevitable gaps between the individual detector modules and is somewhat mitigated by a rotating gantry.
Usually, frames with a length of 20-60 min are recorded and later subdivided or merged into frames of appropriate length.Process monitoring is then accomplished with sequential three-dimensional images.The minimum frame rate depends on the number of events which is required for image reconstruction.The minimum frame rate of 60 s, the gantry rotation period of ClearPET, is a typical limit.The maximum observation period is limited by the decay of the radiotracer to about 8 × T 1/2 .
Currently, the orientation of the gantry axis is horizontal; therefore a horizontal orientation of the samples is preferred.In the near future the scanner will be upgraded with a tilting device to allow vertical orientation of the samples.A similar approach was taken by ClearPET Neuro (Weber et al., 2006).
Frequently the samples are formatted as cylinders and cast in epoxy with inlet and outlet ports at the end planes, which are connected with plastic fittings and tubing.The welldefined geometry facilitates the computation of the attenuation and scatter corrections.If pressure vessels are required, low-density material should be used, like carbon fibre reinforced polymers (CFRP), aluminium, or titanium.According to their Compton cross section, steel fittings or casings should be avoided or minimized.However, even with a point source in a steel pressure vessel with 5 mm thick walls, we experienced satisfactory image quality.

Data processing
The original reconstruction and correction software is based upon the open-source STIR library (Thielemans et al., 2012), amended with procedures provided by the manufacturer which are based on the work of Weber et al. (2006).We tuned these procedures in order to consider strong attenuation and scattering.Abundant data storage ("listmode" files of single events) and the use of this open-source image reconstruction system ease adaptation to our particular requirements, but also exacerbate operating comfort.
The OSEM reconstruction algorithm (Jacobson et al., 2000) implemented in STIR is an iterative procedure (iteration step j ) for computing the activity a in one voxel v from the count number n in a set of bins b described by The OSEM algorithm is an accelerated version of the standard EM algorithm (Shepp and Vardi, 1982), which divides the problem into m subsets with index k.When these subsets are balanced, i.e. the sensitivity image is independent of the subsets, Equation ( 7) reduces to which is implemented as the OSEM procedure in STIR.
Applying ClearPET, we have to consider mainly azimuthal gaps between the single detectors, which, due to the rotating gantry, cause strip patterns of the projections (sinograms) (Fig. 2).This limits the validity of Eq. ( 9) to m ≤ 2. The choice of two subsets still accelerates the reconstruction procedure with respect to the EM algorithm (m = 1) and accords widely with the condition Eq. ( 9).
The azimuthal gaps also cause strong variations of the sensitivity p(v) in the denominator of Eq. ( 10), which reduces the robustness of the algorithm.We observed that errors are accumulated in low-sensitivity zones and destabilize the iteration process.
In order to calibrate the image a(v) (activity a per voxel v) in concentration units, scalar factors for frame length, radioactive decay, scanner dead time, and a calibration factor for the total sensitivity have to be applied.The dead time correction considers the efficiency decrease of the detectors with increasing activity and is implemented in the reconstruction algorithm as polygon fit of the count rate vs. activity curve.The decay correction should consider the frame length t fr , when it is long with respect to the half-life T 1/2 : with the start time of the image (t img ) and experiment (t exp ).
Solid Earth, 7, 1217-1231, 2016 The calibration of the total sensitivity could be computed with a predetermined parameter from a calibrated phantom.This neglects the quantitative effects of the corrections and therefore is a rough estimation.A better method is the use of calibration sources together with the sample in the FOV, which are recommendable in any case as fiducial marks for later co-registration with other imaging modalities.

Reconstruction parameters and corrections
The effect of unequal detector efficiencies p n is determined with normalization measurements on homogeneous phantoms; therefore, measuring errors have to be considered.A normalization tomogram is shown in Fig. 3.The noise level was reduced considering symmetries.Additionally, we developed gap-filling strategies in order to reduce the impact of void bins in the projections.Nevertheless, we have to consider zones with low sensitivity, which are prone to higher error levels.
Another component of the sensitivity that has to be determined is attenuation p ρ .In order to avoid additional noise from low-level measurements, it is recommendable to compute it from a geometrical attenuation model.The LORs are compositions of the pathways of the positron pair; therefore they traverse the complete sample, and attenuation is given by the projection of an attenuation image of the sample: with I c,0 and I c the source and recorded intensity, respectively, and the integration variable r running over the connecting line between both detectors at R 1 and R 2 .The attenuation image µ(r) can be constructed geometrically or computed from a low-resolution CT image with estimates of the density and the mass attenuation coefficients, and the attenu-ation probability is In STIR, scatter is considered as an additional process that adds coherent noise and therefore must be subtracted from the measured projections.We apply the single scatter simulation (SSS) algorithm from STIR 3.0 for correction of the scatter effect (Tsoumpas, 2004).This algorithm was verified (Zakhnini et al., 2013;Kulenkampff et al., 2016) for dense material with Monte Carlo (MC) simulations of the PET measurements with OpenGATE (Jan et al., 2004), which quantitatively yield magnitude and distribution of all events.Estimated source distributions, as well as material properties and noise, were considered in the MC simulations in order to compute the global scatter fraction quantitatively.Then, the SSS results were calibrated with this reference value.This solves both the quantification issue of the SSS, which is due to the large number of void bins outside the source distribution, and reduces the enormous computational effort for MC.

Error estimation
In contrast to CT, which is commonly segmented into a binary image, PET reflects the concentration of the propagating tracer as a continuous variable.Here, spatial resolution is the parameter of the response function which controls the smoothing of the image.More meaningful characteristics of PET images than spatial resolution are quantification errors, noise level, and the detection threshold.
A rigorous error analysis of the PET imaging procedure is an intricate task (Prekeges, 2013;Kirov, 2012;Meikle and Badawi, 2005).The analysis starts on the detector level.Although the total number of recorded events is large (e.g. 10 7 ), the matrix of projections is sparse, and the count numbers stochastically respond to the binomial distribution.The reconstruction procedure represents a complicated surjective www.solid-earth.net/7/1217/2016/Solid Earth, 7, 1217-1231, 2016 function which maps data from the projection space of measurements to the image space.This relation can be expressed as a huge system matrix, which requires additional parameters (normalization, attenuation correction, scatter correction) that are also subject to errors.Currently, a strict error propagation analysis of this complex system is not feasible.
A rough error estimate can be calculated from the number of counts N that are projected onto one voxel.This number is large; therefore we can approximately consider the Poisson distribution with a relative error σ = 1/ √ N .It is typically on the order of ±10 % or less and can be minimized by increasing activity or frame duration.
According to Strother et al. (1990), we can compute the signal-to-noise ratio (SNR) for each voxel from the weighted variance of all bins contributing to the particular voxel according to When the numbers of true (tr), random (rnd), and scattered (sc) counts for each bin are known, the variance term can be written as backprojection of all these events: With this variance term we can estimate the spatial distribution of the relative error σ as The expected value tr(v) is the result of the iterative image reconstruction procedure, which applies normalization, and is therefore not as equally scaled as the backprojection, which is not normalized.This is the reason why we refer to the backprojection of tr, not the iterative reconstruction.
However, we also have to consider systematic errors that are caused by deficiencies of the reconstruction algorithm and the parameter values.These are identifiable as imaging artefacts occurring in zones where the sensitivity, i.e. the magnitude of the response function of the mapping from the projections to the image, is low (small amplitudes in Fig. 3).Eventually, these artefacts include zeros (e.g. from scatterovercorrection) or large errors on the order of a factor of 2.
In summary, the error depends on the number of counts per voxel and is typically on the order of 10 %.In zones with lower sensitivity it eventually diverges and produces ring artefacts.

Detection threshold
Apart from the relative error, we have to consider a detection threshold of activity below which the source is not reliably detected.A general threshold is caused by the number of events required to reconstruct one voxel.
With decreasing activity, the number of detected coincidences falls below this threshold where the source cannot be reconstructed.Practically, this threshold depends on many parameters and conditions, in particular, noise level, sensitivity, reconstruction procedure, scatter, and attenuation.We estimated the detection threshold with the help of a Monte Carlo simulation of five point sources with a diameter of 1 mm each, which were distributed along an axial profile in the Opalinus sample, and one with a radial offset.Figure 4 is a reconstructed image of the initial situation, showing some scattering around four point sources that are visible in this cross section.The initial activity of 82 kBq was reduced stepwise by a factor of 2 down to 1.2 Bq.The simulated recording time was 18 min.It should be noted that the maximum count number per bin decreases from initially 81 to 1 at the seventh step (factor of 64).
Figure 5 shows the peak amplitudes of the reconstructed image at the point source locations against the total activity of the sources.One peak disappears at the fourteenth step, representing a total activity of 50 Bq, or 10 Bq of the single point source.At this point, the other sources are still detectable, but with a high quantification error, and are strongly deviating from the proportionality line.Therefore, we define the activity concentration of 10 Bq/voxel as best-case detection threshold, which corresponds to a total of 50 counts per voxel at the given frame length.
The total background coincidence count rate, as deduced from blank measurements, is on the order of 500 cps, or 0.1 cps per detector crystal.Therefore, the detection threshold is just at the order of the background signal.
Above this threshold, the response is linear over at least 4 orders.The significance increases with the number of recorded events and therefore with the frame length, as long as the count rate exceeds the background count rate.The upper bound of the frame rate is a function of the propagation velocity of the tracer and -if possible -should not exceed the transit time through one or two voxels.

Granite fracture (clinical scanner)
As an early example for the application of a clinical PET scanner, we review an experiment from Kulenkampff et al. (2008) which was conducted with a Siemens ECAT HR+ scanner in collaboration with the Clinic and Polyclinic for Nuclear Medicine of Leipzig University Hospital.It was con-ducted on a granitic drill core (diameter 5 cm, length 15 cm) with one large axial fracture, which originates from the Äspö HRL.The core was glued into a plastic cylinder and the fluid was injected through grooves in the end pieces.A µ-CT image by Enzmann and Kersten (2006) yields a mean aperture of the fracture of 0.5 mm, with local enlargements up to 2 mm. 5 mL of a 0.01 M [ 18 F]KF solution was injected stepwise into the sample with a flow rate of 0.1 mL min −1 .Six frames with a length of 15 min were taken with the clinical scanner at stopped flow conditions.
Figure 6 shows the segmented fracture from the µ-CT and the propagation of the tracer cloud that was recorded with the PET scanner and reconstructed with the standard iterative three-dimensional reconstruction method of the scanner.The tracer cloud is displayed as transparent isosurfaces, in order to show the intensity distribution.
The tracer distribution appears as smooth cloud, which expands up to 1 cm in a perpendicular direction from the fracture surface.This is due to the comparatively low spatial resolution of roughly 5 mm of this particular scanner and reconstruction method.However, it is obvious that the effective volume for the flow process is reduced to a preferential pathway along the fracture plane.A later experiment with a slower flow rate of 0.001 mL min −1 and with the longer living PET nuclide 124 I yielded some broadening of the flow pattern within the fracture plane.The spatial resolution was not sufficient to confirm dispersion of the tracer into the material adjacent to the fracture.
Practically, the conditions in a nuclear medical department are suitable for short-term experiments with short-living nuclides, as far as practicable with a mobile experimental setup.They are inappropriate for long-lasting procedures, specific analytical online techniques, and application of long-  living radionuclides.Additionally, the non-transparent proprietary reconstruction software could restrain adjustment to the conditions of strongly absorbing and scattering dense material.In our example, we had no control of the absorption and scatter correction procedure, which might have caused unnecessary blurring of the images.Consequently, since then we use a high-resolution scanner with which we produced the following examples.

Staßfurt sandstone (high-resolution scanner)
A complete sandstone drill core (diameter 100 mm) from the bunter sandstone in Staßfurt (NE Germany) with a diagonal fracture system parallel to the bedding was cast in epoxy (Fig. 7).The fractures were sealed with adhesive tape to prevent the resin entering.The core was cut at the inlet side in order to prepare a smooth surface at the inlet with a minimum gap space(< 0.1 mm) between the sample and the end cap because preceding experiments had shown that the active tracer in a large gap volume acts as a persistent source impairing image quality by scattered radiation.
Saline formation water was prepared and flushed through the sample with a flow rate of 0.02 mL min −1 for 1 day.5 mL of this carrier solution was labelled with 0.01 M [ 18 F]KF with an activity of 150 MBq and injected with the same flow rate.Then, after 250 min, unlabelled carrier solution was injected again.The experiment is reported in detail in Wolf (2010).
The input zone of the sample was positioned during and after the injection of the tracer pulse in the FOV of the PET scanner.38 frames were recorded during 10 h of continuous injection, with a frame length increasing from 4 to 30 min.The increasing frame length accounts for the tracer decay.
In Fig. 8, the isosurfaces at the 3 % maximum value of six frames are shown.The frame time, length, and injected volume are given in Table 3.
The tracer propagation pattern appears to follow selected pathways across the fracture surface of the single fracture  7; isosurfaces of the 3 % maximum value for six frames (see Table 3).that intersects the sample cross section.These pathways form a complex network of diverging and recombining patterns.At the zone where the pattern reaches the sample surface, 50 mm from the injection plane, the tracer disappears.This is both caused by the decay and the dispersion into a larger volume because both effects reduce the activity concentration and thus decimate the number of counts per voxel below the noise level and the detection threshold.The tomograms of Fig. 8 were produced with an enhanced version of the original reconstruction procedure.Corrections for random coincidences, dead time, attenuation, scatter, and image normalization were performed.In particular, the scatter correction strongly enhances image quality.The scattered radiation was determined using the STIR SSS algorithm.These scatter data were scaled to fit the scatter fraction of 45 %, which was estimated from Monte Carlo simulations on similar material.As an example, the impact of the scatter correction is shown in Fig. 9 as cross sections of frame 5 from Fig. 8.A spatial amplitude distribution of the same frame is given in Fig. 10.
The relative stochastic error is considered according to Eq. ( 16).The variance was calculated from the backprojection of the uncorrected data, and the reference tr(v) was es-  timated from the backprojection of the random-and scattercorrected data.Figure 11 shows that the maximum error is 10 % in the responsive regions where significant numbers of counts were recorded.It depends on the geometric sensitivity, which is included as the backprojection coefficient, and the total number of counts per voxel.The major portion of the stochastic error in these responsive regions is caused by scatter, while the typical random rate is about 1 %.Outside of these regions, the random rate, and thus the error and SNR, is due to background radiation, scattered events, and coinciding singles from different decay events.In zones adjacent to active regions, the impact of scattering degrades the SNR and the detection threshold.Far from active regions, scattering becomes negligible, and the detection threshold only depends on the low random rate.
Figure 11 also shows the impact of less sensitive zones on the stochastic error.In these zones, the total count rate is reduced because they respond to a smaller number of pro-  16), and isosurface of the 3 % maximum value, computed for frame 5 of Fig. 8. jection bins, and therefore the error increases.Unfortunately, currently the OSMAPOSL reconstruction algorithm preferentially projects inconsistent events (i.e.noise) into these less sensitive zones.This causes additional systematical errors which are recognizable as circular artefacts.Here, we could eliminate these artefacts to a large extent by optimizing the normalization data, although they are still perceptible in the cross sections of Fig. 9 as stripe patterns parallel to the axis.

Diffusion in Opalinus Clay
Imaging with PET is a particularly convenient method when diffusion processes are considered, where the only observable is concentration of a chemical species.We established PET as a quantitative method for observing spatially resolved diffusion patterns.From these patterns, we derive anisotropic diffusion coefficients and information on heterogeneity on the scale of drill cores.
With measurements and Monte Carlo model simulations, we found that by its long half-life of 2.603 years and its well-characterized chemical properties, 22 Na is a particularly suited tracer for determining the material impact on tracer diffusion.The propagation velocity of the tracer front is on the order of 0.5-2 mm d −1 (Kulenkampff et al., 2016), which Clay, after 13 and 27 days.Bottom: axial projections with fitted FWHM ellipsoids (full width half maximum).Deviations from the elliptical shape indicate heterogeneous effects that are not considered as homogeneous anisotropy (Kulenkampff et al., 2016).can be derived from the elliptical spreading of the tracer cloud (Fig. 12).Thus, this effect is hardly observable with clinical PET scanners with a resolution above 3 mm.Anisotropic diffusion coefficients have been derived using an optimized finite element model (Lippmann-Pipke et al., 2016).These were in accordance with laboratory results of through-diffusion experiments.Departures from elliptical spreading are indications of heterogeneous effects that could be caused by fine layering.

Conclusions
GeoPET is proven to be applicable for tomographic process monitoring in geological materials.In fact, its molecular sensitivity, excellent selectivity for the decay of the particular radionuclide, and comparably high robustness with regard to material effects make it an outstanding option for quantitative tomography of tracer transport, which is unrivalled by other tomographic modalities.The method should be applied when heterogeneous processes are considered.
It is highly recommendable to apply high-resolution scanners, like our ClearPET, because these provide the maximum achievable spatial resolution and a higher sensitivity than large clinical scanners.This resolution on the order of 1 mm is the integration volume for the tracer concentration.It is a reasonable scale for process observations on the core scale, in the range of 10 cm, which also is the order of the representative elementary volume for a large group of inhomogeneous materials.The resolution of clinical PET scanners (3-5 mm) is rather poor, compared to the maximum sample diameter of 10 cm, and tends to equilibrate the major features of inhomogeneous transport patterns.
However, high-resolution PET scanners afford more workrelated expenses and are not yet mature technology.There is still need for further development, concerning reconstruction software, detection hardware, and the availability of tracers and labelling methods.
The progressing applications of geochemical simulations on the basis of pore-space models from µ-CT images call for experimental verification.It appears that PET is a unique method which is capable of providing such data.Apart from that, the examples demonstrate that it is able to deliver experimental process parameters (for example, flow path distribution and flow velocity) that are hardly derivable with other methods and that are of fundamental relevance for reactive transport.The availability of such parameters offers prospects on realistic intermediate-scale model simulations of reactive transport which inherently incorporate the effects of conservative flow.
PET is the potential gold standard for geoscientific transport process tomography in laboratory studies, as it is for functional studies in biomedical research and clinical applications.This is due to its unrivalled sensitivity for tracer concentrations in combination with an appropriate spatiotemporal resolution and immunity to matrix effects.These benefits are at the cost of high expenses and some limitations with respect to available tracers and labelling methods.

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.
) with the probability p(v, b) that one coincidence originating from the voxel v in the three-dimensional volume space V is detected and assigned to the bin b in the four-dimensional projection space B, which is divided into subsets B k .The detection probability matrix is composed of p n (b) : efficiency of the detector pair attributed to b, p ρ (b): probability that the event is not attenuated in the medium, and p g (v, b): geometrical probability, according to p (v, b) = p g (vb) × p n (b) × p ρ (b).

Figure 3 .
Figure 3. Normalization tomogram of ClearPET: backprojection of the measurement of a homogeneous cylinder source.

Figure 4 .
Figure 4. Reconstructed image of the Monte Carlo simulation of four point sources in Opalinus Clay material.

Figure 5 .
Figure 5. Peak amplitudes of the four point sources in Fig. 4 with dependence on the total activity.

Figure 6 .
Figure 6.Top: µ-CT image of a granite core with a single fracture (recorded at the Federal Institute for Material Research and Testing BAM, Berlin, and processed by F. Enzmann, Johannes Gutenberg University, Mainz).Below: six PET frames recorded with a clinical PET scanner at the Clinic and Polyclinic for Nuclear Medicine of Leipzig University Hospital.

Figure 7 .
Figure 7. Bunter sandstone drill core.Adhesive tape, preventing intrusion of the resin, indicates the direction of the fractures (image taken from Wolf, 2010).

Figure 8 .
Figure 8. Result of the measurements with ClearPET during advection experiments on the sample in Fig.7; isosurfaces of the 3 % maximum value for six frames (see Table3).

Figure 9 .
Figure 9. Vertical slice through frame 5 from Fig. 8; left: without scatter correction, right: with scatter correction.The colour scale refers to the maximum value of the tomogram.

Figure 11 .
Figure 11.Sections of the relative stochastic error according to Eq. (16), and isosurface of the 3 % maximum value, computed for frame 5 of Fig.8.

Figure 12 .
Figure12.Top: horizontal slices and axial projection of the tracer concentration tomograms during diffusion of 22 Na into Opalinus Clay, after 13 and 27 days.Bottom: axial projections with fitted FWHM ellipsoids (full width half maximum).Deviations from the elliptical shape indicate heterogeneous effects that are not considered as homogeneous anisotropy(Kulenkampff et al., 2016).