Distinguishing Capillary Fringe Reﬂection in a GPR Proﬁle for Precise Water Table Depth Estimation in a Boreal Podzolic Soil Field

: Relative permittivity and soil moisture are highly correlated; therefore, the top boundary of saturated soil gives strong reﬂections in ground-penetrating radar (GPR) proﬁles. Conventionally in shallow groundwater systems, the ﬁrst dominant reﬂection comes from the capillary fringe, followed by the actual water table. The objective of this study was to calibrate and validate a site-speciﬁc relationship between GPR-estimated depth to the capillary fringe (D CF ) and measured water table depth (WTD m ). Common midpoint (CMP) GPR surveys were carried out in order to estimate the average radar velocity, and common o ﬀ set (CO) surveys were carried out to map the water table variability in the 2017 and 2018 growing seasons. Also, GPR sampling volume geometry with radar velocities in di ﬀ erent soil layers was considered to support the CMP estimations. The regression model (R 2 = 0.9778) between D CF and WTD m , developed for the site in 2017, was validated using data from 2018. A regression analysis between D CF and WTD m for the two growing seasons suggested an average capillary height of 0.741 m (R 2 = 0.911, n = 16), which is compatible with the existing literature under similar soil conditions. The described method should be further developed over several growing seasons to encompass wider water table variability.


Introduction
Knowledge of the water table depth (WTD) is essential for understanding many environmental scenarios, including water management in agriculture. WTD, where hydraulic pressure equals atmospheric pressure, demarcates the saturated-unsaturated soil boundary. Naturally, a capillary fringe, which is quasi-saturated but still has negative water pressure, occurs above the water table [1]. The capillary fringe, sometimes called the transition zone, in the vadose zone mediates space, water, and nutrients for plants and soil organisms [2,3]. Both WTD and the depth to the top of the capillary fringe (D CF ) are subject to seasonal fluctuations, and thus, can affect agricultural water management practices throughout the growing season. WTD can be measured through a borehole, but this is generally unfeasible at larger scales, since it provides point scale measurements and is destructive.
High-resolution subsurface images of ground penetrating radar (GPR) can be used to measure shallow WTD, especially in coarse-grain soils [4][5][6]. Various GPR field techniques for water table studies have been developed over the decades. Annan et al. [7] stated that it was essential to have a sharp boundary between saturated and unsaturated zones in a GPR profile for precise WTD estimation. Larger wavelengths (i.e., lower frequencies) are recommended, even though the resolution of the radar profile decreases with increasing wavelength [7]. Loeffler and Bano [8] also found that GPR frequencies higher than 900 MHz do not identify the top of the saturated zone (water table) due to the effect from the capillary fringe and the transition zone. Therefore, most WTD studies have been carried out using GPR frequencies lower than 250 MHz [9]. Nakashima et al. [10] and Takeshita et al. [11] used common midpoint (CMP) data acquisition of GPR to explain the multiple reflections from the water table. CMP data can be used to estimate the variation of relative permittivity (ε r ), and therefore, to calculate the radar wave velocity (v) along the soil profile [10,12]. However, the common offset (CO) survey method using 100 MHz is the most commonly used method in water table studies [5].
The GPR technique was employed during pumping tests to measure the temporal fluctuation of WTD [13][14][15], and has been used for various groundwater studies [16][17][18][19][20][21][22][23][24][25][26][27][28][29]. However, it is challenging to interpret a shallow water table depth associated with closely spaced soil horizons with only the GPR outputs [30]. Advances in sophisticated sensor technology such as real-time WTD and soil moisture data would help to improve the GPR outputs [4]. Still, site-specific GPR data validation is needed to distinguish the accurate water table reflection from the capillary fringe or the transition zone. The objective of this study was to calibrate and validate a site-specific relationship between GPR-estimated D CF and measured WTD data in an agriculturally managed podzolic soil under boreal climate conditions. We carried out GPR surveys along a 42 m transect over two growing seasons through the use of a borehole with a real-time WTD sensor and soil moisture probes installed in shallow soil.
There are three characteristics to be considered when interpreting WTD from a GPR profile ( Figure 1). Firstly, saturation decreases from the bottom to the top of the capillary fringe (within the transition zone), consequently increasing the radar wave velocity [53]. The thickness of this zone depends on the amount, size (diameter), and interconnectivity of soil pores [4,54]. Secondly, oscillations of the reflected radar pulse due to the transition zone result in a series of bands representing the water table in the radar profile [4]. The top of the capillary fringe gives an earlier and more robust reflection than the actual water table [55,56]. The wetting front also has the same phase reflection as the water table [57]. Thirdly, the two-way travel time (TWTT) correspondence with the maximum absolute amplitude of the airwave (t air ) is the opposite of that of a reflection event (t reflect ) [55] (Figure 1). In addition to the above characteristics, the maximum absolute amplitude occurs at the second half of the respective GPR wavelet [58].
Equation (1) gives the radar signal velocity (v) of a nonmagnetic and low-loss geological material: where ε r is the relative permittivity and c is the electromagnetic wave propagation velocity in free space (=0.3 m/ns) [25,31,59]. In addition, the depth to a known reflector method [34] can be used to calculate the GPR reflection wave velocity (v rw ) of a monostatic antenna using Equation (2): where t is the TWTT of the reflected GPR wave from a reflection boundary and D is the depth to the boundary (ASTM D6432-11). ε r of the soil just below the ground surface can be calculated using the GPR direct groundwave method, using Equation (3) [32]: where t GW and t AW are the direct groundwave and airwave arrival times, respectively, from the Tx to the Rx antenna, and x is the antenna separation. ε r can be derived from the measured volumetric soil moisture content, θ v , using an empirical model. Topp et al. [60] suggested the first model as given in Equation (4)

