Temporal Monitoring of the Soil Freeze-Thaw Cycles over a Snow-Covered Surface by Using Air-Launched Ground-Penetrating Radar

Abstract

We tested an off-ground ground-penetrating radar (GPR) system at a fixed location over a bare agricultural field to monitor the soil freeze-thaw cycles over a snow-covered surface. The GPR system consisted of a monostatic horn antenna combined with a vector network analyzer, providing an ultra-wideband stepped-frequency continuous-wave radar. An antenna calibration experiment was performed to filter antenna and back scattered effects from the raw GPR data. Near the GPR setup, sensors were installed in the soil to monitor the dynamics of soil temperature and dielectric permittivity at different depths. The soil permittivity was retrieved via inversion of time domain GPR data focused on the surface reflection. Significant effects of soil dynamics were observed in the time-lapse GPR, temperature and dielectric permittivity measurements. In particular, five freeze and thaw events were clearly detectable, indicating that the GPR signals respond to the contrast between the dielectric permittivity of frozen and thawed soil. The GPR-derived permittivity was in good agreement with sensor observations. Overall, the off-ground nature of the GPR system permits non-invasive time-lapse observation of the soil freeze-thaw dynamics without disturbing the structure of the snow cover. The proposed method shows promise for the real-time mapping and monitoring of the shallow frozen layer at the field scale.

1. Introduction

The interactions between seasonal soil frost and snow cover have a substantial impact on the soil moisture balance and energy exchange processes between the soil and atmosphere, as the hydraulic and thermal properties of unfrozen and frozen soils are quite different. The influence of snow cover on the ground thermal regime depends on the accumulation, duration, timing and thickness of the seasonal snow layer, as well as melting processes and geographical location [1]. The soil frost depth is inversely correlated with snow cover depth, as snow acts as an insulating medium between colder air temperature and the soil surface [2]. The snow cover retained by vegetation on the ground surface can enhance the insulation of the soil, which may reduce periodic variation in soil temperature, slowing the freezing and thawing processes [3] and even altering the density of the snow layers [4]. An excellent review is given by Zhang et al. [5] to address the effects of seasonal freeze-thaw and snow cover processes on the ground thermal regime.

Traditionally, investigation of the freeze-thaw cycles are mainly characterized by using invasive in situ point-measurements (e.g., frost tubes, thermistors, dielectric permittivity sensors), which provide limited local scale information and disturb the natural soil system during installation. For instance, Wagner et al. [6] used thermistors to monitor the temperature variation along a cylindrical sample placed in an environmental cell. Frost tubes (plastic fluorescein dye tubes) are usually filled with a solution of methylene blue that undergoes a color change (blue to clear) as a result of freezing. Frost tube readings are taken manually and therefore significantly disturb the subsurface structure during installation, which may influence local soil moisture and local heat exchange. Time domain reflectometry (TDR) has also been used to monitor the dynamics of seasonally-frozen soils and meltwater infiltration due to its ability to measure the soil dielectric permittivity [7,8]. Nevertheless, the TDR method requires the installation of several vertically- and horizontally-positioned probes at fixed soil depths, thereby restricting its applicability in larger spatial investigations.

In the last few decades, satellite remote sensing techniques have been used to measure the surface soil freeze-thaw condition at both regional and global scales [9,10]. Frozen surface and soil conditions have been monitored by both active and passive microwave sensors [11,12,13,14,15,16], as well as using optical sensors [17]. The large contrast between the dielectric properties of ice and water makes it suitable to monitor the freeze-thaw cycle of soil using microwave measurements. Current approaches employing satellite microwave remote sensing to, e.g., to classify landscape freeze-thaw states, involve temporal change detection schemes applied to time series of backscatter (e.g., from ESA’s ERS-1 satellite [11]) or applied to brightness temperatures (e.g., from the special sensor microwave imager (SSM/I) [18] or from the advanced microwave scanning radiometer - earth observing system (AMSR-E) [19]). The retrieval of information on snow and soil information from remote active microwave data (backscatter) is less developed in comparison to the corresponding retrievals using passive microwave data (brightness temperatures) [20]. Additionally, most remote sensing products are either coarse in spatial (especially passive sensors) or in temporal resolution [21,22]. In this respect, ground-penetrating radar (GPR) has been used as a non-invasive method to monitor the spatial and temporal dynamics of the shallow freeze and thaw processes at field scales [23,24,25].

