Freely available online through the author-supported open access option.

Abstract

Magnetic resonance imaging (MRI) was applied to the study of flow processes in model and natural soil cores. Flow velocities in soils are mostly too slow to be monitored directly by MRI flow velocity imaging. Therefore, we used for the first time diethylenetriaminepentaacetate in the form Gd-DTPA2− as a tracer in spin echo multislice imaging protocols with strong weighting of longitudinal relaxation time T1 for probing slow flow velocities in soils. Apart from its chemical stability, the main advantage of Gd-DTPA2− is the anionic net charge in neutral aqueous solution. We showed that this property hinders adsorption at soil mineral surfaces and therefore retardation. We found that Gd-DTPA2− is a very convenient conservative tracer for the investigation of flow processes in model and natural soil cores. With respect to the flow processes in the coaxial model soil column and the natural soil column, we observed totally different flow patterns. In the first case, the tracer plume moved quite homogeneously in the inner highly conductive core only and the migration into the outer fine material was very limited. A numerical forward simulation based on independently obtained parameters showed good agreement between experiment and simulation and thus proves the convenience of Gd-DTPA as a tracer in MRI for soil physical investigations. The natural soil core, in contrast, showed a flow pattern characterized by preferential paths, avoiding dense regions and preferring loose structures. In the case of the simpler model column, the local flow velocities were also calculated by applying a peak tracking algorithm.

Flow processes in soils are generally too slow to be monitored directly by magnetic resonance flow velocity imaging. This study showed the convenience of Gd-diethylenetriaminepentaacetate as a conservative tracer in soil cores and validates the results with a numerical, three-dimensional, forward simulation. The method allows direct observation of different flow patterns such as matrix and preferential flow.

Flow processes in soils control the water and nutrient supply in plants, contaminant transport, and the recharging of aquifers. Natural soils are heterogeneous media with locally varying grain size, pore size, and pore geometry. In addition, highly conductive structures like wormholes and decomposed roots are present. All these factors lead to so-called preferential pathways, along which water flows faster vertically than expected from the average flow velocity vav [(volume flow rate/area)/water content: vav = (Q/A)/θ = Jw/θ] but from which water can also diffuse laterally into neighboring small-pore-size regions. All these processes can strongly affect water redistribution in soils, its bioavailability for plants, and transport of contaminants (Cislerova et al., 1997; Clothier et al., 2008). Classical investigations of such processes have been invasive and not temporally resolved. For example, soil cores are infiltrated with dyes and afterward sliced and photographed, so that only the static result of flow processes can be observed (Wehrhahn et al., 2007). The direct observation of such processes in situ is difficult due to the opaque nature of the systems. Therefore, during the past two decades, tomographic methods have become popular for the investigation of soils and soil–root systems. For example, x-ray tomography (also termed computed tomography) has been used to develop images of the physical density of soil (Cislerova and Votrubova, 2002; Wildenschild et al., 2002). Magnetic resonance imaging has been especially useful because it directly probes the substance of interest: water.

Basics of Magnetic Resonance Imaging

The basics of MRI are treated extensively in textbooks and review articles. Therefore, only the most important features are mentioned here (Callaghan, 1991; Nitz, 1999; Nestle et al., 2002; Lascialfari and Corti, 2007). The basis of MRI is the nuclear magnetic resonance effect of atomic nuclei. The most abundant nucleus for imaging is the 1H nucleus (proton) in H2O, although other nuclei can also be used for imaging purposes. In a magnetic field B0 oriented in the z direction, proton spins exist in two states: precession parallel and antiparallel with respect to the direction of B0, with the Larmor frequency ν0= |B0H/2π, where γH is the gyromagnetic ratio of 1H. For an ensemble of spins in thermal equilibrium, the parallel state is slightly more strongly populated than the antiparallel state according to Boltzmann's law, resulting in macroscopic magnetization in the z direction This ensemble can be excited by the absorption of an electromagnetic radiation pulse with the frequency ν0, which is often termed the radiofrequency pulse (rf pulse) because for usual magnetic field strengths, ν0 is in the range of some tens to hundreds of megahertz. If the pulse is sufficiently long (90° pulse), the macroscopic magnetization is turned into the xy plane and subsequently the spin ensemble relaxes via two different processes: longitudinal (T1) and transversal (T2) relaxation. The latter can be subdivided into at least two components: relaxation due to dephasing in external magnetic field inhomogeneities and relaxation due to random dephasing. In MRI scanners, the first effect strongly dominates the overall relaxation, but this effect can be eliminated by the creation of a so-called spin echo, where the dephased magnetization in the xy plane is rephased by a 180° pulse after a time period τ, which creates an echo signal after an additional time τ that can be detected and further analyzed. The time between the original 90° pulse and the echo is known as the echo time, tE = 2τ. The system relaxes further back to thermal equilibrium with the characteristic time T1, however, so that after about 5T1, the full original magnetization is restored and a new cycle can be started. The time period between two subsequent excitations with a 90° pulse is termed the repetition time, tR. It should be noted that both tE and tR are experimental parameters that can be used to control the signal intensity and contrast in an MRI experiment.