Study Area
The experimental site included a silage-cornfield and a grass field at the Pynn's Brook Research Station (PBRS), managed by the Department of Fisheries and Land Resources of the Government of Newfoundland and Labrador, located in Pasadena (49.073 N, 57.561 W), Newfoundland and Labrador (NL), Canada. The area is characterized by a 2-5% slope, and the depth to the bedrock is >1 m from the surface [62]. Details of the observed shallow soil profile are given in Table 1. The topsoil (ε r1 , d 1 ) is an organic soil layer with gravel. Immediately below the top layer is the Ap horizon (ε r2 , d 2 ), classified as loamy sand (sand = 82.0 ± 3.4%; silt = 11.6 ± 2.4%; clay = 6.4 ± 1.2%) [63]. The average bulk density and porosity of the loamy sand layer (n = 28) are 1.31 g cm −3 (±0.07) and 51% (±0.03), respectively [63]. A well-sorted sandy soil layer (ε r3 , d 3 ) was observed between the depths of 0.35 and 3.47 m, which is likely to continue further down and act as a shallow aquifer. The average annual precipitation at the site is 1113 mm per year, with 410 mm falling as snow, and the annual mean temperature is 4 • C. Both parameters are as recorded at the nearest weather station, Deer Lake, NL, for last 30 years [64].

