Unsaturated flow phenomena, such as unstable wetting fronts and preferential flow, cannot be investigated using small-scale sampling. Dye tracer experiments can help visualize the dynamics of water flow but are destructive and therefore irreproducible. We investigated the applicability of high-resolution ground penetrating radar (GPR) for nondestructive visualization of unsaturated flow patterns arising from a forced infiltration experiment. Synthetic studies using a reflection GPR two-dimensional finite difference time domain modeling code indicated that differences in water content caused by preferential flow and fingering could be resolved. Moisture content contrasts down to approximately 2.5% within the top 2 m were detectable, but with increasing degrees of heterogeneity in the subsurface it becomes difficult to distinguish these moisture content changes. We conducted a field experiment in which 100 mm (900 L) of Brilliant Blue dyed water was infiltrated across a 3- by 3-m area in relatively homogenous and undisturbed sandy alluvial sediments. Reflection GPR data were collected before and after infiltration. Dye-staining patterns, revealed by excavating a 2-m-deep trench through the infiltration area, were compared with changes in the GPR data. Reflection amplitude changes as well as reflection delay revealed significant differences within the dye-stained area. The GPR data provided information about the unsaturated flow below the extent of the dye staining, and the results of the synthetic GPR modeling, as well as the observed changes in the real GPR data set, underline the potential of reflection GPR as a nondestructive method to map unsaturated flow phenomena.

Infiltration of dyed water into sandy alluvial sediments revealed that preferential flow occurs at several scales and dominates the flow field. The infiltration flow field was imaged by reflection ground penetrating radar (GPR), the results of which demonstrated the ability of reflection GPR to nondestructively image the distribution of moisture in preferential flow fields.

The quality and quantity of groundwater resources depend highly on the flow and transport properties of the unsaturated zone. Traditionally, unsaturated hydraulic parameters are estimated using retention and hydraulic conductivity experiments performed on small soil samples in the laboratory. These small-scale analyses are not able to describe preferential flow paths or unstable wetting fronts that are observed at the field scale. Prediction and characterization of flow in the unsaturated zone was found to be difficult and highly complex in several studies using dye tracers for visualization of flow paths (Flury et al., 1994; Schmalz et al., 2002; Weiler and Flühler, 2004). Large-scale experiments on both homogeneous and heterogeneous soil columns have shown that fingering and preferential flow exist and dominate the flow field for soil types ranging from clayey soils to unstructured sandy sediments (Flury et al., 1994; Wildenschild et al., 1994, Schmalz et al., 2002). Šimůnek et al. (2003) and Kung (1990) furthermore showed that the infiltration of water seemed to be dependent on the dip of geologic layers and that the dipping layers could also induce fingered flow.

Dye tracing infiltration experiments are powerful techniques for visualizing the dynamics of unsaturated water flow and have been widely used for many purposes (Flury and Wai, 2003). The popularity of Brilliant Blue as a dye tracer is attributed to its excellent visibility in most soil types, its nontoxicity, and its solubility in water (Flury and Wai, 2003). Brilliant Blue has been found to be subject to retardation because of nonlinear absorption, and, in consequence, it cannot be considered an ideal, conservative tracer (Kasteel et al., 2002). Thus Brilliant Blue is not suitable for the estimation of water flow properties but is ideal for illustrating water flow pathways in soils.

A major drawback in using dye tracing for unsaturated flow characterization is that the method is destructive and unrepeatable at the same location because the survey area needs to be excavated and the method therefore only provides a snapshot of the flow. Alternatively, borehole sampling and other coarsely sampled point measurements do not provide the same detailed information about small-scale unsaturated flow patterns. Dense point sampling is effectively as destructive as excavation of a profile and furthermore such sampling may potentially generate artificial flow paths.

Reflection GPR is a nondestructive geophysical method that is sensitive to changes in moisture content in the unsaturated zone (Topp et al., 1980; Davis and Annan, 1989; Reynolds, 1997) and has therefore been used for soil moisture estimation in several studies (Trinks et al., 2001; Huisman et al., 2003; Truss et al., 2007). Trinks et al. (2001) used time-lapse GPR to monitor the movement of an induced plume of water in dry sand in the laboratory. By collecting densely sampled GPR lines every 30 min and comparing the time-lapse data sets to the baseline signal, they were able to follow the infiltration plume in three dimensions with time and distinguish preferential flow. Truss et al. (2007) conducted a field experiment in which they used reflection GPR to monitor both natural and forced infiltration of water in oolithic limestone containing sand-filled dissolution sinks at the surface. For natural infiltration, bypass flow occurred that was highly dependent on the sediment type and geologic boundaries, while a forced point injection of water in the dissolution sinks allowed three-dimensional monitoring of the downward water migration (Truss et al., 2007).

The aim of this study was to assess the applicability of densely sampled reflection GPR surveys for quantitative monitoring of water movement in unsaturated sandy alluvial sediments after forced infiltration of water. The addition of Brilliant Blue dye to the infiltrating water and subsequent excavation of a trench allowed visual assessment of the movement of the wetting front, and these results were then compared with the high-resolution GPR data. To investigate the sensitivity of the recorded GPR signals to various infiltration patterns, forward modeling of synthetic subsurface scenarios was performed. The simulation results demonstrate how changes in moisture content affect the reflected GPR signal and provide insights into how it is possible to resolve small changes in moisture content in the subsurface using reflection GPR.

Materials and Methods

Ground Penetrating Radar

Ground penetrating radar wave velocity and signal attenuation depend on the relative dielectric permittivity, ε, magnetic permeability, μ, and electrical conductivity, σ, of the subsurface as well as the frequency of the signal (Davis and Annan, 1989). The magnetic permittivity is often disregarded in nonmagnetic materials such as sandy sediments.

