Introduction

Soil is so important to continuous existence of human life and needs to be assessed because of its non-renewable characteristics. Moreover, it is also a common knowledge that the Earth’s population is increasing geometrically and basic products and necessity of life have their roots in the soil. (Weil and Brady 2017) postulate that World population may be stabilized at 9–10 billion in the Twenty-first century. However, real extent of soil available to meet the demands will not increase. Hence, it is pertinent to maintain and manage the available ones. Human activities particularly in the area of agricultural activities have significantly caused degradation to soil due to quest for more crop yields. Mechanization which has helped in achieving more yields is continuously affecting the natural settings of the soil. For instance, the exposures of soil to vehicular traffic impact affect its physical properties such as bulk density, porosity and strength and consequently influence the chemical properties (Nawaz et al. 2013). This has been found also to have direct bearing on the fertility of the soil. In order to understand soil characteristics for its effective utilization, its assessment cannot be overemphasized.

Traditional soil investigations involve soil sampling using various devices such as hand auger, coring, and logging through sections as road cuts and erosional sides. Thus, collected samples are analyzed at the laboratory for the evaluation of soil physical properties (Bai and Jin 2016). Moreover, some soil chemical properties can be estimated using core samples analyzed in the laboratory. Infiltrometer, borehole logging and neutron probes techniques data may also revealed soil information. Also time domain reflectometry techniques can be used to determine the water contents of soil horizons (Hubbard et al. 2002). These approaches can provide the soil attributes but only at discrete intervals which are interpolated. Essentially, they rely on interpolation of point information that has resulted in inadequacies of the methods. Moreover, the concealing nature of the soil horizons makes it difficult to delineate changes in the physical and chemical properties over time. However, to circumvent the shortcomings, geophysical methods can be used particularly with the advent of advances in technology and instrumentations (Everett 2013). The ground penetrating radar (GPR) geophysical method was selected in this study to assess soil horizons and its compaction effects caused by traffic. Noninvasiveness which promote in-situ evaluation, fast and repeated measurement possibilities are what informed the choice of the method. Muniz et al. (2016) studied the compaction effects on soil of a grazing field using GPR and found out that the impact was significant on top horizon. Hard pan—an indurated sublayer of the soil—was delineated within a section of a test site using GPR as reported in the work of Raper et al. (1990). Characterization of tillage effects on soil using GPR and electromagnetic geophysical methods shows the effective applicability of GPR in identifying soil horizons affected by different tillage practices (Jonard et al. 2013). The GPR method has ability for continuous spatial data acquisition which gives better resolutions. Thus, interpolation of point data that characterized conventional methods is ameliorated.

The study utilizes GPR geophysical method for assessment of subsurface parts of the soil. Main objectives are delineation of the soil horizons and compaction effects on its matrix arrangements due to traffics. It is expected that such findings provide basic information on physical properties and characteristics settings of soil that may guide its utilization and management.

Method of study

This section discusses the approach adopted and brief underlying principles of the data measurement method used in the study. It also highlights the test location and its condition in order to give better perspectives on which the study was carried out.

Ground penetrating radar (GPR)

GPR is a geophysical method of investigation that utilizes electromagnetic pulse energy in the frequency range of 10 MHz to 3 GHz (Forte et al. 2014). The transmitter component of the GPR system allows the passage of generated pulse energy. The energy then propagates through the material of investigation, and interactions with the material response are sensed by the receiver component. Typical GPR system assemblage and mode of operation are as illustrated in Fig. 1.

Fig. 1
figure 1

Illustration of GPR system units and EM pulse energy paths

Full size image

Responses from the receiver on which the subsurface imaging of GPR technique is based is a function of the reflectivity of the geological medium or soil components under investigation (Everett 2013). For instance, GPR signals that propagate via a non-magnetic subsurface media of contrasting dielectric constant soil layers will give reflectivity response according to the degree of contrast. Details of principles and theory of GPR are found in (Daniels 2004; Jol 2008; Annan 2003).