Data Collection and Basic Proccessing
The following materials and instruments were used for data acquisition, processing, and interpretation in this study:  (Figure 3a). An additional temperature sensor was installed with the soil moisture probe at a depth of 0.2 m. GPR surveys were carried out between the soil moisture probes and the borehole (Figure 3b). The soil moisture probes were oriented perpendicular to the GPR survey direction.
Background GPR surveys were carried out: (i) before construction of the borehole, (ii) after the construction of the borehole, but before installation of the water level sensor, and (iii) after installation of the water level sensor. Sixteen 250 MHz GPR CO surveys (42 m  Three basic GPR data processing steps were applied using the EKKO Project V3 R1 Software (Sensors and Software Inc., Mississauga, ON, Canada), as listed below.

•
Edit the first break (time-zero correction) • Apply dewow and SEC2 (Spreading and Exponential Compensation) gain • Background subtraction-applied to the full length of the trace The dataset did not require extra processing such as bandpass filters, and the processing was intentionally simplified to be applicable for every GPR profile in the dataset due to the low noise level observed. After completing this basic processing, GPR files were exported to the IcePicker V3 R7 software (Sensors and Software Inc., Mississauga, ON, Canada) for automatic time picking.

Data Analysis
CMP surveys were carried out using the borehole as the common midpoint at every field campaign to calculate the average v rw from the surface down to the top of the capillary fringe. GPR CO survey traces nearest to the borehole were only considered for the D CF estimation. Accordingly, the mean TWTT to the capillary fringe reflection (mean t CF ) was obtained from twelve GPR traces. The subjective error was high when picking the leading edge of the wavelet [39]. Therefore, TWTT related to the absolute maximum amplitude of the airwave (t air ) and the reflection event (t reflect ) were picked. This procedure was similar to the direct groundwave analysis by Grote et al. [65]. The actual t CF of every single trace was determined using Equation (5).
D CF was then calculated from the mean t CF and the average v rw using Equation (2). D CF of eight GPR surveys in 2017 were plotted in a linear regression model against the measured WTD (WTD m ) using the water level sensor at the same time as the GPR surveys. The regression equation obtained from 2017 data was used to predict WTD (WTD p ) for eight survey days in 2018. The WTD p and WTD m for 2018 were compared using a 1:1 plot and root mean square error (RMSE). In a second step, WTD m and D CF for all 16-survey days were plotted in a linear regression plot to estimate the average capillary height. The slope and the intercept of the regression line were compared statistically with those of the 1:1 line. In the same manner, the prediction line was statistically compared with the 1:1 line.

GPR Survey Outputs
GPR profiles from 250 MHz provided high-resolution radar images; therefore, multiple reflections near the water table can be clearly observed ( Figure 4). The elevations have not been corrected; therefore, the reflections in Figure 4 can be observed as undulating boundaries. The GPR velocity is independent of the frequency (for frequencies above 100 MHz to 1500 MHz) and dependent only on the relative permittivity and the magnetic permeability [67]. Thus, the same velocity derived from 250 MHz could be used to analyze 100 MHz data on the same day; 100 MHz gives low-resolution radar images with a relatively flat and clear boundary for the water table zone, and with fewer multiple reflections (Figure 5a). The results from this study indicate that 250 MHz is suitable to examine shallow D CF , whereas 100 MHz is more suitable to examine deeper WTD.

Estimation of D CF
The D CF derived using Equation (2), and CMP derived radar velocities for all 16 GPR survey days, are given in Table 2. Two examples of semblance analyses for CMP velocity estimation are illustrated in Figure 6 . These data imply that the 2017 growing season was relatively dry, which was also confirmed by weather data collected onsite (Appendix A).    Table 2. Ground penetrating radar (GPR) estimated depth to the top of the capillary fringe (D CF ) which is derived from the mean two-way travel time (TWTT) to the capillary fringe reflection (t CF ) and common mid point (CMP) derived radar wave velocity (v rw ) using Equation (2)

Site-Specific Relationship for WTD m vs. D CF
The WTD m at the same time as the GPR survey and corresponding D CF obtained using GPR data are plotted in Figure 7a. A strong linear regression with an R 2 of 0.98 was found (Equation (6)) between WTD m and D CF for GPR surveys in 2017: The WTD p for 2018 based on the GPR measured D CF and the regression model (Equation (6)) was plotted against the WTD m in a 1:1 plot (Figure 7b). The slope of the prediction line (1.6) and that of the 1:1 line (1.0) are significantly different at α = 0.05 (p-value = 0.004, df = 12, t critical = 2.179 < t calculated = 3.536) (Appendix B, Table A1). The error of WTD prediction was high during the wet survey days and overestimated from the 1:1 line (Figure 7b). This behaviour could be due to the fact that the capillary fringe would not fluctuate uniformly with the WTD fluctuation. As Bentley and Trenholm [55] stated, the capillary height is greater when the WTD is increasing (during discharging) and lower when it is decreasing (during recharging). However, the present study determines the capillary height based on the WTD m through a regression equation (Equation (6)). The regression equation is a generalized form of all measured data throughout a growing season. Therefore, the regression model might not be suitable when there is a sudden decrease in WTD, like during heavy or long duration rain events. The rain event at the end of the growing season of 2018 was unexpectedly heavy (Appendix A); consequently, it produced a maximum error in the WTD prediction on 31 October 2018 ( Figure 8). However, it is worth noting that the RMSE of 0.194 m between WTD m and WTD p can be considered acceptable for the scale of application in most agricultural practices. In general, the spatial variation of D CF cannot be measured directly under heterogeneous field conditions [68]. The proposed method provides a noninvasive approach to estimate D CF , which is more beneficial in agricultural fields during the growing season since the capillary zone provides the best soil moisture conditions for plants. The advantage of the proposed method is that both WTD and D CF can be estimated in real time. The results would have been improved if a broader range of measured data were available under different soil moisture and variable water table conditions. As seen in Figure 9, WTD m and D CF for all survey days have a linear relationship (WTD m = 1.0123 D CF + 0.741) with an R 2 of 0.91. The slopes of the regression line (1.01) and the slope of the 1:1 line (1.00) are not significantly different at α = 0.05 (p-value = 0.89, df = 28, t critical = 2.048 > t calculated = 0.146) (Appendix B, Table A2). Therefore, the intercept of the regression line (0.74 m, n = 16) can be considered as the average capillary height within the growing seasons in 2017 and 2018. As a result, an average capillary height of 0.74 m was suggested for the particular site throughout the growing season. The average value agrees with the value of~0.70 m capillary height for similar soil conditions described by Liu et al. [69]. However, in contrast, Saintenoy and Hopmans argued that D CF cannot be estimated directly because the maximum reflection is not related to the capillary fringe, but rather, to the inflection point of the transition curve [70]. It is worth noting that in the present study, we believe that transition zone reflections have no effect due to their relatively larger wavelength (i.e., low frequency) when compared to the transition zone thickness, as described by Bano [9].

Calculation of ε r
There are three soil layers in the experimental plot, as described in Section 2.2. The first layer occurs at 0-0.05 m soil depth (d 1 = 0.05 m, ε r1 ), the second layer between 0.05-0.35 m (d 2 = 0.30 m, ε r2 ), and the third from 0.35 m to the top of the capillary fringe (d 3 , ε r3 ). The average ε r from the surface to the top of the capillary fringe (ε r−avg ) can be obtained from Equation (1), because v rw is known from the CMP surveys. It is worth noting that CMP analysis might be suitable for layer velocity estimation as well; however, in the present study, the resolution was not sufficient to extract the layer velocities separately, due to low soil layer thickness. Instead, ε r1 can be calculated using the GPR direct groundwave arrival time (Equation (3)), while ε r2 can be calculated using soil moisture data (Equation (4)). Once ε r−avg , ε r1 , and ε r2 are known, ε r3 , which has the key contribution to the ε r−avg , can be calculated.
The calculated ε r3 was compared with the literature and onsite weather data to evaluate the results of this study. The average WTD m in all GPR survey days (n = 16) was 2.55 m. When the average capillary height was approximated as 0.70 m (according to [69] and the present study), the average d 3 was calculated as 1.50 m (2.55 − 0.70 − (0.05 + 0.30)). Figure 10 illustrates the overlapping of these three soil layers with the sampling volume geometry of the GPR CO survey, as suggested by Illawathure et al. [71]. Twelve GPR traces collected near the borehole and the soil moisture probes (i.e., the same as in the calculation of mean t CF ) were considered for the illustration and calculations below.  Figure 2b). Note that GPR wave paths are assumed to be straight.
The soil moisture data logging interval was 60 min. Therefore, each daily mean soil moisture data point had 24 replicated measurements. Each soil moisture probe had a cylindrical sampling volume (radius~0.05 m, volume = 0.715 L) that covered 0.05 m soil heights, both above and below the respective probe (Figure 3b) [72,73]. Daily mean soil moisture values at three depths were converted to daily mean ε r values using Equation (4); an average ε r for the soil layer between 0.05-0.35 m (ε r2 ) was obtained from those ε r values.
If the weighted average of ε r1 , ε r2 , and ε r3 was considered as ε r−avg , the percentage sample area of each layer should be calculated with respect to the total GPR sample area ( Figure 10 and Table 3). If w is the weight and x the data number, then the weighted average (x) equals: e.g., estimation of the ε r3 on 6/23/2017 (first data row in Table 4)  Table 3. Percentage sample area of each soil layer out of the total ground penetrating radar (GPR) sample area related to twelve GPR traces collected from the top of the soil moisture probes on the vertical plane of the GPR survey (refer to Figure 10).