Recent studies have illustrated that GPR can be useful to characterize a wide range of glacial and frozen conditions [26,27,28], such as monitoring the spatial-temporal distribution of frozen ground, mapping thaw fronts or assessing freeze-thaw dynamics, the thickness of the active layer above the permafrost and the depth to the permafrost table [29,30]. Furthermore, rapid deployment of GPR along series of parallel transects can be used to estimate frozen depth over a broad area. Hinkel et al. [31] reported that through the use of the GPR system, the depth of thaw within the active layer could be mapped with a precision on the order of 10 cm vertically, while Forte et al. [32] proposed an empirical equation for common offset GPR measurements to estimate the density of frozen media. A high frequency GPR direct ground wave has also been tested to monitor the dynamics of the seasonal development of the thin soil frozen layer [33]. Recently, Butnor et al. [34] measured soil frost depth in a forest ecosystem using GPR and reported that deep snow and shallow frost can make frost detection with GPR difficult, due to signal attenuation. Since commercially-available GPR systems are mostly ground coupled, the structure of snow cover or frozen soil can be significantly disturbed by dragging the GPR antenna over the surface. As a result, precise information on shallow freeze-thaw dynamics is difficult to monitor via time-lapse on-ground GPR measurements, without disturbing the snow cover substantially.

In this research, we tested an off-ground GPR system to monitor the seasonal freeze-thaw dynamics of soil over a snow-covered surface. A horn antenna was used in the frequency range of 200–800 MHz and fixed at a height of about 110 cm above the ground surface in the field. Within the footprint of the horn antenna, the surface of the bare soil was exposed to the natural environmental conditions, such as precipitation, snowfall and evaporation. Near the GPR setup, sensors were installed in snow and at different depths in a trench to monitor the vertical and lateral dynamics of soil temperature and dielectric permittivity. Time-lapse GPR, soil temperature and dielectric permittivity measurements were recorded for nine days, with five freeze-thaw cycles occurring during that period. GPR-derived permittivity was compared to the permittivity obtained by the sensors.

2. Material and Methods

2.1. GPR System

The GPR system consisted of an ultra-wideband stepped-frequency continuous-wave (SFCW) radar combined with a monostatic linearly-polarized double-ridged broadband horn antenna (BBHA 9120-F, Schwarzbeck Mess-Elektronik, Schnau, Germany). A vector network analyzer (VNA; ZVRE, Rohde & Schwarz, Munich, Germany) was used to set up the SFCW radar. The dimensions of the BBHA 9120-F are $68\times 95$ cm${}^2$ aperture area and 96 cm in height, with an isotropic gain range from 6–18 dBi and a nominal frequency range of 0.2–2 GHz. The relatively small 3-dB beamwidth of the antenna at 45${}^{\circ }$ in both the E- and H-plane and a frequency $f=1$ GHz make it suitable for using off ground. A high quality N-type 50-Ohm impedance coaxial cable of 2.5 m in length (Sucoflex 104PEA, Huber+Suhner AG, Herisau, Switzerland) was used to connect the antenna to the reflection port of the VNA. An Open-Short-Match reference calibration kit was used to calibrate the VNA at the connection between the antenna feed point and the cable. An interface between the VNA and a field laptop was developed to automatically communicate with the GPR system and to save the time-lapse GPR measurements on a field laptop.

2.2. GPR Forward Modeling

The electromagnetic model used to describe the GPR wave propagation in the antenna-air-soil system is considered as a multilayered medium, details of which can be found in Lambot et al. [35]. The model is applied to an off-ground horn antenna connected to the SFCW radar. In the model, the ratio between the backscattered electromagnetic field $b\left(\text{ω}\right)$ and the incident electromagnetic field $a\left(\text{ω}\right)$ is expressed in the frequency-dependent complex ratio $S_{11}\left(\text{ω}\right)$, where ω is the angular frequency ($\text{ω}=2\pi f$). In the far-field assumption, the parameters of the antenna are assumed independent of the scatterer, i.e., the plane wave field approximation is distributed over the antenna aperture. The complex, frequency-dependent global transmission and reflection coefficients can be used to model the GPR system precisely via Maxwell’s equations [35]:

$$S_{11}\left(\text{ω}\right)=\frac{b\left(\text{ω}\right)}{a\left(\text{ω}\right)}=H_i\left(\text{ω}\right)+\frac{H_t\left(\text{ω}\right)H_r\left(\text{ω}\right)G_{xx}^{\uparrow }\left(\text{ω}\right)}{1-H_f\left(\text{ω}\right)G_{xx}^{\uparrow }\left(\text{ω}\right)}$$

(1)

where $H_f\left(\text{ω}\right)$, $H_t\left(\text{ω}\right)$, $H_r\left(\text{ω}\right)$ and $H_i\left(\text{ω}\right)$ are the characteristics of the antenna for the feedback loss, transmitting, receiving and the complex return loss transfer functions of the antenna, respectively. The $G_{xx}^{\uparrow }\left(\text{ω}\right)$ term is the Green function impulse reflection of the air-soil system modeled as a multilayered medium. $H_i\left(\text{ω}\right)$ represents the variation of impedance between the antenna aperture, feed point of the antenna and multiple wave reflections occurring within the antenna, which are independent of the backscattered electromagnetic field $G_{xx}^{\uparrow }\left(\text{ω}\right)$. The $H_t\left(\text{ω}\right)$ and $H_r\left(\text{ω}\right)$ transfer functions correspond to the phase delay and antenna gain between the source-receiver virtual point and the measurement point. To reduce the number of transfer functions, Equation (1) can be defined as $H\left(\text{ω}\right)=$$H_t\left(\text{ω}\right)$$H_r\left(\text{ω}\right)$. The characteristics of the antenna transfer functions ($H_f\left(\text{ω}\right)$, $H_i\left(\text{ω}\right)$ and $H\left(\text{ω}\right)$) can be estimated by solving a system of equations as Equation (1) for these three unknown functions, by performing $S_{11}$ measurements above a perfect electrical conductor (e.g., a copper sheet), for which the Green functions $G_{xx}^{\uparrow }$ can be simulated [36,37,38]. The antenna effects can be removed from the raw radar measurements, and once the antenna transfer functions are known together with Green’s function, the response from the air-soil medium can be derived using Equation (1).

The earth medium properties govern the electromagnetic wave propagation and include the electric conductivity σ (S·m${}^{-1}$), electric permittivity ε (F·m${}^{-1}$) and magnetic permeability μ (H·m${}^{-1}$). In this study, μ is assumed to be equal to the permeability of free space, ${\text{μ}}_0=4\pi \times {10}^{-7}$ H·m${}^{-1}$, which is valid for non-magnetic soil materials, such as natural soil. The subsurface is modeled as a 3D multilayered medium consisting of N horizontal layers separated by $N-1$ interfaces. The characteristics of the n-th layer is considered homogeneous and expressed as ${\text{ε}}_n$, ${\sigma }_n$ with a thickness $h_n$. The electromagnetic waves propagating in a multilayered media are defined by Maxwell’s equations. The antenna is monostatic and modeled as a point source and receiver, and hence, the radiation pattern is emulated by that of a point dipole. The analytic formulation for the zero-offset Green’s function in the spectral domain (2-D spatial Fourier domain) can be expressed as [35]:

$$\tilde{G}\left(k_{\rho }\right)=\frac{1}{8\pi }\left(\frac{{\Gamma }_1{R}_1^{TM}}{{\eta }_1}-\frac{{\zeta }_1{R}_1^{TE}}{{\Gamma }_1}\right)exp\left(-2{\Gamma }_1{h}_1\right)$$

(2)

In the equation, $R^{TE}$ and $R^{TM}$, corresponds to the transverse electric and the transverse magnetic global reflection coefficients accounting for all reflections in the multilayered medium, respectively [39]. Γ is the vertical wavenumber and determined as $\Gamma =\sqrt{k_{\rho }^2-k^2}$, whilst $k^2={\text{ω}}^2\text{μ}(\text{ε}-\frac{j\sigma }{\text{ω}})$ with ω being the angular frequency. For the free-space Layer 1, $k_1^2={\left(\frac{\text{ω}}{c}\right)}^2$, where c is the free-space wave velocity.