Generally, the relaxation times T1 and T2 are characteristic for the system, and in porous media they depend on the physical and chemical environment of the water molecules, such as pore size, pore geometry, filling factor, and chemical composition of the pore walls and the fluid. Changes in one of these factors will therefore change the signal intensity.

Imaging is performed by additional switching of external magnetic field gradients in the three Cartesian coordinates. One of the most used imaging sequences is the spin-echo multislice (SEMS) method, where one slice is first selected by a slice-selective 90° rf pulse. The second and third dimensions are then addressed by the so-called phase- and frequency-encoding gradients. After the application of the 180° pulse, the created echo is recorded, during which a frequency-encoding gradient is reapplied. The whole sequence of exciting and rephasing rf pulses, slice, phase- and frequency-encoding gradients can be visualized in a pulse diagram. The image is obtained in a space with reciprocal length dimensions (k-space), from which the image in real space is reconstructed by multidimensional Fourier transformation. In addition to this basic spin-echo sequence, a huge variety of further sequences exist, partly for very special imaging purposes (Callaghan, 1991; Nitz, 1999; Blümich, 2000; Lascialfari and Corti, 2007) and also for the imaging of flow velocities. These methods are mostly inconvenient for investigating flow processes in soils, however, because the lower limit of these methods is on the order of some hundred micrometers per second, and flow velocities in soils are mostly slower. Second, the major obstacle for using MRI in soils is the inherent fast transverse relaxation, which is in the range between 0.1 and 10 ms (Hall et al., 1997; Votrubova et al., 2000; Haber-Pohlmeier et al., 2010). The reasons for this are the presence of paramagnetic impurities and small pore sizes, because T2 relaxation is caused by the dephasing of local magnetic moments in local magnetic field gradients, which are created at solid–liquid interfaces.

Thus, to investigate flow processes in soils, the use of tracers is favorable to modify the MRI signals appropriately. One strategy is to reduce the spin density by replacing 1H2O with 2D2O, which does not resonate at the resonance frequency of ordinary water (Pohlmeier et al., 2008). This requires relatively high concentrations of D2O, which enhance the physical density of water and may be toxic for plants. A second strategy is the use of paramagnetic compounds, which mainly reduce the longitudinal relaxation time of water based on the molecular mechanism of the exchange of water molecules in the hydration shell of the ion (McConnell, 1987; Hermann et al., 2008). The signal intensity in the spin-echo MRI sequence that was used in this study is given by 
\begin{eqnarray*}&&S(x,y,z)(){=}S_{o}(x,y,z)\left[1{-}\mathrm{exp}\left({-}\frac{t_{R}}{T_{1}(x,y,z)}\right)\right]\\&&{\times}\mathrm{exp}\left({-}\frac{n_{E}t_{E}}{T_{2}(x,y,z)}\right)\end{eqnarray*}
[1]
where S is the measured signal intensity as a function of space (x,y,z), S0 is the extrapolated signal intensity for tE → 0 and for tRT1, which is proportional to the proton spin density, T1 and T2 are the effective longitudinal and transversal relaxation times, and tR, nE, and tE are experimentally adjustable parameters (repetition time, number of echoes, and echo time, respectively) (Callaghan, 1991; Blümich, 2000). Therefore, if T1 is locally reduced by the presence of a paramagnetic contrast agent, the local intensity and the contrast can be enhanced by selecting a short tR. The functional dependence of voxel intensity on tR is shown schematically in Fig. 1 for two cases: the solid line represents the signal intensity in porous material filled with pure water and the dashed line in the presence of a contrast agent (tracer). For sufficiently short tR, the contrast between regions with and without tracer is optimal, as indicated by the vertical line.

In soils with small values of T2, caused by small pores, paramagnetic ions (Hall et al., 1997), and high susceptibility contrasts (Bray and Hornak, 2009), tE should be set as small as possible. Obviously, there are technical limits due to the necessary switching of gradients. It should be noted here (although it was not a subject of this study) that reducing tE is a major task to enable the future use of MRI for soil systems. Because medical scanners cannot allow short pulses, dedicated scanners will be required with robust, strong gradient systems and ultra-short-pulse sequences will need to be developed and adapted (Robson et al., 2003; Marica et al., 2006; Pohlmeier et al., 2007).