Field test measurements were taken at a farmland near Krakow that is composed of loam sandy soil with relatively flat topography. Compaction was induced on parts of the soil using movement of tractors (average weight of 1750 kg) from 2 to 10 passes simulating farming activities. Subsequently, 14 profiles were established in a rectilinear geometry of 70 m by 10 m such that it cut across and along the tractor’s movement tracks (Fig. 2a). Both constant offset (profiling) and wide angle reflection and refraction (WARR) techniques were applied in the data acquisition. Sketch of the profile points, WARR measurements and excavation for soil profile are shown in Fig. 2c. The measurements were taken using ProEx model system unit (Fig. 2b) manufactured by MALA GeoScience Sweden with antenna of central frequency 800 MHz. The choice of this frequency was based on its wavelength that improves resolution and penetration depth to the unsaturated zone of soil which was the focus of this study. All measurements were made along the established profiles (Fig. 2c). After the GPR system unit set up and configuration, measurement input parameters (Table 1) were entered. Signal quality was assured and confirmed. Then, the entire system arrangement was mounted and data taken by pulling wheeled antenna along the profiles at walking speed to ensure adequate data acquisition and stacking.

Fig. 2
figure 2

a GPR data test site and profiles. b Wheel 800 MHz antenna in constant offset mode. c Profiles layout

Full size image
Table 1 Field data acquisition parameters
Full size table

WARR mode measurements were taken for estimation of EM wave velocity used in time-depth conversion of the field profiling data. This is to ensure accurate deduction of depth of delineated targets on the GPR section plots. Basically, WARR measurement is a multi-offset GPR technique that involves varying antennae spacing. At a fixed location, transmitter antenna is kept at the spot while receiver antenna is moved away at increase interval. Thus, measuring the two-travel time of the EM energy at the different offset can be represented as time–distance (offset) curves as shown in Fig. 3b.

Fig. 3
figure 3

1-D velocity model calibration from WARR GPR data. a One-dimensional velocity distribution with the Vrms, (b) WARR data plot with the direct waves (air and ground) and reflection hyperbola (black lines). c Velocity/time spectrum evaluated via semblance computation

Full size image

The field profiling data were edited and presented as radargrams which represent the stacked section of the signal amplitude along profile distance against two-way travel time. The plotting was made in REFLEXW version 8.5- software developed by Sandmeier Inc. Germany (Sandmeier 2017a, b) as shown in Fig. 4. Subsequently, the data were subjected to post acquisition processing. This was performed in order to enhance signal-to-noise ratio and distinguish clutters from signals and thus improves targets clarity and better interpretations. GPR data like other geophysical field data are usually accompanied by noise and clutters (Daniels 2004). 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 (Sandmeier 2017a, b) otherwise, depth estimation may not be accurate. Processing that enhanced signal-to-noise ratio includes “dewow”- a filtering technique that suppresses low frequency; “subtract DC”-shift filter enhances proper alignment of signal shift from central axis due to equipment. Others include background removal- a 2D filtering process that removes elements of common environment to show clearly the anomaly present. And time gain 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.

Fig. 4
figure 4

GPR radargram of some profiles from the test site (a) Profile 1 and (b) Profile 2

Full size image

Results and discussion

The results of the WARR datasets analysis were used to compute the velocity of the EM wave and aided the time-depth conversion. Figure 3 shows the semblance velocity analysis and 1D model used in the time-depth conversion. The velocity data were computed based on the principle of correlation-like technique of the arrangement of the traces and DIX equation (Taner and Koehler 1969; Yilmaz 2001). It was executed in the CMP module of the REFLEXW software. The DIX equation is an expression that shows the interval velocity as a function of the stacking velocity (root mean square velocity) and the time intervals as:

$$V_{{\text{int}}} = \frac{{\left( {V_{{{\text{rms}}{.}n}} } \right)^{2} t_{n} - \left( {V_{{{\text{rms}}.n - 1}} } \right)^{2} t_{n - 1} }}{{\left( {t_{n} - t_{n - 1} } \right)}}^{1/2}$$
(1)

(Reynolds 2011), where Vint is the interval velocity, Vrms.n, and Vrms.n-1 are root mean square velocities and tn and tn-1 are reflected ray two way travel times for the nth and (n-1)th reflectors, respectively.

The DIX analysis enhances layer-by-layer evaluations and determination of both thickness and velocity for individual horizons.

Semblance plot is a measure of the similarity between large several trial hyperbolas and the real hyperbola on the multi-offset data (Everett 2013). It represents the tracking velocity versus the apparent depth (Fig. 3c) where the best tracking velocity at a depth relates to the highest semblance value. The EM wave different paths to and from reflector in the subsurface as antennae offset changes resulted in the different curves in Fig. 3b. The straight line that passes through the origin is the direct airwave while the hyperbolae are the reflected wave. If the hyperbola reflection is fitted, a depth–velocity model may be built to give information about the velocity distribution (Takaharshi et al. 2012).

