Nondestructive and noninvasive visualization and quantification of dynamic subsurface hydrologic processes are needed. Using a ground-penetrating-radar (GPR) antenna array, time-lapse common-offset gather and common midpoint (CMP) data can be collected by fixing the antenna at a given location when scanning the subsurface. We have determined wetting front depths continuously during a field infiltration experiment by estimating electromagnetic (EM) wave velocities at given elapsed times using GPR antenna array data. A surface GPR antenna array system, consisting of 10 transmitters (Tx) and 11 receivers (Rx), that can scan each Tx-Rx combination in 10 s at a millisecond scale has been used to acquire all 110 Rx-Tx combinations in approximately 1.5 s. The field infiltration experiment has been conducted at an experimental field near the Tottori Sand Dunes in Japan. Using the estimated EM wave velocity from the CMP data, the depth to the wetting front is computed every minute. The estimated wetting front arrival time agrees with the time at which a sudden increase in the moisture sensor output is observed at a depth from 20 cm and below. We determine that time-lapsed CMP data collected with the GPR antenna array system could be used to estimate EM wave velocities continuously during the infiltration. The GPR antenna array is capable of accurate and quantitative tracking of the wetting front.

A deeper understanding of subsurface hydrologic and physical processes in the vadose zone is important for many applications in agricultural and environmental disciplines, such as optimal water usage for crop cultivation, effective removal of soil contaminants, or accurate predictions of water and heat exchange at the soil-atmosphere interface. To accurately characterize and predict water flow in the unsaturated zone, in situ hydrologic parameters, such as the hydraulic conductivity and water retention curves, must be determined. Although these parameters are traditionally measured using soil samples collected from experimental field sites, recent developments in geophysical techniques have allowed the monitoring of hydrologic processes in the unsaturated zone without disturbance. Once quantitative information regarding subsurface hydrologic processes is derived from geophysical measurements, such as temporal changes in soil water content, hydrologic models can be used or coupled with geophysical models to estimate soil hydraulic parameters (e.g., Lambot et al., 2004; Kowalsky et al., 2005; Busch et al., 2013; Léger et al., 2014). Among the numerous geophysical techniques used, ground-penetrating radar (GPR) has been the most widely used method for monitoring and characterizing soil water dynamics because the soil water content directly characterizes the propagation of electromagnetic (EM) waves emitted and received by GPR antennae (Huisman et al., 2003; Lambot et al., 2006; Binley et al., 2015; Vereecken et al., 2015). When conducting analyses in the time domain, the traveltime of the EM wave from the transmitter to the receiver is measured with a variety of configurations, such as cross-borehole GPR, off-ground GPR, and on-ground surface GPR. On-ground surface GPR (referred to as surface GPR in the remainder of this paper) is the most suitable for scanning in situ unsaturated zones because it does not require the drilling of boreholes (unlike cross-borehole GPR) and it is less sensitive to surface conditions compared to off-ground GPR. To collect common-offset gathers (COG) for subsurface visualization, surface GPR is often used with a constant antenna separation (offset). Surface GPR also can be used with different antenna offsets to collect multioffset gathers (MOGs). This feature is particularly attractive when a subsurface EM wave velocity structure must be determined to further convert traveltime data to depth.

Surface GPR has been used to visualize and quantify field infiltration processes. In a pioneering study, Vellidis et al. (1990) use surface GPR to visualize wetting front movement in sandy soil during irrigation. In this study, in addition to obtaining COG, a known buried object was used to determine the average EM wave velocity and wetting front depth. Trinks et al. (2001) acquire time-lapsed COG (every 30 min) using surface GPR to monitor the movement of water in a sandbox. The average EM wave velocity was estimated using known depths to the bottom of the sandbox and to electrodes buried prior to the experiment. Saintenoy et al. (2008) and Léger et al. (2014) use water injected from a known fixed depth and maintain the GPR antennae at a fixed location on the surface, such that the distance between the antenna and injection point remained constant during the measurement. This allowed computation of an average EM wave velocity and therefore an estimate of the average dielectric constant. The soil dielectric constants can be later converted to soil water content using petrophysical relationships (e.g., Topp et al., 1980). However, this is not always applicable, especially when analyzing soil water dynamics in in situ soils because there is not always a reflector at a known depth for velocity estimation.