The dielectric constant can be related to the electromagnetic wave velocity, v, in the following manner for low-loss geologic materials (Davis and Annan, 1989): 
where c is the speed of light. If the electromagnetic wave velocity or dielectric properties are known, the moisture content of a given subsurface can be calculated using the empirical relationship of Topp, which relates moisture content θ to the dielectric constant (Topp et al., 1980): 
In this study, we used the simpler, approximate form of Topp's equation developed by Ferré et al. (1996): 
which has been shown to provide results that differ <5% from Topp's original equation for volumetric water contents between 5 and 40% (Ferré et al., 1996).

Contrasts in dielectric properties of the subsurface cause parts of the emitted electromagnetic signal to reflect back to the receiver at the surface (Davis and Annan, 1989). Following Eq. [1] and [2], water in a porous material lowers the electromagnetic wave velocity by increasing the dielectric constant. Abrupt changes in moisture content produce high reflection coefficients because of large GPR wave velocity differences, and these show up in the resulting radargram as reflections.

The vertical resolution depends on frequency and bandwidth, and can be approximated by one quarter of the pulse length (Reynolds, 1997). The horizontal resolution depends on the antenna radiation pattern and signal frequency, and decreases with depth as the radius of the first Fresnel zone increases and the signal moves downward in a conical beam of increasing width (Reynolds, 1997). Migration can improve the horizontal resolution, however, by focusing the signal and restoring true reflection dips (Yilmaz, 1987; Grasmueck et al., 2005).

Field Site

For this study, we chose a heath area, Hjelm Heath, situated in the northern part of Jutland in Denmark (see Fig. 1 ) as the field site. The area was formed as a glacial outwash plain during the last phases of the Weichselian glaciation when large amounts of meltwater were released from the melting glaciers and deposited tens of meters of meltwater sediments, predominantly sand (Houmark-Nielsen et al., 2005). A geologic setting such as this is very common for the western part of Denmark, where most of the sandy soils are used for agriculture.

The vegetation cover at the site consists of heather [Calluna vulgaris (L.) Hull], wavy hair grass [Deschampsia flexuosa (L.) Trin.], and black berry (Empetrum nigrum L.). The root zone extends to a depth of 60 cm (Ladekarl et al., 2005), and the soil development is typical of a podzol (Dalsgaard, 1998). A hardpan layer, formed by the precipitation of organic matter, is located 20 to 30 cm below the surface. The top 6 m of soil below the podzol consists mainly of medium to coarse sand (68–99% w/w). The texture varies with depth, some layers being rich in silt and clay and some layers rich in gravel and coarse sand (Ladekarl et al., 2005). With sand of different grain sizes being the dominant texture of the subsurface at Hjelm Heath, the area is suitable for reflection GPR surveys (Davis and Annan, 1989; Reynolds, 1997). The area around Hjelm Heath has not been subject to significant deformation after deposition, and the sandy sediments are assumed to be fairly undisturbed. The groundwater table is found 21 m below the surface (Ladekarl et al., 2005).

A well was installed in 1992 equipped with 11 time domain reflectometry (TDR) probes to a depth of 7 m. Some of the probes were placed at the same depth but on opposite sides of the well (Ladekarl et al., 2005). The distance from the TDR well to the survey area is approximately 20 m. Continuous measurements of moisture content with time have shown that duplicate probes at the same depth do not exhibit similar moisture content development. Wetting fronts were observed to bypass certain probes, and observed responses to infiltrations were in some circumstances recorded at deeper probes before being observed closer to the surface (Ladekarl et al., 2005). These observations clearly suggested that for all soil moisture conditions, water movement arising from natural precipitation was highly nonuniform and heterogeneous.

Infiltration Experiment

The topsoil (upper 10–15 cm) and vegetation including roots were removed before data collection as well as an additional 30 cm of soil. The decision to do so was made after a test infiltration in which we applied 19 mm of water across an area of 12 m2 of vegetation-free soil and found that this amount of water could not infiltrate as quickly as we had hoped. Thus, the final infiltration took place approximately 45 cm below the original surface and directly into the more permeable alluvial sand beneath the podzol horizon. A Brilliant Blue (BB) solution with a concentration of 5 g BB/L was prepared. The conductivity of the solution was 0.287 S/m, which is considerably higher than fresh water (∼0.001 S/m) and what we expected to represent the natural moisture in the subsurface at Hjelm Heath.

The infiltration experiment was inspired by the works of Flury et al. (1994), Javaux et al. (2006), and Kung (1990), who performed infiltration experiments using BB as a tracer. Flury et al. (1994) infiltrated 40 mm of BB-dyed water (concentration 4 g/L) in 14 different soil types and found that after 8 h, the penetration depth was in many cases considerably <1 m but that for some soils extreme bypassing was seen. Javaux et al. (2006) infiltrated 180 mm of BB-dyed water (concentration: 1 g/L) into a 1-m sand column and within 8 h the dyed water had not reached the maximum depth.

We chose to infiltrate a total amount of 100 mm (900 L) of BB solution during 2 h across an area of 3 by 3 m. Watering cans equipped with simple sprinklers for more uniform spreading of the water were used for water application. After infiltration, the area was left uncovered, and a small rain shower delivered a couple of millimeters of precipitation during the night, which was negligible compared with the infiltrated 100 mm. Evapotranspiration was of no significance because the infiltration experiment was conducted in early April, at which time evapotranspiration is low. One day after infiltration, a trench was excavated through the infiltration area to a depth of 2 m for visual assessment of the dye tracing patterns outlined by the BB dye (see Fig. 1 for location and extent of trench excavation). Table 1 summarizes the timing of events in the infiltration experiment.