The 2-D Fourier inverse transformation is applied to transform Equation (2) from the spectral domain to the spatial domain:

$$G={\int }_0^{+\infty }\widetilde{G_{xx}^{\uparrow }}\left(k_{\rho }\right)k_{\rho }d{k}_{\rho }$$

(3)

which reduces to a single integral in view of the invariance of the electromagnetic properties along the x and y coordinates. Lambot et al. [40] proposed an optimal procedure to properly evaluate the integral in Equation (3), which has singularities. Singularities (branch points and poles) are avoided by deforming the integration path in the complex $k_{\rho }$ plane. Furthermore, the integration becomes faster as the oscillations are minimized by defining an optimal path. To define $k_{\rho }$ as the complex number ($x+jy$), the following relationship was used for the constant phase integration path:

$$y\left(x\right)=\frac{x}{\sqrt{{\left(\frac{xc}{\text{ω}}\right)}^2+1}}$$

(4)

where c is the free-space electromagnetic wave velocity.

2.3. Time Domain Inversion Method

The surface dielectric permittivity of the soil is estimated via the time domain inversion of Green’s function, with a focus on a time window containing the ground reflected first arrival [41]. The inverse Fourier transform is used to transform the measured and modeled frequency domain Green’s functions into the time domain. The inversion procedure is defined as a classical, non-linear least-squares problem, and the objective function to be minimized is formulated as follows:

$$\varphi \left(\mathbf{p}\right)={\left({\mathbf{g}}_{\mathbf{xx}}^{\uparrow *}-{\mathbf{g}}_{\mathbf{xx}}^{\uparrow }\right)}^T\left({\mathbf{g}}_{\mathbf{xx}}^{\uparrow *}-{\mathbf{g}}_{\mathbf{xx}}^{\uparrow }\right)$$

(5)

where:

$${\mathbf{g}}_{\mathbf{xx}}^{\uparrow *}={\left.g_{xx}^{\uparrow *}\left(t\right)\right|}_{t_{min}}^{t_{max}}\text{and}{\mathbf{g}}_{\mathbf{xx}}^{\uparrow }={\left.g_{xx}^{\uparrow }\left(t\right)\right|}_{t_{min}}^{t_{max}}$$

(6)

are the vectors representing the observed and simulated time domain windowed Green’s functions, respectively (see Figure 1); $\mathbf{p}=[{\text{ε}}_r,h]$ is the parameter vector to be inversely estimated; ${\text{ε}}_r$ is the soil surface relative dielectric permittivity and h the distance between the soil surface and the antenna phase center. In the far-field region of the antenna, the spherical divergence of the radiated field is initiated from the virtual source point, which is the antenna source point. Lambot et al. [35] proposed a procedure to estimate the phase center of an off-ground monostatic antenna. Furthermore, Jadoon et al. [43] reported that the frequency-dependent phase center position of the antenna is inherently accounted for in the antenna transfer functions of the GPR forward model. Although the magnetic permeability, electrical conductivity and soil layering can be taken into consideration in the inversion process, in this study, we assume them to have a negligible effect on the estimation of ${\text{ε}}_r$, which is a commonly-adopted assumption [41].

Figure 1. Measured Green’s function in the time domain, where $t_{min}$ and $t_{max}$ define the time window focused on the surface reflection in which inversion is performed. Time $t_i$ corresponds approximately to the surface interface in the space domain, and the maximum amplitude of Green’s function is defined by the peak-to-peak $PtP$.

Figure 1. Measured Green’s function in the time domain, where $t_{min}$ and $t_{max}$ define the time window focused on the surface reflection in which inversion is performed. Time $t_i$ corresponds approximately to the surface interface in the space domain, and the maximum amplitude of Green’s function is defined by the peak-to-peak $PtP$.

The misfit between the measured and modeled Green’s function (5) is minimized by the Levenberg-Marquardt optimization algorithm. The surface reflection arrival time $t_i$ (see Figure 1) is detected automatically and used to drive the initial guess for the h. A fast evaluation of Green’s function was used [40], and model inversion is achieved in less than one second, which permits real-time processing applications.

