- Soil Science Society of America
Borehole ground penetrating radar (BGPR) operated in zero offset profiling (ZOP) mode has promise for monitoring rapidly changing water contents within the subsurface. However, the coexistence of multiple travel paths through the subsurface can give rise to measurement errors. Specifically, in layered systems with sharp changes in water content with depth, critically refracted waves may arrive before direct waves at some depths. Velocity profiles are determined based on analyses of the travel time of the first-arriving energy at each depth. Therefore, correct velocity analysis requires that these travel times be classified according to the path followed by the first-arriving energy. We establish criteria that can be used to identify first-arriving critically refracted waves from travel time profiles. Through hypothetical examples and a field experiment, we demonstrate that these criteria allow for more accurate determination of the water content profile. However, these corrections are limited if thin, high water content layers are present in the subsurface.
- ALM, alternating layered model
- BGPR, borehole ground penetrating radar
- EM, electromagnetic
- MOG, multioffset gathering
- ZOP, zero offset profiling
Within the past 10 years, near surface geophysical exploration has become invaluable to hydrologists as an inexpensive and rapid means of intermediate-scale subsurface monitoring, with support scales ranging from 0.5 to 30 m3 (van Overmeeren et al., 1997; Weiler et al., 1998; Huisman et al., 2001). Hydrologic site characterization with traditional tools, such as pumping–injection tests, are expensive and time-consuming. Additionally, these tools typically offer low spatial resolution. Other traditional hydrologic methods, such as coring, allow very high spatial resolution with very small sample volumes (several cubic centimeters). However, many cores are needed to characterize a site, and upscaling of small-scale measurements can be problematic (Wen and Gomez-Hernandez, 1996; Sanchez-Vila et al., 1996). Borehole ground penetrating radar, an electromagnetic geophysical instrument, has become popular among hydrologists because of its ease of use and simple interpretation. Recent documented uses of BGPR in hydrology include water content profiling (Binley et al. 2001; Alumbaugh et al., 2002), fracture delineation (Tsoflias et al., 2001; Sato and Takeshita, 2000), water table monitoring, facies delineation (Turner et al., 2000; Wanstedt et al., 2000; Bellefleur and Chouteau, 2001), and determination of hydraulic conductivity (Hubbard et al., 2001; Binley et al. 2002).
Borehole ground penetrating radar can be used in many modes, such as multioffset gathering (MOG) and ZOP. Multioffset gathering offers multidimensional imaging through high-resolution tomography, but is relatively slow due to the large number of measurements composing a survey. For example, Alumbaugh et al. (2002) reported a 6-h acquisition time for four profiles to 7 m with a 0.25-m vertical sampling interval. Hydrologic processes that are at steady state or very slow may be characterized with MOG. In contrast, the water content profile can be measured to the 15-m depth with a 0.25-m vertical sampling interval in 3 min in ZOP mode (Rucker and Ferré, 2002). As a result, ZOP is ideal for monitoring transient processes such as the change in water content beneath an ephemeral stream during flow or beneath a field receiving applied irrigation. When using BGPR in ZOP depth profiling mode, two antennae are lowered through parallel vertical access tubes of known separation. The centers of the antennae are located at a common depth for each measurement. An electromagnetic (EM) wave is then propagated from the transmitting antenna and measured with the receiving antenna. The procedure is repeated at a series of depths to produce a velocity profile of the subsurface. The velocity is related to the apparent dielectric permittivity of the medium between the access tubes. The permittivity is then related to the volumetric water content using a standard calibration (Topp et al., 1980) or a medium specific calibration (e.g., Masbruch and Ferré, 2003).
Generally, the EM wave velocity is determined with the explicit assumption that the first-arriving energy travels along a direct path from the transmitter to the receiver. Several investigators have pointed out that if refracted waves are first to arrive, this assumption can give rise to erroneous velocity measurements (Wang and McMechan, 2002; Hammon et al. 2002; Rucker and Ferré, 2003). For example, Rucker and Ferré (2003) demonstrated that critical refraction at the air-ground interface can result in an underestimation of the water content of the shallow subsurface. Critical refraction at the air-soil interface is a specific example of a more general case. Critical refraction of EM waves always occurs when an EM wave crosses a boundary from a low-velocity layer into a higher velocity layer (at any angle other than normal to the boundary). For example, soil profiles commonly have high water content layers adjacent to low water content layers. As BGPR antennae are lowered within the high water content layer toward the lower water content layer, critically refracted energy begins to arrive before energy traveling directly through the high water content layer (Fig. 1) . Similarly, if the antennae are lowered from a low water content layer into a high water content layer, critically refracted waves will continue to arrive before direct waves until the antennae are some distance below the layer boundary. The vertical distance from the boundary within which critically refracted arrivals are first to arrive has been referred to as the refraction termination depth, zrtd (Rucker and Ferré, 2003).
In this investigation, we extend earlier work (Rucker and Ferré, 2003), which identified the potential for BGPR water content measurement errors near the ground surface due to critical refraction. Specifically, we examine the effects of critical refraction throughout a discretely layered system. The effects of first-arriving critical refractions are quantified by comparing the total length of water in the soil profile, determined assuming that all travel paths are direct, with the known length of water in the profile. In this paper, we define the total length of water (Lw) as the integrated volumetric water content over the entire length of the profile, that is 1where zL is the total length of the profile. We then show how critical refractions can be identified and accounted for under some conditions. The specific objectives of this investigation are
To develop criteria to distinguish direct arrivals from critically refracted arrivals on a BGPR travel time profile
To demonstrate how to analyze BGPR travel time profiles to account for critically refracted arrivals
To identify those conditions that are not amenable to this correction
To show, theoretically and experimentally, the consequences of not considering critically refracted arrivals when analyzing BGPR travel time profiles
The travel paths of acoustic and electromagnetic energy in a layered earth structure can be described by assuming that the wave travels in a straight line through a medium and then reflects or refracts at boundaries according to Snell's Law (Telford et al., 1990). Examples of reflected and refracted arrivals on radargrams can be seen in Young and Sun (1999) and Fisher et al. (1992). Because ZOP BGPR analysis generally makes use of only the first-arriving energy, correct velocity analysis requires knowledge of the travel path of the first-arriving energy. It is impossible to distinguish these events immediately from direct observation of a recorded trace. It is only through the simultaneous analysis of measurements made at various offset distances or depths that different travel paths can be identified. For example, reflected and diffracted arrivals show distinct hyperbolic patterns with depth on a radargram. The travel time of a critically refracted wave increases linearly with depth.
Ray tracing analyses (Cai and McMechan, 1999) show that reflected arrivals will always arrive later than a critically refracted or direct arrival when conducting measurements with BGPR in ZOP mode (Ellefsen, 1999). Therefore, reflected paths can be ignored. The travel time of critically refracted waves can be shorter than that of a direct wave if the measurement is made close to the boundary between two layers of contrasting dielectric permittivity. The refraction termination depth, zrtd, defines the distance from the boundary within which the critically refracted wave arrives before the direct wave (Rucker and Ferré, 2003). The value of zrtd depends on the velocity contrast between the layers and the antennae separation distance.
Although it is generally less accurate than determining subsurface velocities through numerical modeling of the wave equation (Ellefsen, 1999; Hollinger and Bergmann, 2002) or wavefield extrapolations (Cai et al., 1996), ray tracing analysis is commonly applied to ZOP BGPR because it is much more simple and less computationally demanding. One major simplifying assumption in conducting a ray-trace analysis of ZOP BGPR profiles is that there are no lateral changes in dielectric permittivity between the boreholes. This is appropriate if the changes in water content are larger than one-quarter wavelength of the propagating EM wave (Chan and Knight, 1999, 2001; Schaap et al., 2003). (The wavelength of a 100-MHz signal in the air is approximately 3 m, decreasing to as little as 0.3 m in wet soil.) The following analysis assumes that there are no lateral changes in water content; the analysis only considers vertical changes in water content across distinct boundaries. Zero offset profiling first-arrival travel time profiles do not contain enough information to identify lateral changes in water content. If strong lateral changes are expected, then MOG measurements should be performed.
CONDITIONS GIVING RISE TO FIRST-ARRIVING DIRECT WAVES
The travel time of a direct wave through a high water content (low-velocity) layer is 2where vlow is the velocity of propagation through a low-velocity layer and x is the antennae separation distance. The travel time of a critically refracted wave is 3where vhigh (Fig. 1) is the velocity of propagation through an immediately adjacent higher velocity layer, z is the vertical distance from the antennae measurement depth to the interface between the two layers, and ic is the critical angle. Critically refracted energy arrives at the receiver through secondary waves generated at the boundary between the layers (Parkin et al., 2000), where a head wave is created. The critical angle depends entirely on the velocities of the layers: 4
The critically refracted energy will always be the first-arriving energy if z is small. That is, as z approaches zero, trefr approaches x/vhigh, which is smaller than tdirect. As z increases, trefr increases linearly. The refraction termination depth is defined such that, for z > zrtd, in a simple homogeneous half-space, direct waves will arrive before critically refracted waves: 5
In Rucker and Ferré (2003) the concept of a refraction termination depth was applied to the determination of near surface water contents. However, Eq.  shows that any two adjacent layers with differing velocities will give rise to a zrtd. In fact, a subsurface layer that is surrounded by higher velocity layers will have two ztrd values, associated with the boundaries above and below the layer. For example, consider a clay layer surrounded by sand layers. Borehole ground penetrating radar measurements conducted within the clay will refract critically along the upper and lower boundaries, giving rise to two refraction termination depths. Similarly, during the advance of a wetting front with no surface ponding, measurements made in the wetted region may refract at the ground surface or at the wetting front. In either case, if the thickness of the low-velocity layer is larger than the sum of the two zrtd values, then a direct arrival will occur within the low-velocity layer. If, on the other hand, the low-velocity layer is thinner than this sum, then the direct arrival will never be the first to arrive and first-arriving critical refraction along the upper boundary will give way immediately to first-arriving critical refraction along the lower boundary. Therefore, the minimum thickness of a low-velocity layer that will give rise to direct arrivals, hmin, is 6awhere v1 and v3 are the velocities of the overlying and underlying layers, and v2 is the velocity of the low-velocity layer. Equation [6a] can also be expressed in terms of water content: 6bwhere the velocities of Eq. [6a] were replaced by water contents derived from a linearized form of the Topp et al. (1980) equation (Rucker and Ferré, 2003): 7
Consider two three-layer examples (Fig. 2) . The thickness of the middle layer is 2.5 m for Case A and 1.75 m for Case B. The top layer in both cases represents an extremely dry layer with θ1 = 0.05 cm3 cm−3. The middle soil layer has a volumetric water content, θ2, of 0.52 cm3 cm−3. The lowest layer is a semi-infinite half-space with a θ3 of 0.17 cm3 cm−3. The travel times of both critically refracted and direct waves are shown to the right in Fig. 2. The thick black line represents the first-arriving energy on a radargram with a vertical sampling interval, dz, of 0.25 m. The thin dashed green lines, marked with trefr1 and trefr2, are the critically refracted travel times from above and below the high water content layer, respectively. The thin blue line represents the direct arrival time for the high water content layer. For Case A, the zrtd distances from above and below the layer do not overlap within the low-velocity layer. As a result, direct arrivals appear first for some section of the travel time profile within the high water content layer. However, for Case B direct arrivals within the high water content layer are never first arriving. This is consistent with Eq. , which predicts that the minimum thickness of a low-velocity layer that would have first-arriving direct arrival is 1.927 m. For Case B the low-velocity layer is hidden and its water content will be underestimated.
RESOLVING VELOCITIES OF THIN LAYERS
Alluvial or lucastrine deposits commonly exhibit a layered structure, whereby high-velocity coarse-grained material contains lenses or continuous layers of finer grained material. The fine-grained layers can have a thickness on the order of a few centimeters to a meter or more (Boggs, 1995; Zheng and Gorelick, 2003). A repeating, alternating layered model (ALM) of high–low–high velocity can be used to approximate this type of layered structure for ZOP BGPR analysis. The first high-velocity layer for the ALM is air, with a known velocity of 0.3 (m ns−1). The middle, low-velocity layer represents the uppermost sediments near the ground surface and may contain several distinct layers of low-velocity material. These low-velocity sediments can be grouped into a single low-velocity layer for this analysis. At some depth below the low-velocity layer, a high-velocity layer will be intersected by the BGPR. The high-velocity layer is easily identified by its lower travel time compared with the travel times of the layer above (Fig. 1 and 2). The ALM will begin again at the second high-velocity layer, and can repeat throughout the profile by considering the lowest layer as the uppermost layer of the next high–low–high ALM. Once the pattern of layering is established, the travel time profile can be used to find areas where critical refraction is occurring.
Rucker and Ferré (2003) identified two methods to obtain the near-surface velocity of propagation. The more robust method relied on the slope of the travel time profile: 8where vsoil (m ns−1) is the EM velocity in soil, vair (m ns−1) is the EM velocity in air, and Δ (ns m−1) is the inverse slope of the travel time versus depth profile, dt/dz. Equation  can be used to calculate the soil velocity of deeper sediments, where vair is replaced by vhigh (v1 or v3 in Fig. 2). Equation  can also be applied to the slopes approaching either the upper or lower boundary. Note that although the slope approaching the lower boundary is negative, Eq.  depends on the square of the slope, so the sign of the slope is irrelevant.
The first step in resolving the soil velocity in areas of critical refraction using the slope analysis is to identify those travel times that could be from a direct arrival (Fig. 3a , Step 2). The inversion of the direct arrival travel time, with a known antennae separation, is used as vhigh in Eq. . Two criteria for identifying direct arrivals are (i) travel times that are lower than the travel time above and below the measurement point or (ii) equal travel times for two or more adjacent measurements. The first criterion positively identifies a direct arrival because critically refracted waves can only originate from within a low-velocity layer that is adjacent to a high-velocity layer. The second criterion may incorrectly identify a direct arrival. For example, the two largest travel times shown in Case B of Fig. 2 appear to be equal (within the accuracy of travel time measurements) and therefore could be from a direct arrival. However, both of these first arrivals are from critically refracted first arrivals. Other adjacent first arrival travel times that are equal within this figure are direct arrivals.
The slope equation (Eq. ; Fig. 3a, Step 3) is applied to the remaining travel times, which may be associated with critical refractions. The inversion of the travel time produces a first-guess reconstructed velocity profile (Step 4). Then, at each depth, the velocity inferred assuming that the wave followed a critically refracted path, is compared with the velocity assuming that the first arrival is direct (Step 5). The slower of these velocities is the correct velocity (Step 6). That is, if the inferred critically refracted velocity is higher than the direct velocity, then the measurement is a direct arrival. If it were not, then the velocity required to give rise to a critically refracted first arrival in the measured time would give rise to an even faster direct arrival. In this manner, some measurements will be reclassified as direct arrivals. Finally, some of the arrivals classified as critically refracted first arrivals will be reclassified as direct if they are farther than zrtd from the nearest boundary, based on the inferred velocity profile. Once the measurements have been classified as direct or critically refracted, the water content profile can be determined using Eq.  (Step 7). Figure 3b shows how a spreadsheet may be used to help in the reconstruction of the velocity profile. Note that Measurement 9, located at 2 m below ground surface was originally classified as a critically refracted arrival. However, after comparison with the direct arrival, it was reclassified.
To demonstrate the reconstruction of a velocity profile using this approach, consider two soil profiles, a coarse profile in which all layers have dimensions that are >2dz (Fig. 4) and a fine profile that contains a few layers with dimensions that are <2dz (Fig. 5) . Five soil types are distributed throughout the profile. The van Genuchten water retention functions (van Genuchten, 1980) for these soils were defined based on the database of Carsel and Parrish (1988). The two profiles are similar, except for the section from the 2.85- to 5.4-m depth. In the coarse profile, this middle section contains three soil layers, whereas the fine soil profile contains seven layers. The profile is taken to be under hydrostatic conditions, above the water table. The direct wave velocity profiles were defined based on the water content profiles derived using Eq. , and the first-arriving velocity profile was determined using raypath analysis. To simulate a typical ZOP BGPR survey, the first-arriving travel time profile was sampled with a dz of 0.25 m.
The direct arrival at the soil surface can be identified easily because it has a velocity equal to 0.3 m ns−1. Other possible direct arrivals, based on the two criteria listed above, are marked on the travel time profiles with filled yellow squares. Slopes were fitted through linearly increasing (or decreasing) sections of the travel time profile to calculate vsoil (or vlow) to produce the first-guess velocity profile. Some examples of slope values obtained by fitting are given on the travel time profiles. The direct-only velocity profile and first-guess velocity profile were compared, and the lowest velocity was taken at each depth. Based on this final velocity profile, the water content profile was calculated using Eq. .
The calculated water content profiles for the coarse and fine layered profiles are plotted on Fig. 6 . The thin solid blue line shows the actual water content profile. The dashed green line shows the water content profile assuming all first arrivals are direct. The heavy red line is the estimated water content profile based on the final velocity profile that considers critically refracted first arrivals. The actual total length of water (Lw) in the profile was 1.11 m for the coarse profile and 0.86 m for the fine profile. Assuming all first arrivals are direct, the inferred Lw values were 0.78 and 0.51 m, giving relative percentage errors of Lw of 29.5 and 41%. The RMSE in water content is 0.079 and 0.10 cm3 cm−3 for the coarse and fine profiles, respectively. The relative percentage errors reduced significantly to 5.4 and 15% after we accounted for critical refractions, with RMSE values reducing to 0.016 and 0.059 cm3 cm−3. The larger errors seen for the fine profile are due to the relatively large ratio of sampling discretization to the average layer thickness. One layer in particular, a clay loam located at the 3.75-m depth with a thickness less than the vertical sampling interval, was not sampled. A clay loam layer located at the 2.85-m depth was sampled, but its water content was grossly underestimated.
Borehole ground penetrating radar travel time measurements were collected during a field-scale infiltration experiment (Rucker and Ferré, 2003). The infiltration experiment ran for 66 h, with a total of 0.898 m of water applied to an area of 25 m2. A pulseEkko1000 (Sensors and Software, 1999) BGPR system with 100-MHz antennae was used to conduct the measurements in two boreholes spaced 3.1 m apart. The depth of the profile was 15 m, and the vertical sampling interval was 0.25 m. The travel time profiles measured before and after infiltration (Fig. 7) show the effects of added water. The travel time increased significantly in the top 5 m. The bottom scale of Fig. 7 shows the calculated water content using a linearized form of the Topp equation (Ferré et al., 1996) assuming that all first arrivals are direct. Note that the unreasonably low water contents near the ground surface are indicative of critically refracted first arrivals (Rucker and Ferré, 2003).
Assuming that all arrivals are direct, the calculated values of Lw before and after infiltration measurements were 0.59 and 1.17 m, respectively. The difference suggests that 0.58 m of water was applied at the surface, which is an underestimation of 35% compared with the known length of applied water. The underestimated Lw is in reasonable agreement with the 38% underestimation of Lw reported by Rucker and Ferré (2002) on the basis of a comparison of BGPR with the known amount of applied water and with measurements made with time domain reflectometry and with a neutron probe in a similar infiltration experiment. The profile collected before infiltration has two regions that can be positively identified as having critically refracted first arrivals. These regions are associated with the ground surface (0–1.75 m) and with a clay-rich layer (8–9.5 m depth). After infiltration, critically refracted first arrivals are also seen at the wetting front (5–5.75 m depth). After correcting for the effects of critical refraction, the water content profiles in the shallowest 6 m, both before and after infiltration (Fig. 8) , differ from those determined assuming that all first arrivals are direct (Fig. 7). These differences are confined to the regions near the ground surface and near the wetting front (shown with gray shading on Fig. 8). The change in Lw calculated from the corrected profiles is 0.81 m, which is in good agreement with the applied Lw of 0.898 m. The correction for critical refraction reduced the measurement error from 35 to 9.8%. This 9.8% difference between calculated change in Lw and known applied length could be due to either radial flow outside of the infiltration gallery, thin high water content layers that could not be properly assessed, or both. Rucker and Ferré (2003) showed that a TDR located 0.5 m from the north boundary of the infiltration gallery experienced an increase in water content.
In addition to the underestimation of the length of applied water, unidentified critically refracted first arrivals can have other hydrologic implications. For example, the velocity of the wetting front is underestimated because measurements made with the antennae within a distance of zrtd above the wetting front will have critically refracted first arrivals. As a result, the wetting front appears to be more shallow than it actually is at any given time (Fig. 8). For this example, the velocity of the wetting front was calculated based on the rate of movement of the position of a reference water content (Warrick, 2003) of 0.18 cm3 cm−3. These positions are shown in Fig. 8 as black circles. Without considering critical refraction, the position of the wetting front is 5.17 m. After correction for critical refraction, the wetting front position is 5.63 m. The velocities of the wetting front before and after correction are 2.35 × 10−5 and 2.56 × 10−5 m s−1, respectively. Similarly, a sharp wetting front may appear more smoothed due to the effects of critical refraction. This could give rise to overestimation of the effects of capillarity at the wetting front.
DISCUSSION AND CONCLUSIONS
Electromagnetic waves travel along several paths through a layered subsurface. Correct interpretation of the velocity of propagation requires that the path of the first-arriving energy be identified. In this investigation, we consider the effects of first-arriving critically refracted waves on the interpretation of the water content profile from BGPR measurements made in ZOP mode.
Critically refracted first arrivals may occur whenever BGPR antennae are lowered through a high water content layer that is adjacent to a lower water content layer. Although critical refractions cannot be distinguished from direct arrivals on a single BGPR trace, they exhibit identifiable behavior on a travel time profile. Specifically, the travel time of first-arriving critical refractions will decrease linearly as the antennae approach the boundary with the lower water content region. We describe a method whereby first arrivals can be classified as arising from either direct or critically refracted waves. Once classified, the appropriate relationship can be used to determine the velocity, and therefore the water content, from the travel time at each depth. This correction can greatly improve the accuracy of water content profiles in layered systems. However, errors persist if high water content layers are thin compared with the sampling depth interval.
This material is based on work supported by the National Science Foundation under Grant no. 0097171. Additional support was provided by funds from NASA under contract NASA/JPL 1236728.
- Received May 2, 2003.