Ground Penetrating Radar Data Acquisition

Ground penetrating radar lines were collected using a PulseEKKO Pro System (Sensors and Software, Mississauga, ON, Canada) equipped with shielded 250-MHz antennas and an odometer wheel for trace spacing measurements. We used the following acquisition parameters: time window, 500 ns; stacking, 16; sampling interval, 0.4 ns. For precise data collection, we constructed a setup consisting of two hollow polyvinyl chloride slides along which the GPR instrument was pulled. Before carrying out the actual field work, it was verified that the slides did not adversely affect the GPR data collected by the shielded antennas. The slides were 10 by 10 by 600 cm and were manually fixed to the ground by bamboo sticks in order not to obscure the data collection. Careful repositioning of the GPR instrument for each line along the slides assured that the repeated data sets were comparable line to line.

Before and after infiltration, a total of 100 GPR lines were collected with a line spacing of 0.05 m (see Table 1 for timing of data collection). For all lines, the trace spacing was also 0.05 m, thus providing a uniformly sampled three-dimensional data set. Each line consisted of 94 traces, making the final data acquisition area 4.95 by 4.65 m (23.02 m2). Because BB dye was infiltrated over a 3- by 3-m area within the GPR area (see Fig. 1), the GPR data set collected after infiltration contained lines both affected and not affected by direct infiltration. Also, each GPR line crossing the infiltration area had several traces at the beginning and end that were not affected by BB infiltration. For details of the field site layout, see Fig. 1.

The GPR signal was delayed and attenuated by the extra amount of water (Reynolds, 1997; Truss et al., 2007) and the higher conductivity of the infiltrating BB solution further attenuated the signal. In areas of relatively high electrical conductivity, the conductivity of the studied formation may significantly affect the velocity and the validity of the relations based on the low-loss criteria such as Eq. [1–3] (Giroux and Chouteau, 2010). For the amount of conductive BB solution that we applied in this experiment, however, we assumed that these effects related to the bulk conductivity of the studied formation were negligible and thus we assumed that any decrease in velocity was merely due to an increase in moisture content. There is a chance that some reflections were more pronounced due to ponding of water at geologic boundaries, which would cause the reflectivity coefficient to increase.

Ground Penetrating Radar Data Processing

The collected GPR data were processed using the EKKO View Deluxe and MATLAB (MathWorks, Natick, MA) software packages. For all lines, a dewow correction was applied as the first step in the processing to remove low-frequency noise. Pre- and post-infiltration GPR lines were then plotted in MATLAB and the same simple scaling function was applied to all lines.

The individual GPR data lines were migrated (two-dimensional migration) using the EKKO View Deluxe software to remove diffraction hyperbolas and focus the signal. The same migration velocity was applied to the pre- and post-infiltration data sets. For detailed analysis of single reflections in the GPR data, MATLAB was used to pick reflections in all lines for both data sets, thus obtaining the arrival time and amplitude for each reflection in the entire data set. The arrival time differences could then be translated into electromagnetic wave velocity changes and converted into approximate moisture content changes using Eq. [3] and information about the thickness of the geologic layers obtained from the excavated profiles. Analyses of multiple reflections provided a basis for estimating the spatial distribution of the water in the subsurface and could potentially be validated against the BB dye-staining patterns.

Synthetic Studies

Configuration of Synthetic Models

To evaluate the possibilities and limitations of using reflection GPR to monitor unsaturated flow phenomena, a synthetic analysis was performed using the MATLAB-based two-dimensional finite-difference time-domain (FDTD) code developed by Irving and Knight (2005). Synthetic subsurface models are defined by specifying distributions of electrical conductivity, as well as dielectric permittivity and magnetic permeability. The two-dimensional FDTD code uses a source pulse, which is the normalized first derivative of a Blackman–Harris window function (Irving and Knight, 2005). Any type of synthetic model can be created and modeled using the code; however, the computation time increases with the complexity of the chosen subsurface characteristics and model discretization.

We created a synthetic geologic model of 10-m depth and 20-m width for Hjelm Heath based on TDR measurements as well as soil samples extracted during installation of the TDR well at the field site (Ladekarl et al., 2005; Rasmussen, personal communication, 2008). Dielectric permittivity values used in the model were calculated based on soil moisture measurements from the nearby TDR well. The subsurface part of the model consists of four layers that have different values of ε ranging from 4 to 6 corresponding to moisture content values between approximately 0.05 and 0.12, respectively. An air–earth interface was included by adding a layer with relative dielectric permittivity ε = 1 and electrical conductivity σ = 0 S/m to the top of the model. Because, for simplicity, we were mainly interested in the delay of reflections, the conductivity was set to σ = 0.001 S/m and the relative magnetic permeability was represented by its free-space value, μ = 1, throughout the model. During a real forced infiltration experiment, changes in electrical conductivity are expected to occur, caused by both the added amount of water and the potentially different electrical conductivity of the infiltrated water (compared with the incipient soil water). Thus, the amplitude variations and attenuation of the modeled signal in the synthetic tests were not affected as much as during a true infiltration experiment.

We chose simulation parameters closely resembling the real GPR equipment used in the field experiment. The signal had a center frequency of 250 MHz, the spacing between the transmitter and receiver was set to 0.4 m, and the simulations were made with a constant offset. The model discretization was 0.02 m.

To recreate an infiltration event, we included a wetting front from 0 to 1.5 m with increased moisture content (0.025 higher than the underlying layer) in the synthetic model. The middle of the infiltration was shaped like a 2-m-wide and 1-m-deep channel to simulate an unstable wetting front causing fingered flow (Fig. 2A ). Because lateral flow was also likely to occur, another model containing a wedge-shaped area of increased moisture content was also created (Fig. 3A ). Synthetic GPR radargrams were calculated only for a 2-m-wide subsection of the model, which is marked in Fig. 2A and 3A by the vertical dotted lines.

