Assessment of Soil Horizons and Their Matric Potential from Ground-Penetrating Radar Signal Attributes

Abstract

Soil plays significant roles in different phases and in the continuous existence of human life. Its comprehensive knowledge, particularly as related to its physical characteristics, enhances its utilization, conservation, and management. The traditional methods of soil study are characterized with some pitfalls such as much time needed to perform such assessments. There are also issues of invasiveness that affect the soil structures and discrete sampling that may not reflect true spatial attributes in the outcome of such techniques. These problems are largely due to the concealing nature of soil layers that made its thorough evaluation difficult. In this study, an alternative geophysical approach has been adopted. The technique is the ground-penetrating method (GPR) that utilizes electromagnetic pulse energy via its equipment’s sensors, which can allow the investigation of soil properties, even in its concealing state. This study aimed at qualitatively evaluating the soil horizons and the matric potentials using the GPR signal attributes within the unsaturated zone with a view of having insight into the test field’s characterization. Field data measurements were obtained using MALA ProEX GPR equipment with its accessories manufactured by MALA Geosciences, Stockholm, Sweden. Evaluation of the processed field data results and computed attributes show soil characteristics variations with depth that was interpreted as the layers. This can be seen from the GPR data presentation as an image representing the subsurface of the zones of propagation of the pulse energy. Spectral analysis of the GPR signals allows for the delineation of two zones of contrasting features, which were tagged as high and low matric potentials. Although the conventional direct measurement of the matric potential was not made at the time of the study to complement and confirm the veracity of the approach, the results indicate the possibility of the approach towards a quick and in situ technique of soil investigations. Such evaluation may be valuable input in precision agriculture where accurate data are sought for implementation.

1. Introduction

One of the ways to understand the nature of materials with which the subsurface is composed is a thorough examination of the inherent physical properties, and then its characteristics can be elucidated. Soil, being an essential component for the sustenance of human existence, ought to be constantly assessed for its better understanding and utilization. It is part of the Earth’s crust that starts from where the atmosphere meets the top surface and terminates at the bedrock where un-weathered material is encountered. Soil plays significant roles to both plants and animals by supplying nutrients and the base for which plants grow, as well as forming the base on which animals graze, thrive, and co-exist. Human being existence is also sustained by soil as a critical component apart from air and water. Following the significance of soil to humans, plants, and animals, its continuous investigation and evaluation cannot be over-emphasized. This is largely due to the influence of human activities since the beginning of modern time, which has significantly caused its depletion. Moreover, the world population is growing at a geometric rate, and the human need for food and shelter is also increasing. Therefore, if human, plant, and animal existence is to be sustained and the ecosystem preserved, constant soil study is necessary.