The most popular tracers (contrast agents) for MRI in porous media are dissolved paramagnetic metal cations, such as Cu2+(aq), Ni2+(aq), or Gd3+(aq). As mentioned above, the advantage of solvated metal ions is a very effective T1 relaxation enhancement due to the fast exchange of water in the inner solvation shell. For example, Greiner et al. (1997) and Oswald et al. (1997) used CuSO4 to monitor fluxes in artificial porous media. Later, this investigation was extended to the motion of salt water clouds in porous media (Oswald et al., 2007) and the motion of Ni2+ ions in the rhizosphere (Moradi et al., 2008a) as well as in model porous media (Moradi et al., 2008b). The quite fast T2 relaxation was also the reason why these researchers failed to obtain a signal from a natural loamy soil. A further interesting experiment was performed, in which the motion of a Mn2+ tracer cloud was monitored in a macroscopic heterogeneous system consisting of a stochastic arrangement of blocks of different sands. Because the transverse relaxation was sufficiently long, the motion could be observed at quite long tE = 10 ms. The cloud mainly moved through preferential pathways created by coarse regions in the center of the model aquifer (Yoon et al., 2008).

One general drawback to the use of solvated metal ions as tracers is their tendency toward adsorption or surface precipitation as hydroxyl complexes. This can be partly prevented by acidifying the solution, but the risk of retention of the tracer in the medium should generally be checked (Jelínková et al., 2008).

It would be more convenient to use bulky neutral or anionic contrast agents, such as chelated metal ions, which are used in medical MRI (Hermann et al., 2008). An early example for porous media is the inflow of a Mn2EDTA solution into a glass bead pack (Guillot et al., 1991). A common medical contrast agent is Gd-DTPA2−, where the extremely strong paramagnetic Gd3+ ion (seven unpaired electrons in a half-filled f-shell) is complexed by three amino and five carboxylate groups of the strong chelator DTPA5−. The ninth position in the inner hydration sphere of the central Gd3+ ion is occupied by a water molecule, which can easily be exchanged by pore water and thus affects the enhanced T1 relaxation. The huge advantage of this contrast agent is its chemical stability across a wide pH range and its net charge of −2 in neutral solution. Therefore, in contrast to bare metal cations, adsorption is not expected at mineral surfaces.

The aim of this study was to elucidate whether Gd-DTPA2− is convenient for the monitoring of flux processes by MRI. In a preliminary study, the relaxation time of a medium sand saturated with an aqueous solution of Gd-DTPA2− was investigated to obtain optimal concentration ranges for the best contrast. For the further investigation of fluxes, we chose two model systems: (i) a coaxial model sample composed of an inner core of medium sand with high conductivity and an outer core of a sand–silt mixture with a much lower conductivity, for which experimental observations were compared with modeling results to evaluate their consistency with theoryin order to answer the question of whether water moved only through the central core or whether it diffused into the finer porous compartment; and (ii) a natural soil core of a loamy sand, of which the T1 and T2 relaxation time was known (T1 = 80 ms, T2 < 2 ms [Haber-Pohlmeier et al., 2010]) to answer the questions of whether the tracer motion could principally be observed, whether it was retarded, and which pathways were characteristic for preferential flow. Furthermore, batch equilibrium experiments were performed to check whether or not Gd-DTPA2− adsorbed to the materials used.

Materials and Methods

Soil Columns

For the tracer-infiltration experiments, two different soil columns were used: a coaxial cylindrical column with an inner diameter of 8 cm consisting of a core of FH31 medium sand surrounded by a Millisil W3 sand–silt mixture (Quarzwerke GmbH, Frechen, Germany) (see Fig. 2 ), and a soil column with an inner diameter of 3 cm with natural, not repacked, sandy loam soil. The coaxial sample was prepared as follows. First, the inner core was produced by filling a plastic foil tube with FH31 sand, which was saturated to about 90% and frozen. Second, the frozen core was unwrapped and placed on the bottom bed of sand in the column (Fig. 2, left), and the outer core was filled from the top. The whole system was finally saturated from the bottom and closed at the top by a porous glass filter plate (RoBu Glasfilter-Geräte GmbH, Hattert, Germany).With respect to the soil texture, FH31 is a natural medium sand from Frechen, Germany, that underwent pretreatment including washing and fractionating to the desired grain size. The W3 material is a quartz powder produced from quartz sand that was artificially milled in ball mills and subsequently fractionated. The origin of the natural sandy loam soil was the Ap horizon of a Gleyic Cambisol from Kaldenkirchen, Germany. Table 1 summarizes the most important physical and chemical characteristics of the materials. To check whether the chosen tracer, Gd-DTPA2−, adsorbs at the soil material, batch equilibrium experiments were performed at a solid concentration of 50% (w/w) and Gd concentrations between 0 and 2.5 mmol L−1 at a temperature of 20°C and pH of 6.5. The suspensions were shaken for 24 h, then centrifuged, and the supernatant was analyzed for Gd by means of inductively coupled plasma mass spectroscopy. No adsorption of Gd-DTPA2− was observed.