One-dimensional velocity model (Fig. 3a) was built simultaneously corresponding to points of maximum semblance or energy (black curves and black asterisk) on the distance–time and velocity–time spectrum (Fig. 3b and c, respectively) (Yilmaz 2001). This 1-D model represents the composites velocities of all the horizons above each reflecting horizons which was used to compute the 2D velocity data that served as input in the time-depth conversion for the GPR profiles data.

From the depth estimate, the velocity of the EM wave within any identified horizon can be known and is related to the dielectric constant according to Eq. 2 (Annan 2003).

$$V = \frac{c}{\sqrt \varepsilon }$$
(2)

where c is the velocity of light (~ 0.30 m/ns) and ε is the dielectric constant.

The GPR signal response is dependent mostly on the reflection coefficient that is related to the contrasting dielectric properties of the traversing media according to Eq. 3 (Daniels 2004).

$$R = \frac{{\sqrt {\varepsilon_{1} } - \sqrt {\varepsilon_{2} } }}{{\sqrt {\varepsilon_{1} } + \sqrt {\varepsilon_{2} } }}$$
(3)

where ε1 and ε2 are dielectric constants of the different layers of the media.

A 2D display of some of the profiles is shown in Fig. 4. The horizons are not well discernable distinctly due to suspected relatively thin layer. It is possible that the thin layer is subtle to the raw GPR signals but may be more observable via the signal attributes. Signal attributes, particularly the wave form (seismic and EM) types, are computed quantities obtained from the field data using mathematical relationships. They are mainly used for data visualization and integration purposes (dGB Beheer 2019). Akinsunmade et al. (2019 g) and Tomecka-Suchoń and Marcak (2015) suggested in their works that they may be employed for identifying indistinct subsurface features from GPR data as well as delineation of lateral variations within events. Hence, an attribute (instantaneous phase) of the field data was calculated. It enhances delineations of layers and other subsurface anomalies that were shrouded in the raw data (Fig. 5a). Instantaneous phase calculates phase at the same location which emphasizes the spatial continuity/discontinuity by providing a way for weak and strong events to appear with equal strength (Barnes 1991). GPR time trace s(t) can have phase as:

$$\vartheta \left( t \right) = \tan^{ - 1} \left[ {\frac{{s^{*} \left( t \right)}}{s\left( t \right)}} \right]$$
(4)
Fig. 5
figure 5

a Calculated instantaneous phase attribute of Profile 3 with inverted black triangle showing number of traffic passes. b Excavated part of test site showing (i) Soil profile tied to (ii) a GPR section (profile1)

Full size image

(Barnes 1991), where s*(t) is the quadrature (imaginary) trace which is the Hilbert transform of the real trace.

Essentially, instantaneous phase attribute responds mostly to dielectric changes (Reynolds 2011) which can be used to distinguish subsurface material properties. Instantaneous phase has no amplitude information but relate to phase component of the wave propagation (Taner 1992).

The attribute was computed using the complex analysis module in Reflex software.

Three major horizons (Fig. 5b) are distinguishable with the uppermost as the topsoil and the beneath layers as subsoil horizons as discerned from the soil profile log at an excavated portion around the field test site. Excavated portion allows the observation of the soil profile to the top of the subsoil. Observed logged soil profile was tied to the closest GPR section (Profile 1) for correlation (Fig. 5b). The contrasting horizons may be due to the variations in soil matrix arrangement, water saturation, organic contents, porosity and degree of compaction.

The calculated attributes of the field data gave better resolution of the horizons. It can also be delineated that there are variations in the thickness of the identified layer.

It is also visible from the GPR signals (Fig. 6) that there is distortion in the continuity of the horizons at the points of traffic (Inverted Black triangles). This may be compactness effects of the soil layers which are thought to have affected the soil matrix arrangement and thus reduces the volume of voids (porosity). It may have further influenced the volumetric water contents which control the dielectric property. At the few passes of the traffic (1 and 2), traffic effects are not so significant (dotted white eclipse Fig. 6a). However, effects of the traffic are pronounced and thus discernable at the several passes (9 and 10). The subsoil horizons at these points are visibly reworked as shown in Fig. 6b (dotted white rectangle) which virtually obliterates the adjacent horizons.

Fig. 6
figure 6