Detailed comprehension of the characteristics of the soil starts with the delineation of its emplacement and geometry in terms of relationships amongst its subsurface layers. Several physical, chemical, and biological properties of soil may reveal its nature. The present study focuses on the evaluation of soil layers and its matric potential. The choice of these inherent properties in the study is due to the possibility of inferring other physical properties such as porosity, texture, and water contents and the direction of flow from its findings. Several methods, such as the use of tensiometers, electrical resistance blocks, thermocouple psychrometers, and pressure membrane apparatus, are some of the conventional approaches to measuring soil matric potentials [1]. Furthermore, other soil assessment techniques that may be utilized include morphometric techniques that are found in the literature. Morphometric methods use scientific approaches as tools for measuring, mapping, and quantifying soil horizon properties [2]. The methods can be used to represent soil profiles with depth functions. Soil structure and porosity have been determined by the use of multi-stripe laser triangulation in the test conducted by [3]. In the review of [4], they explored the relationship between matric potential and soil water contents. They showed the convectional direct and indirect approaches of matric potential measurement, as well as their limitations, and pointed out the need for a combination of soil water contents and matric potential for soil water status determination. However, these approaches, apart from being labor-intensive and expensive, are point-based assessments, and as such, cannot be effectively used in spatial characterization of soil properties. These setbacks of the traditional methods of assessment are what informed the idea of the choice of another survey (geophysical) method of assessment that may circumvent the pitfalls. The selection of the geophysical approach for the study was due to the non-invasiveness of its technique of operations, fast nature, and less cumbersome operation with high-resolution capability. Many researchers have tested the possibility of deployment of geophysical tools in subsoil layers’ evaluations [5,6]. Galagedara et al. [7] reported the significance of GPR ground wave velocity analysis in the evaluation of soil water content. The study also compared the outcome of water content evaluation using both WARR and constant offset techniques. In the same vein, the integration of hydrogeophysical inversion and time domain reflectometry (TDR) approaches for shallow subsurface hydraulic properties characterization was documented in the study by [8]. Dumont et al. [9] attempted characterization of a landfill to estimate the spatial and vertical extents of the converted site. They tested the applicability of different methods ranging from EM (electromagnetic), to magnetic, to electrical resistivity tomography (ERT) methods, among others, and found out that the ERT was efficient in evaluating the moisture content of the waste in the landfill. In a study conducted by [10], they also attempted the assessment of soil compactness from geophysical datasets and obtained promising results from the use of the approach. Similarly, [11] used the common-offset technique of the GPR (ground-penetrating radar) for soil water contents to ascertain the performance of the empirical petrophysical relationships documented in the literature. Laboratory-scale dynamic hydrologic events assessment in the unsaturated zone using a sand tank was evaluated from the time-lapse multi-offset GPR dataset [12]. Furthermore, the review of [13], on the versatility of GPR, particularly in the area of its application to soil water content evaluation, gave further courage to the choice of the approach in this study. Moreover, Lombardi et al. [14] reported the applicability of GPR as a technique for providing valuable information about the spatial variability of soil within-field properties such as state variables that are relevant to precision agriculture. Reviews by [15,16] have shown the versatility of the GPR techniques in subsurface media characterizations, particularly in the near-surface that hosts the emplacement of the soil. They submitted that the non-invasiveness of the GPR method has allowed the possibility of sensing soil layers and the process of soil water content evaluation, as well as delineating the root mass of plants. Evaluation of soil bulk density from tillage operations using the GPR data was documented in the work of [17]. Hubbard et al. [18] concluded in their research that soil texture influenced shallow soil water content. According to [19,20,21], the characteristics of the environment of EM waves influence the shape and magnitude of the amplitude of the transmitted pulse. Thus, the signature of the GPR signal may reflect the nature of propagating media. Salat and Junge [22] stated in their report that the relative permittivity of soil is linearly related to its porosity and compaction.

The chosen geophysical method, ground-penetrating radar (GPR), utilizes electromagnetic (EM) pulse energy via antennae sensors, which enhances continuous data sampling, hence promoting spatial scanning in field measurement. Soil is made of solid matrix phase, pore spaces, and liquid phase. The volumes of phases of the earth materials (soil inclusive) greatly influence its electrical properties (conductivity and dielectric constant) [23]. This is because water has a high dielectric constant of about 81 [24], and the solid matrix of earth materials dielectric is between 3 and 30 [24]. This connection between components of earth materials and the electrical properties on which the principle of GPR is based suggests the applicability of the method and was thus put to test in this study. It is envisaged that the signal responses of the transmitted GPR pulse may give insight into the difference in the energy level of the water-holding capacity of soil at the vadose zones.

The expectation is that the results will provide a better understanding of the water-holding capacity in the unsaturated zone, and the water movement direction can also be ascertained to a certain extent. This information can be very useful in the soil layer characterization and proper planning by agriculturists, particularly in the area of irrigation and some other engineering applications. Thus, the study outcomes may be a significant leverage to soil assessment and provide meaningful guides to soil utilization, management, and preservation.

2. Materials and Methods