We assumed the spatial variability of moisture content and thus the dielectric permittivity to vary within each layer. Spatial variability was added to the input dielectric permittivity model as normally distributed and spatially uncorrelated random noise with standard deviations of 0.25 and 0.50, respectively. This amount of random noise corresponds to heterogeneity in the moisture content of up to 1.5%. Figures 2B to 2D and 3B to 3D show the variation in dielectric constant with depth for the three input models. We see how the change in dielectric constant, and thus in moisture content, becomes more diffuse when more heterogeneity is added.

It should be noted that all the synthetic analyses in this study were made using a two-dimensional code, the underlying assumption being that the constructed subsurface models continue indefinitely perpendicular to the presented image. Thus, “fingers” are more like trenches than cylinders. Wave scattering is expected to increase in true three-dimensional environments, and consequently the findings presented here are therefore perhaps overly optimistic, yet they give an indication of what changes we can expect to observe.


Figures 2E to 2G show the resulting synthetic radargrams that contain increasing degrees of heterogeneity within the subsurface for models with a channel-shaped increase in moisture content. If the subsurface is assumed not to be subject to internal variation in moisture content within each layer, the resulting synthetic radargram for the infiltration finger will show a distinct bowtie structure arising from its channel-like shape (cf. Yilmaz, 1987). This effect can be seen in Fig. 2E, marked with an arrow. Likewise, the wedge causes diffraction hyperbolas along its bottom (Fig. 3E, marked by an arrow).

In Fig. 2F, we see that for the lowest degree of added heterogeneity the bowtie structure is still visible in the resulting radargram (marked with an arrow). With the addition of more heterogeneity, the bowtie became more diffuse (Fig. 2G). The larger heterogeneity also made the deeper reflections harder to recognize. Figures 3F and 3G show the resulting synthetic radargrams that contain increasing degrees of heterogeneity within the subsurface for models with a wedge-shaped increase in moisture content. Even with a low amount of added heterogeneity, the extent and shape of the wedge becomes hard to distinguish, and as the heterogeneity increases, the wedge disappears almost completely.

For the infiltration finger, we observed delays of up to 2 ns and a changed dip of deeper reflections due to the increase in moisture content in the top of the model (marked with a dotted arrow in Fig. 2 for the lowest reflection). Even with strong heterogeneity that blurred the shape of the water front, these effects did not disappear and can thus provide information about the infiltration, even in cases where the direct effect cannot be seen. It can also be seen that the parts of the deeper reflections that were affected by the increase in moisture content in the infiltration finger had a greater lateral extent than the actual finger. This was due to the radiation pattern and Fresnel resolution of the signal and means that changes in the radargrams at depth are likely to represent actual changes on a smaller scale at shallower depths.

The increase in moisture content of 0.025 used for this analysis was assumed to be low, and for larger increases the infiltration finger would be more easily recognized, even in the case with large heterogeneity. Based on the results of the numerical modeling, we expect to be able to delineate the bulk changes in moisture content arising from a forced infiltration experiment using reflection GPR. In the case of a moisture content increase in large coherent volumes of the subsurface, it should be possible to see attenuation of the signal in the wet areas (not shown). The type of information that can be achieved through the GPR data is, however, dependent on the amount of heterogeneity in the subsurface.

Sedimentary Structures

The main excavation was a trench dug to a maximum depth of 2.2 m through the middle of the BB-dyed area (see Fig. 1). Only one side of the excavation was kept as a vertical profile, which can be seen in its full length and depth in Fig. 4 . The top of the profile was approximately 45 cm below the soil surface due to the removal of the topsoil and hard pan at the site. Figure 5 shows details from both vertical profiles as well as horizontal excavations.

The vertical profile reveals that the uppermost subsurface at Hjelm Heath consists of several different sedimentary packages of sand with varying grain sizes and structures. This concurs with grain size analysis results from the field site reported by Ladekarl et al. (2005).

The top 0.20 m of the excavated profile consists of sand that is highly disturbed by bioturbation and root remnants (Fig. 4, marked A), but it is still possible to distinguish cross- bedding. Ladekarl et al. (2005) reported that the active root zone at Hjelm Heath extends to 60 cm below the surface, i.e., only 15 cm into the excavated profile; however, the vertical profile reveals root structures to a depth of 0.5 m (i.e., 0.95 m below surface). These deeper root remnants are characterized by having black, hardened centers and are up to 5 cm in diameter. Because the active root zone only extends to 0.60 m below the surface and the deeper root remnants have been severely altered, it is believed that they are ancient tree roots >1000 yr old) from before deforestation of Western Jutland began (Andersen, 1994).

From 0.2 to 2 m in the excavated profile, the sediments are planar with a slight dip of approximately 2° toward the southeast. At the 1.4- to 1.5-m depth, we found a horizon with a dark orange color, probably stemming from precipitated Fe oxides (Fig. 4, marked B). This process is known as redoximorphism (Jacobs et al., 2002), which occurs when the groundwater table is both high and fluctuating. As such, this horizon is expected to have originated from the time around the last glaciation, when permafrost forced the groundwater upward during the winter (Jacobs et al., 2002). Due to redoximorphism, Fe precipitates, seen as orange nodules, are much more prevalent below the horizon than above.

Around the 2-m depth, there is a transition from planar to cross-bedded sediments as well as a decrease in grain size. The interface between the two sedimentary packages dips toward the northwest; this is best seen in Fig. 5A, marked 1.

Dye-Staining Patterns