Zoomed-in plots of Instantaneous attribute of Profile 3: (a) section with few traffic passes points (inverted black triangle) and (b) section with several traffic passes points (inverted black triangle)

Full size image

Furthermore, analysis of traces at the traffic paths and surroundings with no traffic paths show the influence of the impacts of the passes on the subsurface horizons particularly the topsoil (Figs. 7a and b). There are signal reflection magnitude variations (largely dependent on reflection coefficient (R) Eq. 2) at the different passes and non-traffic points. Figures 7a and b also show the thickness variations which is more significant at the ten passes than at the one pass of the traffic (dotted rectangles Figs. 7a and b). It is also observable from the trace plots (Figs. 7a and b) that the EM signals traveled deeper at the traffic zones (~ 1.5 m) than the non-traffic areas (~ 1.0 m). This may be attributed to change in the velocity which is a function of the material properties of the propagating media that may have been influenced by traffic. It may be inferred that signals at the traffic zone are relatively less attenuated and thus may suggest a way to recognize compacted soil portion on GPR section.

Fig. 7
figure 7

GPR trace plots at different numbers of traffic passes (a) one pass and (b) ten passes with corresponding point of no traffic

Full size image

Figures 8a and b show profiles of GPR scans along a traffic path and a parallel approximate equal length where there was no traffic. The uppermost horizon in Fig. 8a (area within dotted lines) depicts alterations that makes it appears as distinctive layer of equal thickness. It is may be attributed to the equal vertical exerted impact of the traffic movements on the soil layer at that depth. However, profile shown in figure (b) reveals clutters/hyperbolas (red circles) in the same horizon with irregular continuity of reflections. Although the profiles are adjacent to each other, yet significant different is observable. The discrepancy may be due to compaction effects of the traffic which may have altered the shape of objects (pebble, gravels) within the soil layer that may have given clutters and hyperbolas signatures.

Fig. 8
figure 8

GPR radargrams. a Profile along traffic and (b) profile parallel to (a) without traffic

Full size image

A 3D representation of the recorded GPR profiles was also attempted. It is worth to note that the 3D plot axes have been adjusted automatically by the software program used in order to accommodate all data points in the cube.

The data display in 3D form (Fig. 9) also gives a better planar spatial and temporal variation within the scanned site pointing to the heterogeneity of the delineated soil horizons. Time slice (Fig. 9b) shows spatial variation at a time cut providing information on the heterogeneity of that horizon of the soil. Such variability information may be useful in agricultural practices.

Fig. 9
figure 9

3Dimension plots of the field data. a 3D cube, b depth slice cut at 0.3 m and c XYZ axes cuts showing traffic effects at one pass and ten passes at ~ 7.5 ns

Full size image

Of particular interest is the display of the depth extent of the traffic compaction at the highest and lowest numbers of passes (Fig. 9c). Significant impact of the traffic is seen at timeslice 7.5 ns which is about 0.5 m depth but not noticeable at same depth at the 1 pass point (right and left sides dotted eclipse, respectively, Fig. 9c).

Conclusions

Efforts made in this study were to assess effects of compaction induced by traffic at different passes on identified soil horizons. The results have shown delineation of soil horizons with visible distorted continuity at traffic points. It is thought that the soil matrix arrangements and its contents that primarily control the reflectivity coefficient in the measured data may have been altered at these points. Thus, the dielectric intrinsic property of the soil may have been influenced. Further data analysis shows that densification of the soil layer by the traffic may have caused EM velocity change as it propagates through the media. This is observable as reflections variations on the field data trace amplitude plots at the traffic points.

It thus suggests that GPR scan can sense compacted zones of soil and also delineates pockets of indurated (hardpan) portion that may intercalate the horizons. The hardpans significantly have adverse effects on plants’ roots penetration and overall yields of crops. Delineating such zones within the near-surface subsoil horizons can inform what measure to be taken in order to make effective use of the soil.

Delineation of zones of compactness within soil layers may be recognized on the radargrams as an abrupt change to the lateral continuity which may be different from fractures or buried objects. The zones may tend to give rise to ringing effects on the adjacent beneath layers if degree of compaction is high. Ringing effects occur when overlaying layer reflections result in signal block to adjacent layers (Utsi 2017). However, core samples to the depth of the zones that can be analyzed in the laboratory will be complementary confirmation of presence. Although there might be variations in the material properties of the propagating media, yet GPR can be used to characterize subsurface horizons with associated impacts of traffic.