3. Test Site and Experimental Setup

Measurements were performed at the Terrestrial Environmental Observatories (TERENO) test site in Selhausen, which is located in the Rur River catchment, North Rhine-Westphalia, Germany. The underlying fluvial deposits from the Rhine/Meuse River and the Rur River system are Quaternary sediments, covered by eolian sediments (up to a depth of 1 m) from the Pleistocene and Holocene. The test site has a slope of 4${}^{\circ }$, and the altitude ranges between 100 and 110 m above sea level. A mean annual precipitation of 690 mm and an average air temperature of 9.8 ${}^{\circ }$C are observed at the test site. According to the USDA textural classification, the major soil type is silty loam. Further information along with a detailed description of the test site is given by Ali et al. [44].

Figure 2. Decagon EC-TM capacitance probes were installed near the footprint of the ground-penetrating radar (GPR) antenna in a trench to monitor the vertical and lateral dynamics of soil temperature and dielectric permittivity (a), and the fully-automated off-ground GPR setup with an antenna fixed at 110 cm above the ground (b) (after [42]).

Figure 2. Decagon EC-TM capacitance probes were installed near the footprint of the ground-penetrating radar (GPR) antenna in a trench to monitor the vertical and lateral dynamics of soil temperature and dielectric permittivity (a), and the fully-automated off-ground GPR setup with an antenna fixed at 110 cm above the ground (b) (after [42]).

Figure 2 shows a small trench with Decagon EC -TM probes installed at a 1.5-m distance from the footprint of the GPR antenna and the fully-automated GPR experimental setup. Measurements were carried out on a bare agriculture field. Soil temperature and permittivity were recorded by EC-TM sensors (Decagon Devices Inc., Pullman, WA, USA), and Em50 data loggers (Decagon Devices Inc., Pullman, WA, USA) were used to record and store the measurements. Air temperature was recorded at a weather station located at a distance of 800 m from the setup. Soil and air temperature measurements were carried out with a constant time step of 15 min for nine days. The lateral and vertical dynamics of the soil temperature and dielectric permittivity were recorded at six depths (2, 3, 4, 5, 7, 8 cm). Two additional EC-TM probes were installed in the snow at 2 cm above the ground to monitor the dynamics of snow temperature and permittivity. The stepped-frequency continuous-wave radar system was set up using a vector network analyzer, with a monostatic antenna, thereby producing an air-launched radar. First, the antenna calibration experiment was performed in the laboratory, and later, the GPR system was set up in the field. The double-ridged horn antenna was fixed at approximately 110 cm above the ground, and the surface of the bare soil was fully exposed to natural processes, i.e., precipitation, evaporation and snowfall (see Figure 2b). Time-lapse GPR, temperature and permittivity measurements were performed during days 330.5–339.5 of 2010 to monitor the soil dynamics. Data acquisition was computer controlled and entirely automated over the nine-day measurement period [42].

4. Results and Discussion

Figure 3 shows the time series of air and soil temperature recorded at different depths over the nine-day measurement period. Air temperature was mostly below $0{}^{\circ }\mathrm{C}$ and reached a minimum level of $-10{}^{\circ }\mathrm{C}$ on Day 8 of the measurement campaign. The vertical dynamics of the soil temperature recorded at six depths indicate five freeze-thaw events during the measurement period, which are identified by soil temperatures below $0{}^{\circ }\mathrm{C}$, as indicated by arrows in Figure 3b. The snow cover fluctuated over the period from 3–10 cm due to precipitation, ablation, sublimation, but also wind drift. Ambient temperature fluctuates around $0{}^{\circ }\mathrm{C}$. Seasonal snow cover has a cooling effect on the soil surface when the snow cover is comparatively thin and with a high albedo. This phenomena can be observed in the second freeze-thaw event, where soil temperature drops rapidly to $-1.7{}^{\circ }\mathrm{C}$ in the top 3 cm (Figure 3b). In Figure 3a, the air temperature dropped below $-5{}^{\circ }\mathrm{C}$ during Days 335–338 of the year, and the snow thickness increased, causing an insulating effect of the snow cover. This insulation effect resulted in a warmer soil surface. As a consequence, soil temperature in the top 3 cm rarely drops to below $-0.5{}^{\circ }\mathrm{C}$ (Figure 3b). Irrespective of the low air temperatures over the entire course of the experiment, the soil freeze-thaw cycles were most prevalent in the top 4 cm of the soil. In the cold regions, the snowpack persists for a longer time with a thickness in meters and can undergo compaction due to frozen water layer accumulation. Complexities related to the accumulation of frozen water layers in the snow were not observed during measurements.