Dye-staining patterns arising from the infiltration experiment can be seen in Fig. 4 and 5. In Fig. 4, the blue arrow indicates the extent of the infiltration at the surface. Note that the edges of the profile have not experienced infiltration.

We observed that the top 1 m of the profile below the infiltrated area was almost entirely blue due to dye staining, although areas existed within the top 1 m where bypass flow had occurred, as in Fig. 4, marked C. Below the 1-m depth, the dye-staining pattern became more irregular and uniform matrix flow appeared to no longer be the dominant flow process. In the middle of the profile, we observed three large flow fingers (Fig. 4, marked D) that had facilitated the deepest dye infiltration. The fingers had a maximum width of 0.2 m in the excavated plane; however, their three-dimensional structure was not revealed by the main excavation. Sedimentary structures within the large fingers did not reveal any variation that would make them more prone to water flow; however, they appeared to initiate at the same depth as the redoximorphic horizon. The precipitates probably prevented flow across sections of the horizon, causing the infiltrating water to break through a few permeable windows, forming infiltration fingers.

Toward the edges of the infiltration area, water had infiltrated laterally to the sides (Fig. 4, marked E). Detailed inspection of the dye-staining patterns revealed that BB-stained water had moved along small sedimentary or structural boundaries in the layered sediments. The irregularity of water flow was in part initiated by grain size differences in the sand acting as small capillary barriers, as illustrated in Fig. 5B. Small-scale flow along horizontal or dipping boundaries between changing lithologies, as well as diffusive flow, can thus be responsible for transport of water out of the infiltration area.

With the purpose of investigating the three-dimensional nature of the infiltration patterns, we removed 0.3 m of sediment from the central 1 m of the main profile (for exact location, see Fig. 1 and 4). The 1-m-wide and 2-m-high profile can be seen in Fig. 5A. Brilliant Blue dye staining in the top 1 m of this profile was less coherent compared with Fig. 4, with almost a third of the area completely undyed. The large unstained or bypassed area in the top of the profile was caused by the presence of an ancient root (Fig. 5A, marked 2), which diverted the water flow. In the horizontal view, some of the ancient roots appeared as circular black structures with a periphery of undyed sand around them, suggesting hardening of the surrounding sand due to precipitation of organic matter (see Fig. 5C). Although the black root remnant in Fig. 5A extended 0.25 m into the profile, the surrounding sand continued to be bypassed to a depth of 0.65 m below the profile top.

Isolated patches of dyed sand were found in numerous places in the profile (see Fig. 4, marked F) and further illustrate the heterogeneous nature of the infiltration. Notice also that the three

large infiltration fingers seen in the main profile (Fig. 4, marked D) are not visible at this position. Rather, there is one large wide infiltration path to the right, which seems to extend below the depth of the excavation.

The “floor” of the excavation pit at the 2-m depth revealed large blue stains (not shown here), indicating that there had been percolation of infiltrated water deeper than 2 m into the subsurface, which was not visualized by the dye-staining patterns in the main profile.

The overall infiltration patterns seen in the excavation at Hjelm Heath are similar to the flow type “heterogeneous matrix flow and fingering” as defined by Weiler and Flühler (2004), with most of the fingers being rather large in diameter with a high connectivity between the dyed areas. The bypass of large areas at the top of the profile caused unstable flow, disrupting homogeneous infiltration. This was further accentuated by the slightly hardened redoximorphic horizon that created infiltration fingers.

Ground Penetrating Radar Data

Figure 6A and 6C show the raw and migrated preinfiltration GPR Line 44, respectively. A two-dimensional migration velocity of 0.125 m/ns was chosen based on the average soil moisture measurements from the nearby TDR. The 100 preinfiltration GPR lines show the same overall distribution of reflections in the subsurface, which can be divided into three sections. The top section consists of parallel reflections dipping slightly toward the southeast, which corresponds to the planar bedding found in the top 2 m of the excavated profile. Below the parallel reflections, we find a more heterogeneous collection of less continuous reflections, with several distinct dipping reflections across the entire profile (Section 2). In the unprocessed data, this feature appears as a bowtie (c.f. Yilmaz 1987), yet in the migrated lines it resembles several generations of a buried river channel. Due to the short length of each line, however, it is generally difficult to determine what this relatively deep feature really is.

Several reflections can be recognized with a high degree of connectivity throughout the entire data set. From Section 3 and below, the reflections become less pronounced and the radargrams are mostly dominated by multiples and noise. The last horizontal reflection in the bottom of Section 1 appears to fit the interface between the two textures of sand seen in the excavation at approximately the 2-m depth (Fig. 4 and 5A).

Figures 6B and 6D show the raw and migrated post-infiltration GPR Line 44, respectively, and also here we see that it is possible to distinguish the three sections. Five reflections have been marked in both the pre- and post-infiltration radargrams (Fig. 6C−6D). Furthermore, the outline of the dyed area in the main excavation has been superimposed onto the post-infiltration radargram.

Ground Penetrating Radar Signal Changes

An important observation in comparing unmigrated radargrams from before and after infiltration is that no new diffraction features are apparent after infiltration (Fig. 6B). We would expect large contrasts in moisture content to create diffractions and disturb the image of the subsurface following the synthetic analyses described above. Since we can clearly recognize reflections at depth after infiltration, this is not the case.

A closer look at the single traces of each radargram reveals that amplitudes for traces affected by infiltration are diminished after infiltration while traces from outside the infiltration area do not show the same tendencies. On several traces, especially the ones inside the infiltration area, the two data sets differ somewhat at the very top, i.e., within the first couple of waveforms. This could be attributed to minor changes in the surface characteristics. To collect the GPR data sets during infiltration, it was necessary to walk within the experimental area, which probably affected the uppermost sediments. Furthermore, the surface was raked during infiltration, which surely disturbed the top few centimeters of sand. These effects, as well as the minor uncertainty accompanied by repositioning the GPR equipment along the slides can probably explain the difference in signals at shallow depths.