None of these referenced approaches directly determines the location of the wetting front or the EM wave velocity structure during infiltration processes. Because the wetting front dynamically and continuously changes its location (or depth) during infiltration, time-lapse MOG must be collected in a sufficiently short time interval compared to the speed at which the wetting front moves. With common bistatic surface GPR antennae, collecting time-lapse MOG takes time because it requires physically moving the antennae for each antenna offset needed. To expedite the MOG data collection procedure, some studies have used multiple antennae. For example, Pan et al. (2012) use four antennae (two transmitting and two receiving) to sequentially collect MOG from four antenna combinations before and after a heavy rainfall. With this configuration, three different antenna offsets could be scanned simultaneously. Common midpoint (CMP) data could then be obtained using three antenna offsets for velocity analysis, although the use of only three different offsets is not sufficient for such an analysis. Although multiple antennae allow for collection of MOG with a limited number of antenna offsets over short time intervals, the limited number of antennae makes it difficult to accurately estimate the EM wave velocity using a standard velocity analysis. Mangel et al. (2012) develop an automated system attached to a sandbox to collect wide-angle reflection and refraction (WARR) data, that is, one MOG with 21 traces in every 30 s. In their system, the transmitter was placed at a fixed location, whereas the receiver was mounted to a carriage that moved along a survey line with a motor. Normal moveout (NMO) analysis was conducted for the WARR data to estimate the EM wave velocity and depth to the reflector, which was the wetting front. With the automated data collection system, Mangel et al. (2020) could collect a single CMP data set consisting of 84 traces in 1.8 s. Although these studies showed promising results, the approach was limited to the sandboxes in the lab and was not applied in the field.

Recently, a family of GPR array (e.g., Eide et al., 2014) and a system that allows for very fast WARR-type data collection in a millisecond range (e.g., Kaufmann et al., 2000; Diamanti et al., 2018) have been developed. Essentially, those systems electrically and sequentially switch among predefined antenna combinations after 10 and several milliseconds to scan subsurface processes with minimal effort. This GPR array has mainly been used to collect COG to generate 3D subsurface images, which are useful for infrastructure inspections (Eide and Hjelmstad, 2002), archaeological surveys (Linford et al., 2010), and unexploded ordnance (UXO) and landmine detection (Eide and Hjelmstad, 2004; Sato et al., 2004). By placing such GPR systems at a fixed location, time-lapse data can easily be collected to visualize dynamic subsurface processes, such as infiltration. Time-lapse radargrams collected every 1.5 s with the GPR array system show that the reflection signals from the wetting front are gradually delayed as the wetting front moved downward during the forced vertical infiltration experiment at the Tottori Sand Dunes test site in Japan (Iwasaki et al., 2016). However, CMP data collected with the array GPR system have not been fully investigated for velocity analysis. The main objective of this study was to determine the wetting front depth at a given elapsed time using the radargrams collected with the GPR array and by continuously estimating the EM wave velocity structure during the infiltration experiment. To achieve this objective, time-lapse CMP data were analyzed to obtain the EM wave velocity structure. Then, the depth of the wetting front was computed for the elapsed time from the start of water injection.

GPR array system

In this study, we used a product from 3D-Radar (Eide et al., 2014) as a GPR array system consisting of a linearly aligned array of 10 transmitting (Tx0–Tx9) and 11 receiving (Rx0–Rx10) bow-tie monopole antennae as shown in Figure 1 (DXG1820 antenna). Unlike conventional impulse GPR systems, the GPR array system is a step-frequency continuous waveform radar system operating within the frequency range of 100–3000 MHz. Frequency-domain data were converted to time-domain data using an inverse Fourier transform. Additionally, data processing algorithms, such as bandwidth filtering and gain, can be applied at the same time to enhance data. In this study, the inverse Fourier transform with a window function of Tukey was used to a frequency range of 50–3000 MHz. Then, a gain function was applied to amplify the signals with depth. All of these processes were performed with a software called Examiner Data Processing provided by 3D-Radar. Based on the recommendation for the 3D-Radar, the antenna offset was determined by the actual distance between the tips of the V-shape Rx and Tx antenna, where the vertical offset is 85 mm as shown in Figure 1. The shape of the antenna has been designed to maximize the performance over the frequency band of the GPR system (Eide et al., 2014). According to the manufacturer, the phase center of an antenna element, that is, the theoretical point at which the signal is effectively transmitted and received, is geometrically coincident with its feed point, which is located at the tip of the V-shape. Because the Tx and Rx were not inline in the GPR array system used in this study, the actual separation distances between Tx and Rx were computed for further velocity analysis. Although we assumed that the vertical distance of 85 mm was sufficiently small such that accounting for any 3D effects in the analysis was unnecessary, it is necessary to further investigate whether this assumption is valid in the future.