Flow Experiments

The columns were irrigated with tap water under steady-state conditions by dripping the water onto the center of the top porous plate, which was exposed to atmospheric pressure. The lower boundary condition was a seepage face. Volume fluxes of 3 and 0.5 mL min−1 were applied to the coaxial sample and the soil core, respectively. After achieving steady state, which was checked by comparing outflow and inflow rates, the systems were equilibrated further for 3 h so that, due to the lower boundary, the model soil columns were near saturation (Table 1). The tracer experiments were started by injecting pulses of 2 (coaxial model sample) and 1 mL (natural core) of Gd-DTPA2− solution (concentration cGd = 2 mmol L−1) into the irrigation flux 50 cm above the outlet while keeping the overall flux rate constant. The pulse reached the outlet at time point = 0. From this time point on, the position of the tracer cloud was monitored by MRI until the plume left the field of view.

Magnetic Resonance Imaging Measurements

In this study, the MRI experiments were performed on two scanners: an ultra-wide-bore vertical 4.7-T Varian system (Magnex Scientific, Oxford, UK) with a gradient system of 300 mT m−1 and a 170-mm rf birdcage coil (RAPID Biomedical GmbH, Würzburg, Germany), and a 7-T MRI scanner from Oxford, UK, equipped with a Bruker gradient system of 300 mT m−1 and a 38-mm rf birdcage (Bruker Biospin, Ettlingen, Germany). Both scanners were operated by VNMRJ software running on a Varian console.

The nuclear magnetic resonance (NMR) images of the coaxial model sample were recorded by a SEMS sequence with tE = 4.8 ms, tR = 200 ms, field of view = 128(x)128(z), where z is the vertical direction, 17 by 17 cm, 1.33 mm pixel−1, isotropic, slice thickness = 2 mm. To achieve the best contrast between the pixels with and without tracer, measurements were performed on a model sample consisting of 11 tubes bound together in a ring structure. They were filled with the FH31 medium sand saturated with an aqueous solution of Gd-DTPA2− in a concentration range from 0.5 to 50 mmol L−1 (see Fig. 3 ). The measurements were performed at of tR = 0.2, 0.4, and 1.0 s. Figure 3 shows a schematic picture and axial images at different tR values.

To image the natural soil core, it was necessary to use the shortest possible echo time because the T2 of this soil is very short. Simultaneously, the overall measuring time was also kept short to allow us to monitor fast preferential flow via the Gd-DTPA2− tracer. For the natural soil core, we therefore used a fast SEMS sequence with tE = 1.9 ms and a turbo factor of two (Haacke et al., 1999). With a field of view of 60 by 40 mm at a matrix size of 64 by 64 and a slice thickness of 1 mm, the water transport could be monitored via the Gd-DTPA2− tracer. The repetition time was set to 0.05 s to obtain some saturation because Kaldenkirchen soil already had a T1 of 80 ms. The temporal spacing between single images was 1 min.

Numerical Simulation