What constitutes effective data collection is the systematic method of gathering and measuring information from a variety of sources to obtain comprehensive details of the area/subject of interest. Following the review of literature on the different geophysical methods, particularly on their suitability for various projects, terrain, and targets, confidence in the deployment of the GPR was built. The ground-penetrating radar technique is an electromagnetic (EM) geophysical method that utilizes EM in the range of 10 MHz to 3000 MHz [25] as a source of energy that propagates through subsurface media. The transmitter of the GPR system radiates EM pulse energy into the underlying media of investigation and the reflected and refracted energy responses as a result of encounters with buried objects and or a boundary between media of different permittivity levels (dielectric constant). The responses are intercepted and recorded by the receiving antenna (Figure 1a). Details of the principles and theory of GPR are found in [19,20,26]. Test measurements were made on a relatively flat grassland loam soil at the Kryspinów area on the outskirts of Krakow, Małopolskie, Poland. Field data were obtained at the beginning of spring of the test year during fair weather conditions, which allowed the field measurements in idealized conditions.

Field measurements were conducted using a ProEx model system unit manufactured by MALA Geoscience Sweden (now: ABEM/MALA) with the antenna of central frequency at 800 MHz. The choice of this frequency was based on its short wavelength and thus the resolution it can deliver at a relatively few meters depth of penetration within the unsaturated zone of soil, which was the focus of this study. Field data were acquired in the constant offset (reflection profiling) along the established two transects (A–A1 and B–B1) that are orthogonal (Figure 1b). Field acquisition parameters (

Table 1. Field data acquisition parameters.

S/N	Field Measurement Input Parameters	Title 3
1	Antenna frequency	800 MHz
2	Trace interval/station spacing	0.010 m
3	Sampling frequency	8000 MHz
4	Number of samples	500
5	Stacking times	16
6	Time window	60 nS
7	Antenna spacing	0.14 m

) were input into the GPR system setup, and the entire system unit was mounted, with data acquired by pulling a wheeled antenna along the transect at a walking speed. Measurements were made at different soil conditions, first at their natural state and later when a known volume (25 L and 50 L at different time intervals) of water was added to the soil surface along the transect path, which is relatively flat over a set period. The added water was to simulate the infiltration mechanism within the soil horizons and monitor the infiltration pattern and direction, which theoretically are controlled by potential energy within the matrix of the soil particles. The added water was made at two consecutive intervals of 3 h, and infiltration was allowed for 1 h before GPR measurements were taken. A few multi-offset measurements, namely, wide-angle reflection and refraction (WARR), were also made for EM velocity determination. WARR is a measurement technique that involves keeping the transmitter antenna at a fixed point while the receiver antenna is shifted at varying continuous offsets from the transmitter antenna [19,20]. WARR data allow for the obtaining of the EM signal travel time as a function of the antenna separation offset distance. Evaluated information can then be used to calculate the velocity of the EM at the point that may facilitate depth determination [16]. After the field data acquisition, the data were edited and presented as radargrams, which represent the stacked plot of the signal amplitude against travel time using REFLEXW version 8.5 software developed by Sandmeier Inc., Karlsruhe Germany [27], as shown in Figure 2a.

The field data were subjected to post-acquisition processing to enhance signal-to-noise ratio and distinguish clutters from signals, which aided in providing better interpretations and clarifying targets. GPR data like other geophysical data are usually subjected to processing due to accompanied noise and clutters [26]. With the aid of the REFLEXW software, the following processing steps were carefully performed on the field data. ‘Time zero’ correction was performed to give the start time specified in the file header of the recorded data [27]. This was to ensure accurate signal travel depth estimation. Processing that enhances signal-to-noise ratio includes ‘dewow’, a filtering technique that suppresses the low frequency. The ‘DC-shift’ subtraction filter enhances the proper alignment of the signal shift from the central axis due to equipment. Also, ‘background removal’, a 2D filtering process that removes the elements of the common environment to show clearly the anomaly present, was also performed on the field data. Finally, a ‘time gain’ application was also performed on the data to boost the energy of the signal that may have been attenuated with time as it travels deeper into the subsurface media. Despite the field data processing, inherent information about the propagated media may still be hidden due to the nature of the GPR data, which are non-stationary. The non-stationary nature of the GPR EM wave allows its signals to have varying frequency components with time and hence makes the characterization of the embedded features of the signals’ locations and in the time frame usually unknown [28,29]). Therefore, to minimize ambiguities that may arise from interpreting the data only in the time domain in which it was recorded, the signals’ attributes computation may aid subtle feature delineation according to [30]. Similarly, spectra estimations of the processed data were performed for the analysis of the frequency spectra components of the recorded signals.