Figure 3. (a) Time series of air temperature and (b) soil temperature recorded at shallow depths (2, 3, 4, 5, 7, 8 cm) over a nine-day period. Arrows indicate the freeze-thaw cycles as the soil temperature dynamics is below 0 ${}^{\circ }$C.

Figure 3. (a) Time series of air temperature and (b) soil temperature recorded at shallow depths (2, 3, 4, 5, 7, 8 cm) over a nine-day period. Arrows indicate the freeze-thaw cycles as the soil temperature dynamics is below 0 ${}^{\circ }$C.

Figure 4a,b presents the raw and processed GPR signal for the measurements carried out with the antenna placed at different heights over a perfect electrical conductor (copper sheet). In total, 17 GPR measurements were performed by lifting the antenna to different heights, with the antenna aperture elevation ranging from 0.08 m–1.2 m. The frequency-dependent GPR measurements were selected from 200–800 MHz, as higher frequencies (>800 MHz) were affected by noise arising due to soil roughness scattering [45]. Furthermore, the signal-to-noise ratio in the higher frequency range is lower due to a higher return loss and lower gain [46]. For convenience, an inverse discrete fast Fourier transform was used to convert the measured frequency domain measurements to time domain signals. In the raw b-scan $S_{11}\left(t\right)$ of Figure 4a, we observe the wave reflections occurring at the antenna feed point and within the antenna between 0 and 8 ns, with reflection from the aperture taking place at the time interval between 10 and 15 ns. In the raw signal b-scan, the second order reflection at the copper sheet occurs between 15 and 20 ns due to the feedback loss of the antenna. In the processed GPR signal $g_{xx}^{\uparrow }\left(t\right)$, the noise in the form of feedback losses and the reflections occurring within the antenna and at the antenna feed point are removed quantitatively by the antenna transfer functions, and only the response of the copper sheet is visible between 5 and 12 ns (Figure 4b). It is important to note that the reflection from the copper sheet has been delayed by the noise if we compare the time of reflection in the raw and processed GPR signal. The slight ripples in the time domain signal may be caused by the inverse Fourier transform, as the signal frequencies between 200 and 800 MHz were measured.

Figure 4. Time domain representation (b-scan) of radar measurement performed with the horn antenna at different heights above a perfect electrical conductor (copper sheet). (a) Measured raw GPR signal $S_{11}\left(t\right)$ and (b) the filtered signal $g_{xx}^{\uparrow }\left(t\right)$ using the antenna model.

Figure 4. Time domain representation (b-scan) of radar measurement performed with the horn antenna at different heights above a perfect electrical conductor (copper sheet). (a) Measured raw GPR signal $S_{11}\left(t\right)$ and (b) the filtered signal $g_{xx}^{\uparrow }\left(t\right)$ using the antenna model.

Figure 5a,b shows the time series of Green’s function in the time domain and the maximum peak-to-peak ($PtP$) amplitude of the processed signal, respectively. GPR measurements were carried out with a constant time step of 15 min for nine days, resulting in nearly 860 measurements. In general, GPR propagates the electromagnetic energy from the monostatic horn antenna and records the reflected energy from the subsurface with different electrical properties. In Figure 5a, Green’s function is obtained by filtering out noise effects via the antenna transfer functions in Equation 1. When the soil is in the process of freezing, the presence of ice (instead of liquid water) in some pores reduces the dielectric permittivity of the soil and results in small reflections in the GPR signal. At the time of thawing, the soil dielectric permittivity increases due to more water in the soil system, and the amplitude of Green’s function is relatively large due to the strong reflection at the wet soil surface. The amplitude of Green’s function varies across the nine-day period due to the different reflections received during the soil freeze and thaw process. To provide a better visual interpretation of this, the amplitude of maximum $PtP$ is plotted in Figure 5b. All five of the freeze-thaw cycles, as indicated by arrows, can be clearly observed in the $PtP$ plot of Figure 5b and follow the trend of freeze-thaw cycles observed in temperature measurements detectable in Figure 3b.