The experiment on the coaxial model sample was numerically simulated. We assumed that the water flow in the three-dimensional porous sample under the given conditions could be described locally by the Richards equation: 
\[\frac{{\partial}\mathrm{{\theta}}(h)}{{\partial}t}{=}{\nabla}{\cdot}[K(\mathrm{{\theta}}){\nabla}h]{-}\frac{{\partial}K(\mathrm{{\theta}})}{{\partial}_{z}}\]
[2]
where t is time [T], h is pressure head [L], θ is volumetric water content (dimensionless), K is hydraulic conductivity [L T−1], and z is the vertical coordinate directed downward [L]. Both θ(h) and K(θ) are highly nonlinear functions, here described by the van Genuchten–Mualem parametric expressions (van Genuchten, 1980). The hydraulic parameters of the two materials FH31 and W3 were separately determined using standard laboratory measurements (pressure plate apparatus and multistep outflow experiments) and are given in Table 1; the parameters of the porous plate were taken from Pohlmeier et al. (2008). Richards' equation was numerically solved with a cell-centered finite-volume code developed by O. Ippisch, IWR Heidelberg (e.g., applied in Ippisch et al., 2006; Rossi et al., 2008). The code provides numerical stability and divergence-free flow fields, even for highly heterogeneous porous media with sharp material discontinuities.
For the solute transport simulation, we assumed that solute transport is governed by the advection–dispersion equation: 
\[\frac{{\partial}\mathrm{{\theta}}C_{\mathrm{w}}}{{\partial}t}{=}{\nabla}{\cdot}[\mathrm{{\theta}}\mathbf{\mathrm{u}}C_{\mathrm{w}}{+}\mathrm{{\theta}}\mathbf{\mathrm{D}}{\nabla}C_{\mathrm{w}}]\]
[3]
where Cw is solute concentration in the water phase [M L−3], u is the pore water velocity vector [L T−1], and D is the local-scale dispersion coefficient tensor (L2T−1), usually defined as 
\[\mathbf{\mathrm{D}}{=}(\mathrm{{\alpha}}_{\mathrm{T}}{\vert}\mathbf{\mathrm{u}}{\vert}{+}D_{\mathrm{m},\mathrm{eff}})\mathbf{\mathrm{I}}{+}(\mathrm{{\alpha}}_{\mathrm{L}}{-}\mathrm{{\alpha}}_{\mathrm{T}})\frac{\mathbf{\mathrm{uu}}^{\mathrm{T}}}{{\vert}\mathbf{\mathrm{u}}{\vert}}\]
[4]
where ‖u‖ = √(uuT) is the Euclidean norm of the pore water velocity vector, αL and αT are the longitudinal and transverse dispersivities, respectively, Dm,eff is the effective molecular diffusion coefficient, and I is the identity matrix. We set αL to the mean grain size and αT to 1/10 of αL, which are typical values reported for saturated unconsolidated homogenous porous media (e.g., (Greiner et al., 1997; Yoon et al., 2008; Lamy et al., 2009)). The diffusion coefficient Dm of Gd-DTPA2− in water at 25°C was assumed to be 0.4 × 10−9 m2 s−1 (Osuga and Han, 2004) and was multiplied by the tortuosity factor evaluated as a function of water content using the relationship of Millington and Quirk (1961) to obtain Dm,eff.

The flow model was coupled to the particle tracking code PARTRACE (Agrosphere Institute, Jülich, Germany), which is based on the equivalent stochastic differential version of the advection–dispersion equation and uses an improved algorithm accurately accounting for discontinuous dispersion tensors (Bechtold et al., unpublished data, 2010).

Water flow and solute transport was simulated in three dimensions by a grid of 80 by 80 by 137 cells with a regular grid spacing of 0.1 cm. The numerical domain was divided into five regions, i.e., the sand, the sand–silt mixture, the porous plates at the top and bottom, and a zero-conductivity region around the soil column (see setup in Fig. 2). Irrigation was uniform and constant in time applied to a circular region (diameter = 2 cm) at the upper boundary with Jw = 3.0 mL min−1, which mimicked the single-drop irrigation used in the experiment. The lower boundary condition was the seepage face. The tracer was applied under steady-state flow conditions as in the experiment by setting the concentration of the irrigating water for 40 s to 2 mmol L−1.

Results and Discussion

Relaxation Times

First, we checked the relaxation times of the Gd-DTPA2− solutions in the porous medium. For this purpose, the intensity of the SEMS sequence was investigated for different repetition times, echo times, and concentrations of Gd-DTPA2−. Figure 3 shows the results obtained for the FH31 medium sand for one axial slice through a model sample consisting of a ring of sand-filled cuvettes saturated with Gd-DTPA2− solutions of different cGd. The experiments were performed for different repetition times 0.1 s < tR < 1 s at the shortest possible echo time tE = 4.8 ms. The signal intensities were normalized to the maximum signal intensity S0 and are plotted in Fig. 4 as a function of cGd. Equation [1] was fitted to the data set, where the relaxation rates 1/T1 and 1/T2 were given by (Haacke et al., 1999) 
\[\frac{1}{T_{1}}{=}\frac{1}{T_{1,\mathrm{eff}}}{+}a_{1}c_{\mathrm{Gd}}\]
[5a]
 
\[\frac{1}{T_{2}}{=}\frac{1}{T_{2,\mathrm{eff}}}{+}a_{2}c_{\mathrm{Gd}}\]
[5b]
where T1,eff and T2,eff are the effective relaxation times of water in the porous medium including bulk and surface relaxation. These quantities are independent of the tracer concentration and were determined by separate infrared Carr, Purcell, Meiboom, Gill (CPMG) measurements using tE = 4.8 ms (1/T1,eff = 0.9 s−1, 1/T2,eff = 6.6 s−1). By fitting Eq. [1] and Eq. [5a] and [5b] to the data set (fitting parameters a1,a1, and S0), we obtained a1 = 4 L mmol−1 s−1, a2 = 10 L mmol−1 s−1, and S0 = 0.7, which underlines the clear dominance of the observed signal intensity S by the contrast agent and the minor importance of the influence of longitudinal surface and bulk relaxation. At Gd concentrations >5 mmol L−1, however, the S/S0 curves decreased in a similar manner, meaning that tR did not influence the signal intensity, i.e., the T1 weighting was lost and the signal intensity was T2 weighted. Second, for Gd concentrations <3 mmol L−1 and for tR ≤ 0.2 s, a linear dependence of S/S0 on log(cGd) was expected. As a result, in our experiments, a cGd of 2 mmol L−1 was used, which is in the linearly increasing part of the curves.