The GPR array system electronically and sequentially switched predetermined combinations of Rx-Tx antennae using radio frequency multiplexers (Eide et al., 2014). The maximum number of antenna combinations for this particular array GPR system was 110 because there were 10 Txs and 11 Rxs. In this study, we set up the GPR system to scan 110 antenna combinations by sequentially switching from Tx0 to Tx9 and each Tx was combined with 11 Rx from Rx0 to Rx10. Scanning all 110 antenna combinations requires less than 1.5 s. After conversion to time-domain data, some scanned data were extracted to reconstruct a data set necessary for further analysis. For example, Figure 2a depicts antenna combinations used to create COG because all antenna combinations have the same antenna separation distance of 113 mm. The open dots in Figure 2 indicate the central locations of the Tx and Rx tips that are the measuring points. With this setup, at each measurement, a radargram with 20 Rx-Tx data can be generated, which is equivalent to a 1.425 m COG radargram with a 0.075 m antenna step size. By generating the COG radargram every 1.5 s, a 3D space-time radargram can be generated. When Tx is fixed and combined with all Rxs, WARR, one of the MOG data, can be generated because the antenna offset varies depending on the Rx location. Figure 2b depicts one of the WARR setups, where Tx0 is coupled with Rx0 to Rx10. Again, the open dots indicate central locations between the Tx and Rx tips. Figure 2c shows antenna combinations with common central locations again represented by an open dot. This type of antenna setup is often used for CMP measurements. CMP data can be further used for EM wave velocity analysis. With a conventional bistatic pulse GPR system, obtaining CMP data requires physically moving the antennae for every antenna separation, which can take several minutes depending on the antenna step size. Mangel et al. (2012) develop a system that allowed for CMP data to be obtained in 30 s by mounting the Tx and Rx to a steel rail. However, although this may be effective in the lab, using such as a system in the field would be difficult. To improve the signal-to-noise ratio and obtain a smooth hyperbolic reflection signal, a small antenna step size with hundreds of measurements is usually necessary. With the GPR array system, there are only 10 antenna combinations, with an increase in the antenna separation distance from 113 to 1425 mm. Although a constant antenna step size is used for typical CMP measurement, due to the Tx and Rx alignment (Figure 1), the antenna step size for CMP with the GPR array system shown in Figure 2c is not constant. Depending upon the target depth, the maximum antenna separation should be significantly larger than the depth of investigation needed, such that a clear hyperbola shape can be observed for better velocity analysis. However, the CMP data from the GPR array are limited to the maximum separation distance of 1425 mm. The accuracy of the velocity estimation based on the CMP data collected with the GPR array system may be therefore affected by these shortcomings. With the GPR array system, we can obtain CMP data every 1.5 s, which allows for the application of this technique to monitor subsurface dynamic hydrologic processes, such as water flow in soils (e.g., infiltration).

Infiltration experiment

The infiltration experiment was conducted in a large greenhouse facility at the Arid Land Research Center of Tottori University in August 2015 (Iwasaki et al., 2016). The experimental site was located near the Tottori Sand Dunes (Figure 3) with more than 30 m of sand deposition, resulting in a uniform sand dune layer. The soil surface inside the greenhouse was kept bare during the experiment. Because the water table depth at the site was approximately 7 m, any influence from capillary rise in the groundwater could be ignored during the experiment. During this infiltration experiment, to ensure uniform infiltration beneath the antennae, water was supplied over an area larger than the size of the array antenna using six porous tubes, each 250 cm in length. The porous tube made from rubber could seep out water from the entire tube surface when pressure was applied to water inside the tubes. These tubes are commercially available and are commonly used for water-saving irrigation in horticulture. They were placed parallel to one another on the surface with approximately 15 cm distance between them (Figure 4). The tubes were connected to one inlet to control the flow rate. The GPR array system was then placed on top of a thin wooden plate slightly larger than the antenna covering all of the porous tubes (see Figure 4). A preliminary test confirmed that the effect of the thin wood plate on receiving reflected signal was minimal. After placement of the array antenna, all predefined 110 antenna combinations were scanned continuously for 6 h, until the end of the experiment. Unlike a typical GPR survey in which the antennae are moved during measurement, in this study, the position of the array antenna was fixed to ensure data reproducibility. Prior to applying water, the reproducibility of the data was tested by collecting time-lapse background data, which showed no artifacts or variations. Water was supplied through the porous tubes at a rate of 7000  cm3/min from the inlet for 4 h equivalent to an infiltration intensity of 22.4  cmh1.