Soil Layer
Polygon (Refer to Figure 9

Correlation Analysis
The results of the Pearson's correlation test are shown in Figure 11. The dataset and the detailed correlation analysis are given in Appendix C (Table A3 and Figure A3, respectively). The present results reveal that most of the variables tested were not correlated with WTD (also see Figure A1). However, the WTD had strong negative correlations with ε r−avg (r = −0.9041, p = 0.0001) and ε r3 (r = −0.9019, p = 0.0001). Dry sand had lower ε r than wet sand [74,75]. According to Equations (1) and (2), ε r and WTD are inversely proportional, based on the fact that water table recharging makes unsaturated soil wetter, and water table discharging makes unsaturated soil drier [73]. It was also clear that ε r−avg and ε r3 were strongly positively correlated (r = 0.9961, p = 0.0000), as ε r3 was dominant in the ε r−avg calculation ( Figure 10 and Table 3). WTD had negative (but moderate) correlations with P-E 11 through P-E 20, though the p-value of WTD vs. P-E 12 was slightly above the α-level (p = 0.0534). The highest correlation among P-E vs. WTD was for P-E 16 (r = −0.6285, p= 0.0286). This suggested that the WTD does not quickly respond to variables like daily P, daily E, or daily P-E, and that there is a long time-lag between WTD response and P-E timing ( Figure A2). This could be due to high runoff, relatively thicker aquifer, etc. [76,77]. Since it is a sandy aquifer, low infiltration was not a factor in the long responsive time of the WTD to the P-E [78,79].
GW-EC or GW-Temp were not correlated with any of other variables tested. The Temp20 had a strong negative correlation with P-E 8, P-E 9, P-E 10, and P-E 11, with the maximum at P-E 10 (r = −0.8478, p = 0.0005). Further, Temp20 had moderate correlations with SM20 and SM30 as well.
To support these observations, the collected TDR data were examined. First, the SM10, SM20, and SM30 variables had strong positive correlations with each other (r > 0.95, p = 0.00), since all three values mostly represented the same soil layer. Next, SM10 had no correlation with P, E, or P-E, while both SM20 and SM30 had no correlation with P, E, P-E, P-E 1, or P-E 2. The highest correlations between SM and cumulative P-E were as follows: between SM10 and P-E 6 (r = 0.7121, p = 0.0094), between SM20 and P-E 6 (r = 0.7486, p = 0.0051), and between SM30 and P-E 10 (r = 0.7920, p = 0.0021). From SM10 to SM30, the time lag of P-E increased, as did the strength of the correlation. These results imply that cumulative P-E affected the soil moisture as a function of depth and time. Cumulative P-E also significantly affected the soil temperature at 20 cm, but not the temperature of the groundwater. The same idea was reflected from the correlation analysis between ε r and other variables as well [80,81]. ε r1 had no correlation with any of the variables tested. It is worth noting that the depth of this first soil layer was only 5 cm. It is obvious that ε r2 had a strong positive correlation with SM10, SM20, and SM30 (r > 0.98, p = 0.00), since ε r2 was calculated from TDR data. Also, ε r2 had a negative but moderate correlation with Temp20 (r = −0.6551, p = 0.0208). There are positive moderate correlations between ε r2 and from P-E 3 to P-E 20. ε r3 also had a negative but moderate correlation with Temp20 (r = −0.7226, p = 0.0079), and a moderate positive correlation with P-E 11 (r = 0.6516, p = 0.0217) and beyond, except P-E 12 and P-E 15. This was because part of the cumulative P-E infiltrates through the soil, increasing soil moisture and lowering soil temperature [81,82]. The soil moisture governs the dielectric properties of the soil; therefore, the correlation analysis suggested that the behavior of the estimated and calculated data was supported by the measured and onsite weather data.

Summary
In GPR, the Rx records different amplitudes of the receiving signals with respect to wave travel time. The GPR interpreter observes the reflected or direct radar waves at Rx in a radargram and obtains relevant TWTTs. Without knowing the v rw , it is impossible to derive D CF (e.g., using Equation (2)). Under these circumstances, there are two challenges to estimating the D CF in a GPR radargram: First, picking the TWTT of the capillary fringe reflection correctly; and second, knowing or properly assuming the average v rw from the surface down to the top of the capillary fringe reflection.
The procedures described in the literature were combined to overcome the challenges of picking the correct TWTT. The standard error of the mean t CF was > 1 except in 4 cases in 2017: 1.44 ns (28 July), 1.26 ns (7 November), 1.11 ns (15 September), and 1.02 ns (3 October). The standard deviation (SD) of the mean t CF was above 3 ns for those four survey days. The error associated with the surveys was mostly due to a high signal-to-noise ratio in the data acquisition. The lowest SDs of the mean t CF were 1.70 ns and 1.35 ns for 23 June 2017 and 29 June 2018, respectively. The inconsistency of the GPR reflection amplitude under wet conditions, as observed by Lunt et al. [21], could be a reason for the error of time picking observed in our data set.
The challenge of defining v rw during GPR data analysis is associated with determining the average ε r of the material above the capillary fringe. ε r controls the v rw and the reflection coefficients at the interfaces. The common v rw values suggested in the literature may not be accurate for heterogeneous soil profiles in highly variable field conditions under different agricultural practices. In addition, seasonal fluctuation of WTD and capillary height can significantly change the average ε r . However, v rw can be measured by using multi-offset GPR survey methods such as CMP and WARR (wide-angle reflection and refraction), even though it is time-and labor-consuming to carry out multi-offset surveys in every field campaign [5,32]. In the present study, three different ε r values were assumed for different soil layers above the capillary fringe. A weighted average method using CMP data, ground wave analysis, and TDR data based on the GPR sampling geometry was used to calculate ε r of unsaturated soil layers in order to evaluate the velocity estimations.
These types approaches might help in addressing natural barriers of inherent soil properties and the effect of their variability on agriculture. Maintaining the required soil moisture in the root zone will lead to optimum crop yields, and implementing the best possible water management strategies will reduce groundwater pollution threats due to intensive agriculture. Nevertheless, the average capillary height considered (0.70 m) based on the study of Liu et al. [69] to estimate the ε r3 is closer to the average capillary height obtained from this study (0.74 m). It should be mentioned that taking an average capillary height would be reasonable under static conditions, but this would not always be suitable if the seasonal fluctuation of WTD is significant during the study period. Acknowledgments: Authors wish to acknowledge the Department of Fisheries and Land Resources (FLR) of the Government of Newfoundland and Labrador, NL for providing the field site and other field facilities to carry out this research. Special thanks to Adrian Unc for giving constructive comments for the initial draft. Teamwork, including the data collection support by Kamalesh Sadatcharam, Ivo Arrey, Emmanuel Badewa, and Dinushika Wanniarachchi, is highly acknowledged.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix B. Summary of Statistical Analysis
Comparing the Slopes of Regression Line and Prediction Line with 1:1 Line To test whether the slopes for two independent populations are equal, the following null and alternative hypotheses were tested: Hypothesis 0 (H 0 ). β 1 = β 2 i.e., β 1 − β 2 = 0 The test statistic is ∼ T(n 1 + n 2 − 4) n = sample size; b 1 and b 2 are slopes, s y.x = standard error of predicted y for each x in the regression, s x = standard deviation. If the null hypothesis is true, then  Since t > t-critical and p-value < α, two slopes are significantly different at α = 0.05.