Tracer in the Coaxial Model Sample

The next experiment involved the investigation of tracer transport in the coaxial model sample composed of an inner core of medium sand and an outer core of sand–silt mixture (see Fig. 1). Figure 5 shows the motion of the tracer plume through the system immediately before injection (top left) and 31 images, each separated by 37.2 s. In the image at t = 0, we can differentiate between the inner fine-sand core and the outer sand–silt mixture, which shows slightly less signal intensity. This is due to faster relaxation in the smaller pores of this medium, although the total porosity of ε = 0.41 cm3 cm−3 is comparable to that of the medium sand (ε = 0.39 cm3 cm−3). The homogeneous signal intensity within each material indicated that the system was near saturation and in equilibrium under the given conditions (steady-state flow and a free-drainage lower boundary). Furthermore, in the bottom part consisting of medium sand, an anomaly is visible, probably due to packing at a different density or a slight fracture. This structure developed during the initial accommodation period. At t = 37.2 s (second image in the first row), the tracer plume entered the system, and the plume started to infiltrate. Tracer can be seen to proceed through only the inner core; no migration of the tracer into the outer core was observed under these experimental conditions. Approximately at 9 min, 18 s after injection (time point 15), the plume reached the lower end of the inner core and started to spread laterally into the bottom of the system, which consisted of medium sand. After time point 23, it also infiltrated through the fractured zone and was eventually strongly diluted and hardly visible in the bottom region. The average velocity of the tracer plume was 1.2 cm min−1, nearly identical to the average pore flow velocity of 1.09 cm min−1 calculated from the irrigation rate of 3 cm3 min−1 and assuming that it only moved into the inner core of 3-cm diameter. A significant tailing was not observed. The lack of a pronounced tailing and the consistency between the measured and calculated average pore velocity indicate that the exchange of the tracer with the outer core was very limited. Another characteristic feature of the migrating tracer plume is the change from a convex to a rather concave shape at around time point 10.

A closer look at the images in Fig. 5 reveals some tailing at the interface between the inner and outer cores; however, some leading is also visible in the outer region of the inner core. These contradictory observations could be caused by the preparation procedure. The inner core was prepared by filling a plastic tube with sand, saturating it to about 90%, and freezing it. Because freezing starts from the outside and works its way inside, this could have produced some larger pore sizes in the outer region of the inner core. As a result, the hydraulic conductivity was slightly higher there, causing the leading effect observed in the lower half of the model system. A possible explanation of the slight tailing is that small amounts of tracer diffused laterally across a small distance into the outer core at small concentrations. After the main plume had passed, this tracer moved back and was transported at the interface to the bottom. To check this hypothesis, we performed an analogous experiment for W3 as shown for FH31 in Fig. 3 and 4. The result was that tracer should be detectable at a concentration of 0.3 mmol L−1 if a detection limit of S/S0 > 0.2 is assumed. In the experiments presented here, however, no signal in the outer core was visible, i.e., a higher concentration of tracer was not present. The reason is probably the low hydraulic conductivity of the sand–silt mixture, preventing migration of larger amounts of tracer into the outer core.

Flow Velocities