Figure 5. (a) Radar signal in the time domain for a single measurement and (b) peak-to-peak ($PtP$) maximum amplitude of the time-lapse GPR data, with the effect of the five freezing and thawing events clearly observed as marked by arrows.

Figure 5. (a) Radar signal in the time domain for a single measurement and (b) peak-to-peak ($PtP$) maximum amplitude of the time-lapse GPR data, with the effect of the five freezing and thawing events clearly observed as marked by arrows.

Reflections from snow water films or wet snow were not observed in the GPR signal. In the far-field assumption, the layer thickness of surface water films or wet snow should be more than one tenth the minimal wavelength (${\lambda }_{min}$) to obtain a reflection in the GPR radargram. Considering the ${\text{ε}}_r$ for snow water films to be 80 and for wet snow as 50, the vertical resolution of the GPR signal in the frequency range of 200–800 MHz can be 0.53 cm and 0.41 cm for the surface water films and wet snow, respectively. Reflections in the GPR radargram can be observed in the case of surface water films or wet snow exceeding this calculated layer thicknesses. The electrical conductivity of the snow cover was recorded to be below 0.001 mS/cm, which cannot attenuate the GPR signal in the low frequency range.

Figure 6 shows the relative dielectric permittivity (${\text{ε}}_r$) recorded by the in situ sensors and estimated from time-lapse GPR data. The dynamics of the freeze-thaw cycle can be observed in the ${\text{ε}}_r$ measured by the in situ sensors installed at a 2–4 cm depth over the nine-day measurement period, as indicated by arrows. Below a 4 cm depth, the ${\text{ε}}_r$ variation is less pronounced and indicates that the freezing front did not reach this depth, which is in correspondence with the recorded soil temperature measurements at these deeper depths, which never drop below $0{}^{\circ }\mathrm{C}$ (Figure 3b). As can be seen, the GPR signal is highly sensitive to the variation of ${\text{ε}}_r$ in the subsurface and the contrast between the lower ${\text{ε}}_r$ of the frozen compared to higher ${\text{ε}}_r$ of the unfrozen soil. As already mentioned, the surface reflection in the measured time-lapse GPR data was used to estimate ${\text{ε}}_r$ of the soil by inversion. With respect to the simple topography of the objective function dealing with two unknown parameters ($\mathbf{p}=[{\text{ε}}_r,h]$), the initial guess for the relative dielectric permittivity can be chosen arbitrarily. Here, it was fixed to eight. To illustrate the uniqueness of the inverse solution, tests were performed with randomly-selected initial guesses of ${\text{ε}}_r$ and showed that the solution found by the algorithm was independent of the initial guess. Additionally, the GPR estimated ${\text{ε}}_r$ matches closely to the ${\text{ε}}_r$ recorded by the sensor installed at a 3-cm depth and followed a similar trend, with the effect of the freeze-thaw cycles well reproduced. Jadoon et al. [45] reported that the time domain inversion of surface reflection is sensitive to the soil moisture dynamics in the top 4 cm of soil. Furthermore, Jadoon et al. [46] tested the full-waveform inversion of off-ground GPR measurements to estimate the depth profile of soil moisture in the top 60 cm.

Figure 6. Time series of the relative permittivity (${\text{ε}}_r$) of the soil recorded at shallow depths (2, 3, 4, 5, 7, 8 cm) and estimated from time-lapse GPR data. Five freezing and thawing cycles can be clearly observed over a period of nine days, as indicated by arrows.

Figure 6. Time series of the relative permittivity (${\text{ε}}_r$) of the soil recorded at shallow depths (2, 3, 4, 5, 7, 8 cm) and estimated from time-lapse GPR data. Five freezing and thawing cycles can be clearly observed over a period of nine days, as indicated by arrows.