A 100 cm probe-type capacitance soil moisture sensor (Profile Probe [PP], Delta-T Devices) was installed adjacent to the antenna to independently monitor changes in the soil dielectric constants during the infiltration experiment, corresponding to the depths at which the moisture sensors were mounted (Figure 4). The sensors were mounted at 10, 20, 30, 40, 60, and 100 cm from the top, referred to as PP10, PP20, PP30, PP40, PP60, and PP100, respectively, in the remainder of this paper. The PP has been successfully used in previous research at the Tottori Sand Dunes to accurately measure the soil moisture content of the dune sand (e.g., Inoue et al., 2008).

Velocity analysis and wetting front depth estimation

To convert the EM wave traveltime data into the depth data, a velocity structure or a dielectric permittivity structure needs to be determined. Although there are a variety of approaches, including full-waveform inversion (e.g., Busch et al., 2012), many studies using surface GPR systems used CMP data to construct such profiles to determine soil water content profiles through the NMO analysis (e.g., Mangel et al., 2012; Steelman et al., 2012). One common approach, known as semblance, constructs a velocity spectrum by stacking CMP traces after NMO correction (Yilmaz, 2001). In the semblance approach, the velocity spectrum over a range of possible velocities depicts stacked amplitudes for the two-way traveltime (TWT). Amplitude peaks corresponding to the optimum average velocities where reflections occurred in the CMP data are then selected to determine the apparent root-mean-square velocities (vrms), assuming horizontal reflectivity. Then, vrms was used to determine the interval velocity of a given layer using the Dix (1955) equation. To reliably and accurately estimate the velocity section from the CMP data, the data acquisition spacing needs to be less than half of the shortest wavelength to prevent aliasing (the Nyquist criterion). However, the CMP data collected with the GPR array system used in this study clearly violated this criterion because an average spatial step size was 146 mm. Such sparsity can be overcome by interpolating sparsely collected CMP data, for example, with a Fourier domain reconstruction method (Yi et al., 2015).

In this study, because of severe sparsity in the CMP data, a common semblance approach was not used for the velocity estimation. Because there was only one clearly observed target reflection plane, a different approach was used. For the CMP data, the TWT tc from a single given reflection plane at depth d0 is a function of the antenna separation x as follows:
where vr is the root-mean-square velocity to the reflection plane and t0 is the TWT at zero offset determined from the CMP radargram. On the CMP radargram, the TWT draws a hyperbolic-shaped curve. In this study, vr was estimated by visually fitting equation 1 to the hyperbolic curve in the CMP radargram. Before fitting the curve, time zero was corrected by fitting the airwave traveltime to the CMP radargram. Using the estimated vr, the depth to the reflection plane (the wetting front) d0 could be computed. Although this approach cannot provide any measure of uncertainty to the estimated velocity profile, some uncertainties can be evaluated when determining interval velocities. In this study, the CMP data for every 1 min (60 s) were used for further velocity analysis.

Time-lapse array GPR radargrams