2.1. Signal Attributes

Signal attributes as used in seismic and GPR data are enhancement mathematical computations utilized to improve the anomaly of targets within the subsurface on the data plots. They are found in the instantaneous parameters of the signals. Principles of instantaneous parameters correlate with geometry and the physical property changes of the material by which the GPR signal travels [31]. The basic purpose of the attributes is to reveal subtle features by identifying the same trend or pattern by quantifying particular characteristics. Thus, sought targets are enhanced by removing or hiding some portion of the signal. Most signal attributes are derived through filters, statistics, and transforms [32]. The commonly used attributes in GPR are the instantaneous variation of different parameters, which are computed through Hilbert transform [33]. Instantaneous amplitudes, frequency, and phase are mostly computed by the transformation of the real data by the Hilbert function. What the Hilbert transform does is a filtering function that goes through the amplitude of the spectral components uninterrupted but alters the phases’ spectral components by ninety degrees [28].

Basic relations between instantaneous frequency and seismic signal attenuation have been given in the work of [33]. The author considered that geophysical signal (e.g., seismic or GPR time traces) s(t) could be represented by its envelope a(t) and phase ϑ(t), and thus

s(t) = a(t)·cos(ϑ(t))

(1)

and also the quadrature (imaginary) trace of the same signal may be represented as

s*(t) = a(t)·sin(ϑ(t))

(2)

and the complex (analytical) trace z(t) may then be given by the following equation:

z(t) = s(t) + i·s*(t) = a(t)·ei·ϑ(t)

(3)

The work concluded that the quadrature (imaginary) trace is the Hilbert transform of the real trace. It is obtained by the phase shifting of the recorded trace by 90 degrees. The complex trace comprises the real trace (recorded) and the imaginary trace [34]. Once the quadrature trace is found, the instantaneous amplitude a(t) and the instantaneous phase ϑ(t) can be expressed as

a(t) = [s(t)2 + s*(t)2]1/2

(4)

and, respectively,

ϑ (t) = tan − 1[(s*(t))/(s(t))]

(5)

Further details of other derivative attributes such as frequency are found in [33,35].

Computed instantaneous attributes have different applications depending on the sought targets. For instance, instantaneous amplitude outputs the envelope of the selected data at the same location, which enhances lateral variations within events [36]. The attribute considered for the interpretation of the field data in this study is the instantaneous phase. It calculates phases at the same location, which emphasizes the spatial continuity/discontinuity by providing a way for weak and strong events to appear with equal strength. Essentially, it responds mostly to dielectric changes [24], which can be used to distinguish subsurface material properties. The instantaneous phase has no amplitude information but relates to the phase component of the wave propagation [37]. Hence amplitude suppression due to earth convolution and attenuation may not directly influence phase attributes. It is on this premise that the phase attribute was computed on the field data to delineate the changes that accompany the soil horizons as water infiltrates through its matrix, which can be a good pointer to its potential status. Mathematically, the instantaneous phase is expressed as shown in Equation (5). The attributes were computed using the appropriate module in the REFLEXW software.

2.2. Spectral Analysis of Field Data

Discrete-time data, which are usually in sequence and contain a series of data points, are indexed in time order. Recorded GPR signal data that follow this trend are classified as time series data. They are considered to be taken at successive equally spaced points with the change of time. Thus, their representations are frequently plotted on curves in the time domain such that details culminating in the materials that gave rise to such signals are not known. Although a series of processing such as filtering, correlation, and enhancement may improve the analysis of the data, details about the inherent information may still be concealed. Hence, the concept of spectral analysis may provide insight into the variations in properties of what make the recorded data [38].

Spectral analysis involves several methodologies that evaluate time series data in terms of frequency contents. It commences with the breaking up of the data often referred to as ‘spectral decomposition of time series data’, which is the partitioning of the variance of the time series in terms of frequency [38]. There are various methods used to decompose time series data that are computed on the waveform of the data and, in the case of GPR, on the traces [39]. The most commonly used include Fourier transform, Wavelet transform, Wigner–Ville distribution, and Empirical Mode decomposition. Following the decomposition of the signal, the analysis of the frequency components is made, which may reveal details about the signal. Wavelet transform (WT) was used for the decomposition of the processed field GPR signal data in this study. This was informed by the better performance of such transform in non-stationary signals such as the GPR signal [40].