Inspection of GPR lines from within the infiltration area reveals that some reflections toward the middle changed and are less coherent and pronounced after infiltration. An example of this is pointed out in Fig. 6D (marked 1). Also, the signal around Reflection 4 in the middle of the radargram appears to have become attenuated after infiltration (Fig. 6D, marked 2). Both of these effects were probably due to the presence of infiltrated water.

Below Reflection 3 at approximately the 2.5-m depth, an otherwise incoherent reflection became more pronounced after infiltration. It is best seen in the unmigrated data, Fig. 6B, marked 3. This could be attributed to ponding of water on one side of a geologic boundary, thus causing the reflection coefficient to increase due to the higher contrast in GPR velocity between the two media. Similar observations were not as pronounced in lines outside the infiltration area, which further indicates that the observed changes were caused by the infiltrated water.

In Fig. 6E, the five marked reflections from before and after infiltration are shown together. We see clearly how they were delayed, especially within the infiltration area, which is marked by dotted lines. As mentioned above, such a delay must be attributed to a change in moisture content between the pre- and post-infiltration data acquisitions. If the same comparison is made for reflections from lines outside the infiltration area, no delay is found, suggesting that the moisture content was unchanged.

The delays are most prominent for Reflections 3, 4, and 5, with Reflection 4 being the most delayed. This can partially be explained by the delay accumulating through the radargram but also indicates that the bulk amount of water had percolated through the top section and was now found at greater depths.

Reflection Delay

The arrival time and amplitude of the five reflections shown in Fig. 6C to 6E were selected out of the entire data set. The difference in arrival time for each reflection was subsequently calculated, and the result is illustrated in Fig. 7 . For Reflection 1, Fig. 7A shows the delay for the entire volume between the surface and the reflection, whereas Fig. 7B to 7E represent the delay between each selected reflection and the previous one. Blue colors indicate delays. Due to the short length of each GPR line, the channel structure that constitutes Reflection 5 is not always well determined and therefore the results in Fig. 7E should be interpreted with caution. In spite of the inherent picking difficulties and apprehensions toward the reflection delay analysis, Fig. 7 shows continuous areas experiencing similar delays for the picked reflections that cannot be dismissed as errors.

The observed delay was mostly confined to the infiltration area, indicated by a black box in Fig. 7A to 7E. The delay of Reflection 1 is mainly found in the southern part of the area with no delay found in the northern parts; Reflections 2 and 3 were delayed across the entire infiltration area, and for Reflection 4 the delay was mostly concentrated in the northeastern part of the infiltration area. The large red areas in Fig. 7E illustrate how Reflection 5 was not further delayed compared with Reflection 4.

The delay also extended outside the infiltration area for the other reflections and could be caused by lateral flow, as was also observed in the dye-staining images (Fig. 4 and 5). It should be noted here that because propagation of the electromagnetic wave is three dimensional, the received signal represents a volume of the subsurface that becomes larger with depth (cf. Reynolds, 1997). A change in moisture content or conductivity within the infiltration area can thus give rise to changes in GPR signals emitted from the surface outside the infiltration area. Migration of the GPR data does not remove this effect entirely, although it will be less for shallow reflections.

The delay distribution between the different reflections and associated moisture content change can be correlated with the observations in the radargrams in terms of dipping reflections, as well as the sedimentary structures revealed in the excavation. Hence, the water movement at Hjelm Heath on this scale seems to be controlled by the dip of sedimentary structures in the subsurface.

In Fig. 7A to 7E, the position of the excavated trench is marked by a dotted line through the experimental area. We found that, based on reflection delays, there was almost no extra water present in the first ∼1.5 m of the profile (i.e., down to Reflection 1) after infiltration. The profile revealed extensive dye staining of the sediments; however, the reflection delays indicated that the infiltrated water must have percolated beyond this depth.

Moisture Content Measurements

Ten soil samples were collected from the main profile for volumetric moisture content measurement (for exact locations, see Fig. 1). The samples were collected every 0.20 m, the first at a depth of 0.20 m from the top of the profile (i.e., 0.65 m below the surface). The results of the moisture content measurements are shown in Fig. 8 . The moisture content increased from θ = 0.05 at the top to θ = 0.12 at a depth of 2.45 m below the soil surface. The samples that were collected from the excavated profile ranged from being completely dyed to undyed; however, no connection could be made between the measured moisture content and the degree of dye staining.

Moisture content profiles obtained from TDR measurements in the nearby well are also depicted in Fig. 8 (Rasmussen, personal communication, 2008). We have included average as well as minimum and maximum measurements for the year 2008 as well as the average moisture content for the week in 2008 during which the field work took place. Note that the depth axis represents meters below the vegetated surface.

The TDR data show an increase in moisture content from 0.06 to 0.11 within the first couple of meters; thereafter it decreased to 0.05 between 2 and 2.5 m. The measurements after infiltration were in general slightly higher than the mean moisture content measured by the TDR probes during the field campaign. This difference could be caused by the forced infiltration experiment but could also arise from changes in soil properties between the two measurement points because there was approximately 20 m from the TDR well to the experimental area.

The largest difference in moisture content was found at the 2.45-m depth, where the samples showed moisture contents around 0.12 while the TDR data showed values of 0.05. Such an increase in moisture content should, according to the synthetic studies, affect the GPR radargrams substantially, which, however, was not the case. The TDR data have here been assumed to represent background values of moisture content for this area; however, it is possible that the background moisture content at the TDR profile was different from that of the field site due to changes in sedimentology. If, instead, an average moisture content increase of 0.02 due to the general tendencies seen in Fig. 8 is assumed, this change would not give rise to major changes in terms of diffraction effects in the GPR data.