To calculate the flow velocities in the z direction, an extended Interactive Data Language routine based on the front-track algorithm (Herrmann et al., 2002) was used. The originally noisy data were fitted by a Gaussian function along the vertical axis (z axis) for all time points: 
\begin{eqnarray*}&&S_{0}(t_{i}){=}A_{0}(t_{i})\mathrm{exp}\left\{{-}\frac{1}{2}\left[\frac{z{-}c(t_{i})}{w(t_{i})}\right]^{2}\right\}\\&&{+}A_{1}{+}A_{2}z{+}A_{3}z^{2}\end{eqnarray*}
[6]
where S0(ti) is a fitted NMR signal at time point ti, A0(t) is the amplitude at time point ti, z is the height from the bottom, c(ti) is the vertical position of the center at time point ti, w is the width, and A1, A2, and A3 are correction coefficients for vertical background inhomogeneities caused by systematic rf-field intensity variation of the MRI probe head. To illustrate the fitting process, Fig. 6 shows the fit results for the central vertical profile at three time points: ti = 3 min, 43 s; 7 min, 49 s; and 10 min, 32 s after the start. The returned 128 by 128 array contains the fitted images and the positions of the peak centers c(ti) in Eq. [6]. Figure 7 (top) displays these fitted images S0(ti) for the first 32 time points (corresponding to ti = 0 s to 19 min, 13 s after injection). For easier recognition, a shape of the coaxial model sample has been drawn in the center by a wire frame. In the original MRI images (Fig. 5), it can be clearly seen that the flow took place initially in the inner core only; after the plume reached the bottom part, it spread out laterally. This is reliably represented in the fitted images. For the calculation of the z components of the flow velocities, the following equation was used: 
\[\mathrm{{\upsilon}}_{z}{=}\frac{{\Delta}z}{{\Delta}t}\]
[7]
with 
\[{\Delta}z{=}c(t_{i}){-}c(t_{i}{-}1)\]
and 
\[{\Delta}t{=}t_{i}{-}t_{i{-}1}\]
where c(ti) is the center of the fitted peak at the ith time point. The resulting velocity was then assigned to the time point ti, leading to an array with the original dimensions of the central slice, which contains the velocities at the vertical positions of the plume (Fig. 8 ). The curved lines correspond to the positions of the center of the tracer plume at different time points. The obtained velocity field therefore represents dynamics under steady-state conditions. We can differentiate among four regions (see also Fig. 7): (i) after entering the inner core, the flow velocities exhibited regular behavior with a curved shape; (ii) in the middle central core, the velocities slowed down on the left side, leading to a slight local retardation of the plume, while on the right side, a region with flow velocities higher than the average velocity are visible; (iii) this changed in the lower part, where the shape no longer bends downward but rather bends upward; and (iv) reaching the bottom part of the column, the velocities decreased strongly until they could not be further resolved. The velocity field image thus reflects the shape of the tracer plume satisfactorily. The average flow velocities correspond to the green part of the color bar, corresponding to an average flow velocity of about 1.1 cm min−1. This agrees well with the average flow velocity of 1.2 cm min−1 calculated from the irrigation rate and the porosity.

Comparison with Modeled Tracer Transport

The tracer experiment in the coaxial model sample was simulated with the three-dimensional numerical model described above to verify the experimental observations. The numerical model could be used because both the structure and the soil properties were well defined, which allowed the application of a distributed model based on Richards' equation and the advection–dispersion equation. The distribution maps of the modeled concentration are plotted at the same time points as for the experiment (Fig. 7, bottom). The color scale (blue to red, 0.05 to 1.0 mmol L−1) was kept constant throughout the sequence. Modeled concentration values that were <0.05 mmol L−1 are shown in black.

The principle features of the tracer displacement observed in the experiment are consistent with the results of the numerical simulation. The tracer was mainly transported in the inner core of the column and very little tracer mass was transported into the outer sand–silt mixture. As in the experiment, a thin line of tailing occurred at the sand and sand–silt mixture interface that is obviously not any artifact but explainable by the advection–dispersion equation. The concentration values penetrating deeper into the outer core were <0.05 mmol L−1 and far below the limit of detection of MRI in this study. At time point 16, the modeled tracer entered into the bottom region of the column fully composed of sand. The modeled data indicate that the dominant convex shape of the plume front present in the upper part was less pronounced in the lower part of the inner core. Before the tracer entered into the bottom region, the plume front approached a flat shape. This change in shape was caused by the hydraulics in the lower part. The cross-sectional area of the highly conductive sand was abruptly increased at the bottom of the inner core, which led to convex-shaped pressure head contour lines and thus to higher flow velocities at the outer part of the inner core relative to the center of the inner core. This flattened the convex shape of the plume front. The development of a concave plume front, as observed in the experiment, however, was not observed in the numerical simulation. This probably indicates additional packing effects, e.g., slightly lower bulk density of the sand in the outer part of the inner core due to sequential freezing or sporadic preferential pathways at the sand and sand–silt mixture interface. Both could cause the development of a concave shape of the plume front. In the bottom region, the model and experiment were consistent again, both indicating a pronounced convex shape of the plume front that was spreading laterally into sand.

Summarizing, it should be emphasized that the simulation is a forward simulation, i.e., it is based on independently obtained soil physical parameters. The observed slight deviations between the experiment and simulation were only due to inhomogeneities in the packing, which were not implemented in the simulation.

Tracer in a Natural Soil Core

The third experiment investigated the tracer motion in a natural soil core. The experiment was run in a manner similar to the second experiment. The soil core was irrigated at a constant flow rate, into which, at time point 0, a pulse of Gd-DTPA2− solution was injected. Figure 9 shows the tracer plume between 0 (top left) and 10 min after injection (bottom right). The average pore flow velocity calculated from the irrigation rate, the cross-sectional area, and the total porosity are indicated by the displacement of the yellow bar on the left-hand side of Fig. 6. At time zero, an S-shaped bright structure and a further area of high intensity are clearly visible at the left wall. Dark areas also appear that are neither stones nor open voids, but rather dense material. In these regions, T2 was less than tE (data from T2 maps not shown here), which means that the greater part of the signal was already relaxed at the time when the echo was created.