2.2.1. Wavelet Transform

Wavelet transform represents the projection of a signal to a group of functions called wavelets, which offer localization of the signal in the frequency domain. Wavelet transform enhances frequency resolution at low frequencies and high time resolution at high frequencies [41]. A function x(t) will have a wavelet transform in the form of [42]

$$X(a,b)={\displaystyle \frac{1}{\sqrt{a}}}{\int }_{-\infty }^{\infty }\overline{\psi }\left({\displaystyle \frac{t-b}{a}}\right)x(t)dt$$

(6)

where $a$ is the scaling factor and $b$ is the time shift factor. Ψ(t) is a continuous function in both the time and frequency domain called the mother wavelet.

The convolve mother wavelet generates a source function that produces a daughter wavelet that is a scaled and translated version of the mother wavelet [42].

2.2.2. Power Spectral Density

The power spectral density (PSD) of the signals was also computed to show the representation of the power contents of the signals against the frequency. It is a means of characterizing the signal based on the average power of the signal sinusoid frequency component variations [43]. A signal x(t) has power spectral density in the form of [44]

$$\mathrm{P}\mathrm{S}\mathrm{D}={\int }_{-\infty }^{\infty }{\left|x\left(t\right)\right|}^{2}dt={\int }_{-\infty }^{\infty }{\left|\widehat{x}\left(f\right)\right|}^{2}df$$

(7)

where $\widehat{x}={\int }_{-\infty }^{\infty }{e}^{-2\pi ift} x\left(t\right)dt$ is the Fourier transform of the signal and f is the frequency (Hz).

Both the wavelet transform and the power spectrum computations were executed in MATLAB software version 2023a developed by Math Work Inc., Natick, MA, USA using the appropriate functions [45].

3. Results

After the field data processing and the computation of the necessary signal attributes, the analysis of the results and the discussion are presented in this section. GPR sections from a profile and corresponding calculated instantaneous attribute (phase) are displayed in Figure 2a,b respectively. Figure 2a is the representation of the recorded GPR pulse EM energy amplitude variations with time along the profile. It is the response of the propagated soil media characteristics to the transmitted EM energy.

3.1. Instantaneous Phase Attribute Result

The result of GPR signal instantaneous attribute computation before and after water infiltration gave clear indications of the parts of the section that have an affinity towards the infiltrated water within the identified horizons (Figure 3). There is an obvious sudden change in the horizons feature after the infiltration of the water with marked sharp contrast (Figure 3, inset).

Traces of instantaneous phase plots at zones tagged ‘high and low’ MP (matric potential) (Figure 4a) were compared. The plots also reveal the discrepancy in the two zones using the calculated attributes. Phase magnitude was higher at the suspected high matric potential. There was an increase in the magnitude of phase at the low MP zone in comparison with at the high MP when water infiltration was conducted.

3.2. Signal Spectral Analysis

Spectral analysis of the signals at the different zones of consideration is presented in Figure 4b. It is a representation of the power spectral density of the GPR signals at the different areas of consideration. It describes the power present in the signal as a function of frequency, per unit frequency [46]. The magnitude of the lower spectral was higher at the suspected high matric potential, while the reverse was the case at the high-frequency spectral.

Time-frequency plots for the GPR signal at various soil conditions are shown in Figure 5. Before water infiltration, low-frequency spectral evolved longer to the time ~22 ns at the suspected high matric potential. In contrast, the same frequency spectral was short-lived to ~11 ns at the suspected low matric potential (Figure 5a). This pointed out the difference in the dielectric properties at the different zones, which may have been influenced by the variation in degree of water saturation. However, with water infiltration, the time-frequency plots show a significant influence of the water presence at the suspected low-potential zones than the high-potential zones (Figure 5b,c).