Amplitude Analysis

An analysis of amplitude changes between the two GPR data sets can provide insight into the overall impact of the infiltration experiment on the GPR wave energy distribution. An increase in moisture content or conductivity near the top will attenuate the signal and cause increased changes with depth. As was the case with the reflection analysis, changes in the amplitude root mean square between the two data sets can be interpreted to represent changes in moisture content or conductivity.

Amplitude analysis was performed for the migrated data set in three intervals of the subsurface and the results are shown in Fig. 9 . Each interval has been marked with different colors in Fig. 9D. The first section encompasses the entire signal above Reflection 1 (Fig. 9A, green in Fig. 9D), whereas the second and third intervals are considered between Reflections 1 and 2 (Fig. 9B, blue in Fig. 9D) and Reflections 2 and 4 (Fig. 9C, red in Fig. 9D), respectively. The first two intervals incorporate all the horizontal reflections in the top of the radargrams, whereas the third amplitude analysis interval contains the upper part of Section 2 (Fig. 6C).

Root mean square values were calculated from the amplitude values of each GPR line in the given time interval, and the difference between the before and after data set was calculated by simple subtraction. We used the arrival time of the specific reflections obtained previously to determine the upper and lower boundaries of the intervals of interest, thus automatically taking into account the time shifts occurring in the radargrams. This method renders results that merely show the change in GPR signal amplitude within the part of the subsurface between the two reflections delineating the analysis interval.

The outline of the infiltration area is shown in Fig. 9A to 9C as well as the position of the excavated profile (dotted line). Blue colors indicate that the signal has been attenuated, whereas reddish colors show areas where the amplitudes have increased after infiltration. For the analyzed intervals of the subsurface, it is evident that the GPR signal was subject to change primarily within the infiltration area but also up to almost 1 m outside the infiltration area. Note that the area influenced by changing amplitudes outside the infiltration area increased between the three sections. As mentioned above, part of this tendency can be attributed to the three-dimensional propagation of the electromagnetic waves, which caused lines outside the infiltrated area to be affected by the changes caused by the infiltration.

For the first interval, Fig. 9A, we found that the change in amplitudes was not particularly significant but that the decrease in amplitude was mostly confined to the infiltration area. This observation is consistent with the findings from the reflection analysis (Fig. 7), which showed that Reflection 1 experienced less delay than the other reflections analyzed. The small change in amplitude within the entire infiltration area of the first interval can be attributed to the increase in conductivity of the soil water caused by the BB dye present in the upper layers (as seen in the excavations, Fig. 4 and 5). It is apparent for the second and third intervals (Fig. 9B and 9C) that the main change in amplitude occurred within the infiltration area. It is also worth noticing that for all intervals the primary change within the infiltration area was toward lower amplitudes, i.e., an attenuation of the signal. An increase in amplitude is seen primarily along the edges of the infiltration area, which can be explained by lateral flow at the edges of the wetting front and subsequent increases in reflectivity at geologic boundaries just outside the infiltration area.


Many previous studies have attempted to describe preferential flow paths and unstable wetting fronts. These investigations have rarely included moisture content measurements and subsequent assessment of the whereabouts of the infiltrated water. In Flury et al. (1994) and other studies, the soils were classified as either initially “dry” or “wet” but the actual moisture contents of these soils were not measured either before or after infiltration (Flury et al., 1994; Weiler and Flühler, 2004). Hence, the results of these studies merely illustrate the movement of the infiltrated dyed water and as such do not shed light on the overall movement of soil water in the surrounding subsurface.

Our results indicate that the excess water from the infiltration was not present in the dyed areas and that displacement flow was responsible for distributing the water farther into the subsurface. Moisture content increases were seen deeper than 4 m into the subsurface although the dye staining in the excavation was most prominent in the top 2 m of the soil. Whether BB dye in selected places had moved as far into the subsurface as the water is unknown because we only had visual assessment of the water and dye movement from the excavation. The observation of large stained areas in the bottom of the excavation, however, could be an indication of deeper infiltration of the dye.

According to the texture analyses performed by Ladekarl et al. (2005), the sand between 3 and 5 m is much coarser than the surrounding geology and contains much gravel compared with the other samples analyzed. The part of the subsurface between 3 and 5 m corresponds roughly to the radar signals between Reflections 4 and 5. Assuming that the geology also represents the experimental site used for the present study, this has two implications: (i) the coarser layer facilitated a much faster infiltration to deeper layers, which cannot be distinguished in the GPR data set beyond Reflection 5 because of noise, multiples, and loss of energy; or (ii) the coarser layer was acting as a large capillary barrier (because the general moisture content was still low) and thus water was moving laterally on top of this layer around Reflection 4. The latter is in agreement with the assumption that the water is likely to move along the geologic boundaries, as the reflection analysis also highlighted.

Compaction has been shown to induce preferential flow in sandy loams (Mooney and Nipattasuk, 2003), and as it was necessary to walk inside the infiltration area while the GPR data sets were collected and during infiltration, this could have caused a slight compaction of the upper sediments and thus influenced the infiltration. Cracks and small deformation could have caused the water to bypass some areas while flow could be facilitated in others, but on excavating the profile we did not find evidence that such compaction factors governed the overall water flow. Instead, we observed how ancient roots induced the bypass flow and that unstable flow and fingering was induced at a cemented horizon.