One minute after injection, the first tracer already appeared in the central S-shaped structure; after 2 to 3 min, the whole plume had entered the core. Most of the tracer moved through the left part of the core near the wall, while the rest entered the right half within 5 to 8 min. It is highly likely that the tracer not only moved through the monitored central slice but also entered neighboring regions. Eventually, the tracer completely left the field of view at the bottom. The images at t = 0 and t = 10 min are identical, i.e., no tracer remained in the monitored area of the core. The average velocity was about three to four times faster than the average pore flow velocity. The observed flow pattern is therefore considered to be preferential flow.

Conclusions and Outlook

This work has shown that Gd-DTPA2− behaves conservatively because no retention of the contrast agent was observed in the model sand or in the natural soil core. The reason for this is its anionic nature and bulky structure. Adsorption at negatively charged soil particle surfaces is therefore highly unlikely, as is adsorption via electrostatic attraction at positive sites due to its large size. Furthermore, anion exchange vs. OH ions, as observed for phosphates, is also unlikely. In agreement with these considerations, we found no adsorption of Gd-DTPA2− on the FH31 medium sand, the W3 fine sand–silt mixture, or the Kaldenkirchen soil in the range between 0 and 3 mmol L−1 at pH = 6.5. This tracer could thus be suitable for imaging further flow processes in natural soils by MRI.

The minimum detectable concentration can be derived from Fig. 4. It depends on the choice of the repetition time tR and the echo time tE, both of which should be as short as possible. There are machine-dependent limits, however, meaning that such information cannot be generalized. In our case, for tR = 0.2 s and tE = 4.8 ms, the lower limit was about 0.3 mmol L−1. Shorter values for tE can be achieved by increasing the slice thickness and the power of the exciting and rephasing pulses. Shorter values for tR can be achieved by decreasing the number of slices and the resolution in the phase-encoding direction.

The flow behavior in the coaxial model soil core was detected only in the inner core without entering the outer core, which consisted of much finer material with an equal entire porosity but smaller pore sizes. The low conductivity hindered significant transport of tracer into this material. A slight tailing was observed, however, which could be caused by diffusion of small amounts of tracer into the first layer of pores in the outer core. In this case, the local concentration would be below the detection limit under the given experimental conditions, and this layer could be too thin to be resolved in the given voxel size of 1.3 mm. After the main plume had passed, tailing was seen as cleaning out that area. The quantification by calculating flow velocities corresponded well with the regional MRI images and is a further step in the direction of a detailed description of flow processes in soils. The heterogeneous structure in the middle of the inner core, visible in the MRI images at higher intensity, can be recognized in the flow velocity image by slower local velocities. Much more pronounced is the deceleration of the flow when the plume reached the bottom, where it could spread laterally, as shown in the original MRI images, because the outer core of fine material was no longer present.

We further demonstrated by numerical simulation of the flow and transport that the MRI-derived details of the tracer displacement were consistent with flow and transport theory, i.e., Richards' equation and the advection–dispersion equation. The agreement between this forward simulation based on independently obtained soil physical parameters and the experiment is very high. The slight deviations in the flow pattern between the experiment and the model can be explained and were due to inhomogeneities occurred during packing, which was not considered in the numerical model. In future flow scenarios, interruptions and changes in irrigation rates should also be introduced into the experimental setup.

In contrast to the coaxial model soil column, the flow pattern in the natural core showed considerable contributions of preferential flow. Dense areas that appear dark in the MRI images were avoided by the plume, whereas loose regions appeared to be strongly conductive with locally high flow velocities. In a more microscopic view, these flow behaviors appear to be similar to that of the coaxial model column: the plume avoided the dense areas and followed loose structures. After the first positive results of the peak tracking in the simple coaxial model core, which possessed only vertical structures, the next step should focus on further quantifying flow velocities to describe the more complex behavior with respect to strong lateral components. Because these components will appear in three dimensions, improvements with respect to multislice or three-dimensional imaging are necessary. A second improvement is necessary concerning the echo time. Stronger gradients and more robust rf coils will allow shorter pulse lengths and therefore higher intensity.

We would like to thank the German Research Foundation (DFG) for financial support (Sta 511/4-1, PO 746/2-1, Transregional SFB TR32). We further thank Horst Hardelauf for support when coupling the finite-volume flow model to the PARTRACE code, and Ms. Wettengl and Ms. Lippert, Forschungszentrum Jülich, for the chemical analysis of Gd in the used soil material.