Figure 5a shows a time-lapse panel diagram of the 110 Tx-Rx combinations acquired every 1.5 s during the first 60 min after the start of water injection in the infiltration experiment. The acquired frequency-domain data were converted to the time domain by inverse Fourier transform. At the same time, high- and low-pass filters were applied with a depth-dependent gain. We assumed that the frequency dependence of the dielectric permittivity is negligible for our study because of the frequency range covered by the GPR system and the type of soil being investigated. Although the frequencies where the effect is stronger are those greater than 1 GHz, the dominant frequencies of the GPR array system are less than 1 GHz. Tottori dune sands are known to have low electrical conductivity. Therefore, the attenuation was expected to be small (the imaginary part of the dielectric permittivity). The panel diagram can be divided into 10 sections depending upon the Tx used. The location of the Tx is depicted by a white triangle. Although some noise in the data was observed especially at early traveltimes, as the time elapsed from the start of the water injection, a strong reflection signal appeared that was delayed with the elapsed time. This originated from a reflection point or plane that changed its location with time, which was the wetting front. From this panel diagram, the COG or CMP can be easily reconstructed.

Figure 5b shows COG radargrams at the elapsed time equal to 0 (te=0  min), 25, and 50 min reconstructed from MOG shown in Figure 5a. Recall that the COG offset was 113 mm. The horizontal axis of the COG radargram in Figure 5b indicates the position of the central point between Rx and Tx (Figure 2), where the horizontal origin corresponds to the location of Rx5. In the initial COG radargram (te=0  s), in addition to the strong direct airwave signals, there are some continuous reflections at traveltimes of 10, 20, and 35 ns, which may be due to layering or differences in soil moisture conditions. There is also a hyperbolic event at 15 ns TWT prior to water injection (te=0  min), which indicates a point reflection in the subsurface. The signal for this reflection started to fade as the time elapsed (te=25 and 50 min). We did not excavate the site to confirm these reflections after the experiment, so their exact cause remains unknown. However, these events did not affect any analyses performed in this study.

Although a reflection from the wetting front was not observed in the initial COG of Figure 5b, reflection signals were slightly delayed as the time elapsed. This was due to a decrease in the EM wave velocity near the surface as the infiltration was enforced. At an elapsed time te=25  min, a signal, different from the direct airwave signal, was observed at 10 ns TWT, which was the reflection from the wetting front. That signal further delayed to 15 ns TWT as the time elapsed to 50 min.

Velocity analysis and wetting front depth estimation

Because the reflection from the wetting front was clearly observed from te=5  min onward, we performed velocity analysis only on the CMP data obtained after te=5  min. For each minute of elapsed time te, the TWT given by equation 1 was visually fitted to the reflection signal from the wetting front by adjusting t0 and vr after correcting for time zero, which was done by computing the airwave with a velocity of 0.3mns1. Figure 6 shows the CMP data for every 5 min from te=5  min to te=50  min, along with the hyperbolic curve visually fitted as solid white curves. As for the airwave, instead of manually fitting a curve to the radargram, a theoretical value of 0.3mns1 was used to compute the traveltime as a function of the separation distance (the dashed white line in Figure 6). All curves were fitted to the positive peaks (the white band) of the reflectors. All curves fit well to the observed reflection curves shown in the CMP radargrams at all separation distances, indicating that t0 and vr effectively represent the field status during the experiment. Figure 7 exhibits the TWT for the antenna separation distance corresponding to 0 m as a function of the elapsed time. For both plots, the color indicates the values of the estimated root-mean-square velocity vr. The curves delay smoothly with the elapsed time, whereas the estimated vr decreases with time.

The estimated GPR wave velocity vr to the wetting front was used to compute its depth d0 at each given elapsed time te. Figure 8 shows the COG data for every 5 min from te=5  min to te=50  min, with TWT computed from vr using equation 1 with x=113  mm and estimated depths to the reflection, d0 depicted by the yellow dashed lines. At earlier times before te=20  min, although the direct airwaves and the reflected waves were clearly distinguished in the CMPs (Figure 6), the reflected waves in the COG data shown in Figure 8 are indistinguishable from the direct airwaves because the reflection point was too shallow for the given velocity. After the elapsed time of 20 min, the direct airwaves and the reflected waves in the COG data become distinguishable (Figure 8). Because the vr was obtained from the CMP, where the CMP is the center of the antenna (i.e., at the position of 0 m), the yellow dashed line generally matched the reflection wave at the position of 0 m. At other positions, the computed TWT underestimated the observed TWT, indicating that the wetting front was not completely flat. Considering the nature of heterogeneity in the hydraulic properties of field soils and naturally occurring preferential flow or unstable flow especially in sandy soil, it is not surprising that the reflection signal from the wetting front was not linear.