Assuming an average moisture content between 0 and 2 m of 0.075 before infiltration (value based on TDR measurements in Fig. 8), the amount of water in this volume would be 1350 L. The addition of 900 L of water through infiltration would increase the moisture content to 0.125, which was only measured in the deepest sample location (2.65-m depth). The reflection delays do not correspond to moisture content increases of this magnitude. If the average delay and depth of Reflection 2 based on the results from Fig. 7 are assumed to be 0.6 ns and 2 m, respectively, the increase in moisture content corresponds to an additional 200 L of water. This means that the remaining 700 L of infiltrated water has either moved laterally out of the area or percolated deeper. Because there was little evidence from either the excavation or the GPR reflection analysis that such large amounts of water had moved to the sides, it must be concluded that the infiltrated water had passed through to the subsurface beneath Reflection 2. This is also evident from the increasing delays of Reflections 3 and 4, which suggests that most of the infiltrated water was located between these reflections. Figure 7 also shows that the delay of Reflection 5 was less than that of Reflection 4, which suggests that there was no extra water present in the subsurface between these two reflections.

By performing a two-dimensional migration, we were correcting the influences introduced by the radiation patterns of the GPR signal; however, these radiation patterns are not restricted to the line along which the migration was performed. The signal is also influenced by off-line effects; however, the processed data shown here clearly indicated that reflections were delayed primarily within the infiltration area. True three-dimensional migration would potentially further improve the analysis (Grasmueck et al., 2005).

It has been documented that BB dye is subject to retardation (Ketelsen and Meyer-Windel, 1999; Kasteel et al., 2002) and thus increased moisture content is expected farther into the subsurface than what was outlined by the tracer. Simulations of water and tracer movement in the subsurface using the HYDRUS-1D model (Šimůnek et al., 2008) and specifying an appropriate retardation coefficient showed that the outlined increase in moisture content to approximately the 4-m depth as well as a tracer migration depth of about 2 m is achievable for a sandy soil and for a high value of saturated hydraulic conductivity (500 cm/d).

Dye staining was clearly visible in the excavations for the top 2 m of sediment (Fig. 4 and 5), but the infiltrated water did not have significant influence on the GPR data for this part of the subsurface, neither as attenuation in the dye-stained area nor as scattering effects, which would be expected based on the GPR forward modeling (Fig. 2 and 3). The reason for this discrepancy between expected and actual GPR signals is that most of the infiltrated water was no longer present in the dye-stained areas, as the moisture content samples also showed. Apparently, the infiltrating water did not create sharp transitions between areas of high and low moisture content, which further diminished the influence on reflection GPR data. We did find several interesting differences, however, between the pre- and post-infiltration data sets. Reflection delay effectively illustrated the changes in moisture content at depths to approximately 5 m, and amplitude changes were also apparent at this depth. These observations were possible due to the presence of the distinct and easily recognized reflections in the subsurface. This was also the case for Truss et al. (2007), who relied on reflections from interfaces between different sedimentary packages for their GPR data analysis.

The infiltration experiment performed in this study did not represent naturally occurring conditions in Denmark because it is highly unlikely that a rain shower would produce 100 mm of precipitation within 2 h (annual precipitation is 781 mm [Danish Meteorological Institute, 2009]). The purpose of this experiment, however, was not to simulate natural conditions but rather to assess whether useful information regarding unsaturated flow can be obtained using reflection GPR surveys. In future investigations of the type presented here, it could be of interest to apply lower infiltration rates and conduct the infiltration for a longer period without removing the topsoil to better mimic natural conditions. The number of GPR surveys should be increased accordingly to follow the temporal development of the moisture profiles. Also, it would be interesting to see the effect of the BB dye staining on the data by conducting a similar experiment without the tracer, thus removing the effect of the increased conductivity of the applied water, as well as the unknown effect of the viscosity and density of the BB solution.

The experiment performed at Hjelm Heath is in some ways similar to the work presented by Truss et al. (2007), who used time-lapse GPR to monitor natural infiltration; however, they did not assess the infiltration patterns on smaller scale using tracers. The preferential flow observed in the excavation at Hjelm Heath was mostly controlled by macropores and soil heterogeneities. The GPR reflection amplitude analysis and the delay assessment suggest that funneling due to sloping layers was mostly responsible for the generation of preferential flow. The dye tracer patterns found at Hjelm Heath are similar to the findings of Kung (1990), whose field experiment was performed in a comparable environment. Because Hjelm Heath represents soil and subsurface conditions found many places in western Denmark, the obtained results have a broad applicability.


The BB staining revealed that at this field site the unsaturated flow regime is controlled by unstable wetting fronts as well as bypass flow, fingering, and lateral flow. High-resolution reflection GPR was found to be a suitable nondestructive method for mapping complex flow phenomena in the unsaturated zone. The quasi-three-dimensional GPR data set provided an excellent means of assessing both reflection delays caused by electromagnetic velocity changes as well as amplitude differences caused by large-scale movement of water. These results provided more knowledge of the flow regime at the field site. The observed effects were mainly confined to the infiltration area, with GPR data from outside being less or not affected at all.

The changes in GPR signal turned out not to be as significant as expected based on the numerical modeling as well as other experiments described in the literature. We did not manage to resolve the small infiltration patterns outlined by the dye tracer using reflection GPR. Brilliant Blue dye staining is helpful for illustrating and highlighting the existing flow paths; however, due to retardation of the dye, it does not shed light on the deep unsaturated flow. With high-resolution GPR it was possible to map the deeper infiltration, which in this case was substantial, with water percolating much deeper than revealed by the dye staining. The method can also be used in time-lapse measurements of soil moisture development, which will allow for even more detailed assessments of the nature of unsaturated flow.

This work was supported by Geocenter Denmark and HOBE–Center for Hydrology (www.hobe.dk), which is funded by the Villum Foundation. Lars Rasmussen provided invaluable assistance during the field work, and Martin Laier, Aarhus University, is thanked for his hard work removing the topsoil.

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