As shown by the experiment, the freezing front and the freeze-thaw cycles could be detected by the off-ground monostatic GPR system, without disturbing the natural snow cover or the frozen soil. Additionally, time-lapse measurements are feasible even if the (sub-)surface may be exposed to natural environmental conditions, such as precipitation in the form of snow, wind, etc. Furthermore, antenna effects can be removed from the raw radar signal due to the accuracy of the radar electromagnetic model. In comparison, commercially available GPR systems used to detect freezing zones and permafrost dynamics in the subsurface are mostly ground coupled. For time-lapse GPR measurements, the structure of the top snow cover or frozen soil, as well as environmental conditions (e.g., precipitation in terms of snow or radiation) can be significantly disturbed and, hence, affect the fidelity of measurements. Additionally, the tested off-ground GPR system was stable with respect to ambient temperature and recorded precise measurements even at fairly low ambient temperature ($<-10{}^{\circ }\mathrm{C}$). The tested time domain inversion approach is based on surface reflection and, therefore, mainly sensitive to the top few centimeters. As such, this approach can be very useful to calibrate different remote sensing products and to deliver precise information that can be used to improve remote sensing retrieval algorithms.

In the setup presented here, soil surface roughness as well as vegetation effects, were excluded by placing the GPR system over a smooth bare soil field. In general, a differentiation between rough and smooth surfaces can be made with respect to the wavelength of the signal based on the Rayleigh criterion ($h_c=\lambda /8$), where λ is the wavelength and $h_c$ is the critical height of the surface protuberances. If the surface is smooth and homogeneous (in terms of electrical properties), most of the reflected energy will be in the specular direction (coherent component), while if the surface is rough, diffuse reflections or scattering (incoherent component) can occur, leading to less energy being recorded in the specular direction. To account for rougher surfaces, Jonard et al. [47] introduced a roughness sub-model that can be incorporated into the full-waveform inversion of the off-ground GPR model. Similar surface roughness models can be included for spatial off-ground GPR measurements operating on a large scale. Furthermore, at high frequencies (800–2000 MHz) the GPR signal has a low signal-to-noise ratio due to the far-field assumption [45]. The far-field approach (placing the antenna above a critical height) has the disadvantage of less energy penetrating into the subsoil and, therefore, a shallow investigation depth. To increase the investigation depth, Lambot et al. [48] recently proposed a near-field modeling approach for 3D electromagnetic wave propagation in planar layered media. In this case, the signal-to-noise ratio at high frequencies can be improved by considering the near-field modeling approach and placing the antenna a few centimeters above the ground without disturbing the structure of the ground surface. Dynamics within the snowpack can be observed in the high frequency range using near-field modeling for GPR measurements. The off-ground GPR system can also be mounted on mobile platforms (e.g., quad or snow mobile) to map the real-time freeze-thaw cycles over a field scale.

5. Conclusions

In this study, a monostatic off-ground GPR system was tested to monitor the temporal development of soil freezing and thawing processes in a natural setting over a snow-covered surface. For comparison of the GPR-derived information and to get detailed insight into the soil freeze-thaw state, soil temperature and relative dielectric permittivity were recorded by in situ sensors installed close to the GPR setup. Over the course of nine days, several freeze-thaw cycles of the shallow subsurface could be detected by the in situ sensors in terms of soil temperature and dielectric permittivity changes. Additionally, the freezing front could be delineated by the sensor setup.

The first results of the off-ground GPR system showed that the measurements were stable at low ambient temperature over the course of the nine days. Additionally, a detectable change in the amplitude of the surface reflection of the time domain GPR signal corresponds to the freezing and thawing cycles observed by the in situ sensors. It was shown that for the GPR setup, the snow cover and changes in the snow cover height and physical properties did not substantially influence the GPR reflection, and therefore, the information of surface reflection is from the soil surface. Finally, shallow soil ${\text{ε}}_r$ was estimated from the inversion of the surface reflection of the GPR signal, which was in good agreement with ${\text{ε}}_r$ measured by in situ sensors. Overall, the results represent a promising step toward the application of off-ground GPR measurements for noninvasive monitoring of the shallow soil freeze-thaw cycles at the field scale and can be used to calibrate the footprints of different remote sensing products.