The estimated wetting-front depths were plotted as a function of te in Figure 9a, along with the PP output voltages as a function of te (Figure 9b). The depth to the wetting front was computed for the following two conditions: one in which the space with the surface porous tubes was completely wet (dielectric constant of the surface layer ε0=10.6) and the other in which it was completely dry (ε0=1.8). The depth to the wetting front lies between these two conditions, which define the upper and lower limits of the depth to the wetting front at each elapsed time. As shown in Figure 9, the estimated wetting front moved downward almost linearly with time, except when the wetting front slowed down between te=10 and 20 min. The wetting front reached a depth of almost 0.9 m after 90 min. In Figure 9b, changes in the PP sensor outputs recorded are plotted for six different depths (PP10 to PP100). Instead of converting the output in millivolts to the volumetric water content, the observed millivolts were directly plotted because we did not have a site-specific calibration for this conversion. The PP data provide quantitative information on how the dielectric constant (i.e., the moisture content) evolved over time during the experiment at each corresponding depth. Although the rates of increase were almost the same for all the depths, the constant values after the increase depended on the depths, where PP10 and PP20 showed slightly lower millivolts compared to the other three depths (PP30, PP40, and PP60). This indicates that basic soil characteristics, such as porosity or bulk density, may be different between soils at depths shallower than 20 cm and those at deeper than 20 cm. At depths of 30, 40, and 60 cm, the GPR estimated wetting front arrival time at that particular depth generally agreed with the time of the sudden increase in the PP reading regardless of the surface condition. The GPR estimated wetting front arrival time at depths of 10 and 20 cm was significantly faster than that demonstrated by PP. Although the reflection signals from the wetting front were distinctly separated from the direct airwave signals in the CMP radargrams during the early stages of infiltration (i.e., te=10, 15, and 20 min in Figure 6), the estimated depth disagreed with the PP data. This may be explained by the fact that, during the early stages of infiltration, the main driving force of the water movement was capillary force, whereas in later stages, the main driving force was gravity (Jury and Horton, 2004). Therefore, in the early stages of infiltration, the wetting zone may have evolved three dimensionally. Then, depending on where the sensors were installed or the location from which the EM wave was reflected, the resulting wetting front depth may be different. When the main driving force was gravity, water movement became more or less 1D. Overall, this study demonstrates that the GPR array system quantitatively estimates the wetting front depth during infiltration into the field soil in a noninvasive and nondestructive manner. Because these types of multioffset data are only recently becoming available in a fast and efficient manner, those obtained in this study constitute one of the first data sets actually showing how the wetting front evolves over time in the field soil. Although we collected data approximately every 1.5 s, only a part of the data was used, due to the fact we used the manual fitting approach, which is time-consuming, mainly because of data sparsity. In the future, more rigorous and fast velocity estimation approaches, that can be further developed for an automatic velocity estimation, need to be developed for sparse CMP data, for example, by coupling the Fourier domain reconstruct method (Yi et al., 2015) with the NMO analysis.

This study aims to determine wetting front depths continuously during a field infiltration experiment by estimating EM wave velocities at given elapsed times using GPR antenna array CMP data. This study reveals the following conclusions.

  1. A GPR antenna array system consisted of 10 Tx and 11 Rx successfully collected time-lapsed COG and CMP during a field infiltration experiment conducted at a Tottori Sand Dune site.

  2. EM wave velocities were estimated for every minute by visually fitting the hyperbola equation for the TWT to the CMP data acquired with the GPR antenna array during the field infiltration experiment.

  3. The depth of the wetting front computed with the estimated EM wave velocity agreed well with the data observed independently with a probe-type soil moisture sensor at a depth of 20 cm or below.

This study was financially supported by JSPS Grant-in-aid Scientific Research Program (20H03097 and 16H02580) and by Joint Research Program of Arid Land Research Center, Tottori University. We would like to thank the four anonymous reviewers and the associate editor for their constructive comments that helped to improve the manuscript. We would also like to thank Editage for the English language editing. Finally, we would like to express our sincere gratitude to GeoFive Co. Ltd. and 3D-RADAR AS for their technical support.

Data associated with this research are confidential and cannot be released.

Biographies and photographs of the authors are not available.

Freely available online through the SEG open-access option.