The centroid frequency plot (Figure 6) also corroborated the delineation of matric potential variations within the considered zones in the subsurface. Essentially, the centroid frequency is a spectral decomposition technique that allows the location of the center mass of the signal spectrum as a function of time with distance [47]. It is computed using the weighted mean of the frequencies present in the signal. It is determined using S-transform, with their magnitudes as the weights [48] to represent changes in the propagation condition of the signals. Thus, the value of centroid frequency decreases/increases as the signal enters a high/low attenuation area. In the present scenario, the centroid frequency decreased with depth at the suspected high matric potential (dotted rectangle labeled as High MP in Figure 6) and increased with depth to a depth of about 0.6 m at the low matric potential.

The centroid frequency computation allows a glance at the state variable of the media through which the pulse EM energy has propagated, and its integration with other forms of data analysis enhances the fidelity of the resulting outcome.

4. Discussion

The GPR sections enhanced the delineation of the uppermost horizons (blue line in Figure 2b) within the scanned distance. The phase attributes enhanced the identification of the soil layers based on both vertical and lateral distortions on the continuity of horizontal patterns in the sections. Integration of the enhancement with the amplitude reflections (Figure 2a) in the processed data facilitated the recognition of the soil layers. A better soil horizon architecture was made possible by the GPR signal phase attributes. The identified soil layers included the topmost horizon, which is thought to be composed of high organic contents, particularly from the field observations. Based on literature records, the adjacent layer is considered loam sandy soil. The concealing nature of the adjacent layer and the dearth of exposed soil profile in the area did not allow us to substantiate this assertion. In addition, the depth to the horizons is also determinable from the GPR sections, which may provide fair information about the thickness of the layers. The zone identified as the low matric potential is marked with a cyan dotted line (Figure 3b). It is suggestive of the zone of water movement direction after the water infiltration, indicating lower potential that triggers the water movement. The phase attribute also enhanced the delineation of adjacent soil horizons suspected to have higher matric potential (white dotted line in Figure 3b). Moreover, the increase in the magnitude of the instantaneous phase computed for the GPR signals is thought to be the effect of the infiltration that informed the affinity of the zone to water. Furthermore, the discrepancy in the magnitude of the lower spectra at the different zones of consideration may indicate the influence of water content variations in the two scenarios. Noticeable power peak shifts from high frequency to low frequency in the high matric potential zone may also be an indication of variations in the volumetric contents of the two zones. In addition, changes observed from the time-frequency plots of the recorded GPR EM signals from the test site may also explain the tendency of movement of infiltrated water towards the low matric potential rather than the high potential zone. Finally, the centroid frequency plot analysis also gave credence to the soil characteristic variations identified from other results analyses. Nevertheless, it is worth pointing out that the non-uniqueness characteristics of geophysical data interpretations may also inform the possible explanation of the obtained results to be stratigraphical changes in the vadose zone. The results of the frequency analysis may also be related to changes in near-surface material and antenna coupling. It should be noted that the GPR technology application in subsurface evaluation is largely site-specific. The soil conditions (clay contents, levels of water saturation, and accessibility) need to be considered before the use of the method. This is to ensure high resolution, accuracy, and depth of penetration of the EM pulse energy due to attenuation.

5. Conclusions

A non-conventional technique for soil investigation was adopted in this study to delineate the soil inherent characteristics of the study field and also to prove the possibility of the adopted approach to circumvent the traditional methods’ pitfalls. Analysis of the results obtained from the field data provided significant insight to map out the soil horizons and the matric potential of the identified layers. The matric potential is the driving force of the amount of water held within the layer, an indication of understanding the rate and direction of water movement within the layer, and a pointer to the possible particle sizes of the grains of the horizons. Results of the complex mathematical manipulation of the field data also assisted in the delineation and distinguishing of the adjacent layers, which are relatively thin and may have been obliterated by the overlying and underlying layers’ characteristics. Although the conventional direct measurement of the matric potential was not made at the time of the study to complement and confirm the veracity of the approach, the research outcomes show that the chosen technique, GPR, can be used to study both subtle and overt soil properties and at the same time allows quick and repeated evaluation while overcoming some problems associated with the traditional soil survey methods. Thus, data from such an evaluation may be valuable input in precision agriculture where exact details are sought.