Next Article in Journal
Ultra-Light Aircraft-Based Hyperspectral and Colour-Infrared Imaging to Identify Deciduous Tree Species in an Urban Environment
Next Article in Special Issue
Detection and Identification of Remnant PFM-1 ‘Butterfly Mines’ with a UAV-Based Thermal-Imaging Protocol
Previous Article in Journal
Modelling the Vertical Distribution of Phytoplankton Biomass in the Mediterranean Sea from Satellite Data: A Neural Network Approach
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Laboratory Measurements of Subsurface Spatial Moisture Content by Ground-Penetrating Radar (GPR) Diffraction and Reflection Imaging of Agricultural Soils

1
Department of Civil Engineering, Faculty of Engineering, Ariel University, Ariel 40700, Israel
2
School of Geosciences, Faculty of Exact Sciences, Tel Aviv University, Tel Aviv 6997801, Israel
3
Geo-Sense Ltd., Netanya 4266010, Israel
4
Leon H. Charney School of Marine Sciences, University of Haifa, Haifa 3498838, Israel
*
Author to whom correspondence should be addressed.
Remote Sens. 2018, 10(10), 1667; https://0-doi-org.brum.beds.ac.uk/10.3390/rs10101667
Submission received: 7 August 2018 / Revised: 3 October 2018 / Accepted: 19 October 2018 / Published: 22 October 2018
(This article belongs to the Special Issue Recent Advances in Subsurface Sensing Technologies)

Abstract

:
Soil moisture content (SMC) down to the root zone is a major factor for the efficient cultivation of agricultural crops, especially in arid and semi-arid regions. Precise SMC can maximize crop yields (both quality and quantity), prevent crop damage, and decrease irrigation expenses and water waste, among other benefits. This study focuses on the subsurface spatial electromagnetic mapping of physical properties, mainly moisture content, using a ground-penetrating radar (GPR). In the laboratory, GPR measurements were carried out using an 800 MHz central-frequency antenna and conducted in soil boxes with loess soil type (calcic haploxeralf) from the northern Negev, hamra soil type (typic rhodoxeralf) from the Sharon coastal plain, and grumusol soil type (typic chromoxerets) from the Jezreel valley, Israel. These measurements enabled highly accurate, close-to-real-time evaluations of physical soil qualities (i.e., wave velocity and dielectric constant) connected to SMC. A mixture model based mainly on soil texture, porosity, and effective dielectric constant (permittivity) was developed to measure the subsurface spatial volumetric soil moisture content (VSMC) for a wide range of moisture contents. The analysis of the travel times for GPR reflection and diffraction waves enabled calculating electromagnetic velocities, effective dielectric constants, and spatial SMC under laboratory conditions, where the required penetration depth is low (root zone). The average VSMC was determined with an average accuracy of ±1.5% and was correlated to a standard oven-drying method, making this spatial method useful for agricultural practice and for the design of irrigation plans for different interfaces.

1. Introduction

Monitoring soil moisture content (SMC, i.e., water content) is important for determining water-use efficiency and planning precision-agriculture programs. Optimizing available water and land resources by adjusting the correct SMC to the specific agricultural crop will maximize crop yields and prevent damage. Optimization serves both the economy and the environment and will increase field productivity and decrease irrigation expenses, which can be enormous, especially in arid and semi-arid zones [1].
Soil moisture serves as a binding factor for soil materials, thereby affecting soil structural stability and strength. Soil moisture is involved in the processes of soil development and degradation with respect to both chemical and biological aspects, and agricultural crop production is influenced mainly by water availability [2].
Today, SMC is monitored by conservative local sampling and testing methods in the field and in the laboratory, which are time-consuming and expensive, making use of neutron gauges, various types of tensiometers, and time-domain reflectometry [3,4].
The subsurface imaging methods using a ground-penetrating radar (GPR) and the frequency domain electromagnetic (FDEM) techniques measure subsurface electrical conductivity and magnetic susceptibility [5]. Combining them with passive tools such as analytical spectral devices (ASD) has made it possible to estimate soil properties in the field rapidly and with high precision [3,6]. The advances in remote-sensing technologies now enable their use to measure SMC, which is responsible for plant and crop growth, and to increase water use effectiveness and efficiency.
The GPR methods currently measure SMC in various ways: multi-offset techniques from reflected waves [7], waveguide inversion of dispersive data [8,9], off-ground measurements [10], full-waveform inversion [11].
These recent studies have developed practical tools to obtain real-time monitoring of changes in spatial SMC up to depths of 3–4 m at low cost. Spatial volumetric soil moisture content (VSMC) measurements can be performed using diffraction hyperbola analysis methods combined with qualitative and quantitative methods.
The aim of this paper was to present the use of diffraction and reflection to estimate the VSMC in the laboratory for various agricultural soils and a wide moisture content range and to establish the first step toward a field-scale GPR application. The final aim of this research is to develop a spatial, non-destructive, non-invasive, cost-effective method based on GPR application for measuring SMC.
Electromagnetic (EM) methods provide an additional advantage in that they enable spatial mapping of the soil structure below the surface. Moisture content measurements using GPR are not limited to a specific location or depth. Spatial measurements performed using GPR systems allow monitoring areas that include thousands of measurement points per hectare, thus providing a picture of the spatial subsurface situation in real time, with precision levels between 1% and 5%. The application of this study’s results will increase the efficiency and effectiveness of the agricultural production and will help decrease water waste in agriculture.

2. Materials and Methods

The chosen research sites represent three main agricultural soil types in Israel, i.e., the northern Negev in the south of Israel, the Jezreel valley in the north of Israel, and the Sharon coastal plain in central Israel (Figure 1).
The soil from the northern Negev area at the Gilat research site (Figure 1) was characterized as light-brown loess (calcic haploxeralf [13]). The soil there is very deep (more than 100 cm depth), with rock fragments mixed with the soil down to a depth of 20 cm [14].
The soil from the Sharon coastal plain area at the Bnei-Dror research site (Figure 1) was characterized as hamra (typic rhodoxeralf [13]). The research site was near the dried-up Dror River area in the eastern Poleg River basin in the Sharon coastal plain in central Israel. The soil was red-brown in color and very deep (more than 100 cm depth), with no rock fragments [14].
The soil from the Jezreel valley area at the Ginegar research site (Figure 1) was characterized as grumusol (typic chromoxerets [13]). The soil was dark brown in color and was very deep (more than 100 cm depth), with no rock fragments [14].
The research was based on GPR measurements and analysis. GPR is a reflection–scattering EM-wave-imaging method; the imaged sections are time-dependent recorded responses to the propagation of EM energy in subsurface materials through a radio-wave frequency range of about 10 MHz to 2 GHz. It has the ability to provide three-dimensional (3D) cross-sectional imaging with high resolution for stratigraphic descriptions of the subsurface in ‘real time’. This makes GPR one of the best tools in shallow geophysical exploration for non-invasive subsurface characterizations [15].
Conventional operation of the GPR system produces reflection profiles in which most returned signals are influenced by the reflections and diffractions from discontinuities in the subsurface. Other types of radio waves can also be created in the system, such as direct waves in the air, refracted waves, and direct ground waves, as expected from ray theory principles and simple geometrical relationships [16,17]. Material permittivity (F/m), ε, material conductivity (S/m), σ , and material permeability (H/m), μ , are the three main components of EM methods and GPR. Dielectric permittivity characterizes the bond charge displacement in the material structure due to the electric field, which causes energy storage in the material. Electrical conductivity σ ˜ characterizes free charge movement (creating an electric current) when an electric field is applied and the resistance to the charge current leads to energy dispersal. Magnetic permeability μ ˜ describes how inherent molecular and atomic magnetic moments respond to the magnetic field. For simple materials, the distortion of inherent magnetic moments leads to energy storage in the material [18].
Permittivity, ε , describes a material’s ability to store and release EM energy in the form of an electric charge; it is related to a capacitor’s ability to store and limit the flow of free charges or the degree of polarization (F/m) demonstrated by a material under the effect of an active electric field. Typically, a dimensionless relative permittivity term ( ε r ) is used [18], also termed a dielectric constant (k). This term can be defined as follows:
  k = ε r = ε ε 0  
Permittivity (dielectric constant) of free air (vacuum) is ε 0 , given as 8.8542 × 10−12 F/m, and ε is the material permittivity. The permittivity of subsurface materials can change significantly, especially in the presence of free or bound water (i.e., it is higher for wetter material, and lower for drier material), and it is usually complex and frequency-dependent, with a real (storage) and imaginary (loss) magnitude. The permittivity value of a material is often simplified to its constant, real part of the complex permittivity at low frequency (static), ignoring the loss expression. It is convenient to calculate an approximation of radar wave velocities and lengths. The relative permittivity of most common subsurface materials changes from ε ave 3 30 .
These parameters—c (~3 × 108 m/s), ω (rad/s), ε , μ , and σ —can be used to estimate velocity (m/s) during experimental measurements and to determine the real velocity characteristics of lossy materials. In most cases, the conductivity constituent is considered to be independent of frequency and, in general, is related only to the ionic conductivity of internal fluids and the surface charge conductivity of clay minerals [19,20,21].
EM waves are influenced by reflections from subsurface elements that are of the approximate size of the transmitted wavelength, and by diffractions created by small discontinuities and objects [15]. In the subsurface medium, diffraction can be detected as a hyperbolic reflection in a temporal section and appears in two cases: when the dominant wavelength, λc, is larger than the source of the diffraction, and when the waves are broken by sharp edges [22,23]. Velocity analysis of hyperbolic reflections in the GPR cross section is a significant factor in measuring volumetric moisture content, because of the relationship between the medium’s wave velocity and the dielectric constant [17]. Calculating velocities from radar cross sections allows obtaining the values of the different dielectric constants in the medium [24]. A practical approximation of the dielectric constant from wave velocity in the soil layers is given by [25]:
  k   = c 2 v 2  
where k is the dielectric constant, c = 0.3 m/ns is the wave velocity in a vacuum, and v is the characteristic radio-wave velocity in the medium.
The volumetric effective relative dielectric constant is affected by soil particles volumes (mineral, liquid, and gas). It can be described approximately as a mean value of all of the dielectric constants of the different soil components [26]:
  k e f f = i = 1 N k i V i  
where keff is the effective relative dielectric constant, ki is the dielectric constant of the specific soil component, Vi is the relative volume of the specific component, and N is the number of soil components.
Soil mineral particles are inorganic materials derived from rocks as a result of weathering, with varying size. The pore space is filled with water and air, with soil moisture and air typically making up about 50% of the soil volume and effective porosity. The difference in the effective dielectric constant of ‘dry’ and ‘wet’ soils is mainly a function of the ratio between the air and water volumes when the total volume is normalized to 1 [26].
  V w = k e f f k s ( 1 n ) n 79  
where Vw is the relative water volume, ks is the soil dielectric constant (values of 3–8), keff is the effective dielectric constant, n is the soil porosity (from 0 to 1).
The laboratory setting (including soil samples from the research sites) consisted of high-resolution GPR measurement experiments in relatively homogeneous soils with a large range of SMC. The soil samples were gathered manually after passing through a 4 mm sifter. GPR measurements included calibration and irrigation experiments in soil boxes, using a shielded 800 MHz antenna operated with a ProEx control unit built by Mala Geoscience AB of Guideline Geo, in a single-offset position. The calibration experiments were performed by static GPR measurement of sections in 21 L soil boxes, while the irrigation experiment was conducted in 270 L and 750 L soil boxes using GPR to monitor the cross sections. The controlled laboratory setting allowed maximum measurement range, resolution, and accuracy.
Soil characteristics such as moisture content (volumetric and gravimetric), porosity, and density, among others, were measured and calculated in the GPR soil laboratory. GPR measurement is volumetric. Therefore, SMC measurements were conducted as such or converted from gravimetric measurements via the soil bulk density, and then correlated to the GPR measurements. A semi analytical balance with accuracy of 0.1 g was used to weigh the soil samples. For ease of use, the term VSMC may be simplified to moisture content.
Soil properties (e.g., density, compaction, moisture content, etc.) differ as a consequence of soil destruction by extraction, sifting, drying, and repacking. The measured soil properties of bulk density and porosity are detailed in Table 1 and Table 2.
The GPR system was calibrated for each soil by performing volumetric measurements at different moisture contents. The measurements were performed at a nominal frequency of 800 MHz in plastic boxes with a volume of 21 L (Figure 2), with a metal plate placed at a known depth at the bottom of the box (Figure 2a) for strong and clear reflection of the GPR signals. The moisture content was gradually raised in steps by 5%. The soil samples were mixed manually before each measurement to rapidly create relative uniformity in the medium (Figure 2b). The antenna was placed above and attached to the soil box surface in a static position for a finite measurement time (Figure 2c). The GPR moisture content measurements were compared to the standard method using 100 mL samples dried in the oven for 48 h at 105 °C [27]. The GPR time-section measurements were conducted in the boxes using a transmitter antenna and a receiving antenna with spatial separation of 14 cm, sampling frequency of 15 GHz, 64 stacks, and time window of 13.8 ns in hamra soil and 19.5 ns in loess and grumusol soils. The GPR measurements were performed shortly after manual mixing.
Laboratory experiments involving irrigation simulations with the geophysical software Reflex-Win (by Sandmeier geophysical research) were performed to predict the SMC measurements in soil boxes with the GPR. The model size was based on the small soil boxes (90 × 50 × 60 cm) with three diffractors at depths of 10, 30, and 50 cm; the distance between the diffractors was 20 cm and 25 cm from the box edges (Figure 3). The EM properties controlling the models were the dielectric constant, electrical conductivity, and electromagnetic wave velocity. The simulation was based on a finite difference (FD) method algorithm for propagating waves in the medium in the time domain. This method allows the simulation of the kinematics (travel time) as well as the dynamics (amplitude, frequency etc.); an explicit method is used, i.e., the wave field is calculated for a certain point of time with the values of the preceding point of time etc. [28]. The section was a single-offset 2D section with parameters of dx = 1 cm (lateral sampling), dt = 0.02 ns (time sampling), and a central frequency of 800 MHz with Kuepper wavelet signal type.
The 2D high-resolution GPR sections for the irrigation experiment were carried out in small soil boxes for loess and hamra soils and in large soil boxes for grumusol. The radar system was moved manually along the section segments. The moisture content was measured at the beginning (air-dried state) and end (saturation state) of the irrigation experiment for each soil type.
The small wood soil box size was 90 × 50 × 60 cm, with a total volume of about 270 L. Inside the boxes, three iron-tube line diffractors were inserted with external diameter of 3 cm, at depths of 10, 30, and 50 cm, placed 20 cm apart and 25 cm from the box edges (Figure 3a). The large metal soil box dimensions were 200 × 60 × 60 cm, with a total volume of about 750 L. Inside the box, there were three iron-bar line diffractors with external diameter of 2 cm at depths of 20, 30, and 40 cm, placed 50 cm apart and 50 cm from the box edges (Figure 3b).
The boxes’ bottom included layers of thin netting, consisting of gravel as the upper level and geotechnical fabric as the upper bottom layer, to prevent soil loss and allow water flow. The boxes were filled to the top with soil (Figure 4). A tap was placed at the box bottom to drain the water after the experiment under saturation conditions. The irrigation system was built using distilled water columns (TREION, Treitel Chemical Engineering) supplying non-saline water through a super-sensitive water meter, specially designed for high accuracy and a wide measuring range. The drip-irrigation system included three lines. In the small boxes, there were five drippers every 15 cm along each line (the same as the distance between the lines), with built-in drippers. The flow rate in this system reached 30 L/h, whereas, in the large box with 36 drippers, it reached 72 L/h.
Much of the general GPR processing was adapted from high-level seismic data analysis [29]. Basic processing can be carried out in a real-time display for on-site interpretation and, therefore, was adopted for the laboratory irrigation experiment, highlighting the diffractions hyperbolic pattern. More complex processing (e.g., migration, de-convolution, etc.) was tested but discarded, since the research focused on diffraction hyperbola analysis close to real time in the field.
The applied filters in the data processing were: (i) de-wow: Ideally, the baseline amplitude of each trace should be zero. However, because of the dynamic range limitations of the instrumentation and the EM induction effects in the conductive ground, there were low-frequency component variations in the baseline amplitude, including bias or constant shift [30]. The temporal filter removed the low-frequency components, and the running mean value of each trace was subtracted from the signal. (ii) Time zero: The first arrival time from the ground can be shifted by thermal or electronic drift, cable-length differences, or changes in antenna height. It needs to be leveled to a common datum (e.g., max phase correction). The subtraction of the first arrival time was also needed for accurate evaluation of later events. (iii) Time gain: Wave propagation in the medium is accompanied by geometrical spreading losses (approximately 1/r2) and signal attenuation. To equalize the signal return from all depths, a time-dependent gain function was applied. An energy-decay function, also known as spherical and exponential compensation (SEC), automatically amplified the signal amplitude to compensate for the energy loss by geometric spreading. The automatic gain control (AGC) function enhanced each trace by creating a temporally equally distributed signal amplitude within a predefined time window. This is a nonlinear process, so other filters will operate differently before and after time gain. It also changes the relative amplitudes and phase relationship [31].

3. Results

3.1. Calibration Experiments

The GPR signal reflection, using an 800 MHz antenna, from the metal bottom in the plastic 21 L soil boxes used for calibration, was identified on a single representative trace and by generic static section analysis (Figure 5). The image in Figure 5b shows a static GPR section from the grumusol soil calibration experiment in an air-dried condition (without added water), processed only with de-wow filtering with a time range window length of 1.25 ns. The horizontal scale is the measured time, while the vertical scale is the two-way travel time of the signal through the medium and back to the receiver. The scale on the right is the voltage amplitude of the signals. The left image is a wiggle trace window presentation of the color-mode profile, chosen from the section to accurately calculate the differences in travel time: t0 is time zero (upper orange dotted line), and t1 is the two-way travel-time (TWT) of the reflection wave from the metal plate at the bottom of the calibration soil box (lower orange dotted line).
A slight but visible change between each section (Figure 6) can be noted. The velocity decreased with increasing moisture content, such that the returning signal from the bottom was delayed because of slower propagation of the EM signal in the medium. These changes were analyzed and used to calculate the GPR-derived moisture content parameters (Table 3). The travel-time analysis was performed at a very high accuracy of ±0.065 ns, derived from the defined system frequency sampling, or by the section time window and the number of samples along the measuring trace. The average error of EM velocities was ±0.003 m/ns, leading to an effective dielectric constant average error of ±0.64, resulting in an average moisture content error of ±0.81% in the grumusol calibration experiment.
Graphs of the calibration experiment results (Figure 7) were constructed by comparing each moisture content increment in the soil boxes. The differences between the moisture content values derived from the GPR and the standard method (volumetric oven-drying) measurements are shown. GPR sections of loess soil included 13 steps with increased water levels, demonstrating an EM velocity range of 0.141 ± 0.010–0.045 ± 0.001 m/ns. The calculated effective dielectric constant ranged from 4.50 ± 0.66 to 45.13 ± 2.09. The moisture content results (Figure 7a), which included 13 iteration steps, ranged from 2.4 ± 0.86 to 53.8 ± 2.7%, with an average step of 4.0%, linear regression y = 3.8181x + 6.4554, and R2 = 0.9223. Standard sample measurements showed a range of 6.1–57.9%, with the same average step of 4.0%, linear regression y = 3.7894x + 11.232, and R2 = 0.9699. The linear differential calibration was y = −0.0287x + 4.7766.
In the hamra soil calibration experiment (Figure 7b), the velocity ranged from 0.129 ± 0.009 to 0.049 ± 0.001 m/ns. The calculated effective dielectric constant range was 5.39 ± 0.72–37.93 ± 1.91, leading to a moisture content range of 3.1 ± 0.9–44.3 ± 2.4%, with an average step of 4.1%, linear regression y = 4.6126x + 2.3558, and R2 = 0.9876. Standard sample measurements ranged from 3.2 to 46.8%, with an average step of 4.4%, linear regression y = 4.4554x + 4.9849, and R2 = 0.9897. The linear differential calibration was y = −0.1572x + 2.6291. The data were obtained in 10 iteration steps with a GPR travel time measurement error of 0.1 ns.
The grumusol soil calibration experiment (Figure 7c) showed a velocity range of 0.010 ± 0.004–0.080 ± 0.002 m/ns, a calculated effective dielectric constant range of 7.43 ± 0.55–13.96 ± 0.75, and moisture content results ranging from 7.5 ± 0.7 to 15.7 ± 0.9%, with an average step of 1.6%, linear regression y = 1.9187x + 6.7781, and R2 = 0.9173. Standard sample moisture content measurements ranged from 14.4 to 31.6%, with an average step of 3.4%, linear regression y = 4.1333x + 12.653, and R2 = 0.9526. The moisture content results included five iteration steps, and the analysis of the GPR data gave a signal travel-time measurement error of ±0.1 ns.
The major advantage of the GPR method are its spatial remote-sensing (non-invasive) capabilities, whereas the standard method requires drilling at local points. The correlation between the different methods is not simple. The standard method samples a local 100 mL volume, whereas the GPR senses larger volumes. There are some moisture content differences throughout a soil because it is never completely homogeneous or isotropic as assumed. The measurement is affected by changes in volume, location, direction, and time. As a result, a comparison of the two methods will show deviations in the results, and a correlation ratio is therefore required for their comparison. The results should not be interpreted only according to the values but also according to their trends and assessed for relative similarity or contrast. Calibration experiments are needed to create a connection between the two methods, especially when an absolute quantitative value is desired rather than just a relative or qualitative one.

3.2. Irrigation Experiments

3.2.1. Simulation

A simulation process was used to improve the predictions from the soil box irrigation experiments, for both small and large boxes. Twenty-one simulations through all SMC ranges, at 5% water content increments, were conducted with the loess soil parameters in the small box model (Figure 8). The dielectric constant was composed of three phases: (i) the solid soil part—sand, silt, and clay, (ii) the liquid part—water, and (iii) the gas part—air. In the experiment, the porosity was set to 50%. The dielectric constant of the solid part was 3.9. The SMC was increased theoretically from 0% to 100% of the pore space by 5% steps, while the air was reduced equally from 100% to 0% of the pore space. A theoretical 0–100% moisture content had a bulk dielectric constant range of 2.45–41.95 and a wave-velocity range of 0.192–0.046 m/ns.
The simulation experiment showed visible quality and quantity differences between the moisture content steps. The diffractors and the reflector were in a constant location in the box (10, 30, 50, and 60 cm), but they appeared deeper in each section because of the larger signal travel time caused by the increased SMC. The increased moisture content was accompanied by an increased dielectric constant and a decreased velocity, shown directly in the GPR sections. A higher water amount in the pore space instead of air decreased the EM velocity because the velocity in water (about 0.033 m/ns) is approximately 1/9 less than in air (0.3 m/ns).
The hyperbola shapes were created by the three diffractors, the shape was defined by the signal travel time, medium parameters, reflectivity, target size, etc., and the strong horizontal reflection was created by the soil box bottom. The inclined reflections were created by the diffractors boundary reflections from the box edges.

3.2.2. Experiments

The irrigation experiments in the small soil boxes showed changes in the soil average wave velocity with increasing SMC. GPR sections in grumusol, hamra, and loess soil boxes are shown in Figure 9, with the locations of the diffractors and reflectors, after processing of time-zero correction (moving start time), de-wow with time range window length of 1–1.5 ns, and gain by AGC, with a window length of 5–15 ns and a multiplication factor of 5–10, or energy decay curve, usually with a multiplication factor of 1; the amplitude scale was adjusted to fit each section. The measured average EM velocity from the soil surface to the diffractors and reflectors was 0.112 ± 0.002 m/ns, 0.156 ± 0.009 m/ns, and 0.139 ± 0.007 m/ns in the grumusol, hamra, and loess soil boxes, respectively. The bulk dielectric constant range calculated from those velocities was 7.15 ± 0.26, 3.71 ± 0.42, and 4.63 ± 0.46, respectively, with moisture content of 7.1 ± 0.3%, 1.0 ± 0.5%, and 2.56 ± 0.58%, respectively. The diffractors and reflectors (from the box bottom) marked in yellow (ellipses and lines) in the figure, were clearly identified; those marked in green were not, and their approximate expected location is indicated. The elliptical shape of the diffractors in the GPR section, which in reality are spherical, was due to the different distance scales of the horizontal and vertical axes.
In the small loess soil box, changes in EM velocity due to increasing water content were observed and calculated from the GPR sections, with some shown in Figure 10. The added distilled water in those sections was: 0 (air-dried), 25.00 ± 0.63 L, 55.00 ± 1.38 L, and 80.00 ± 2.00 L. The average measured EM velocities were 0.139 ± 0.007 m/ns (air-dried), 0.093 ± 0.003 m/ns, 0.074 ± 0.002 m/ns, and 0.064 ± 0.001 m/ns, respectively. The bulk dielectric constants calculated from those velocities were: 4.63 ± 0.46, 10.33 ± 0.68, 16.31 ± 0.86, and 22.11 ± 1.00, respectively, and the moisture contents were 2.6 ± 0.6%, 9.8 ± 0.9%, 17.3 ± 1.1%, and 24.7 ± 1.3%, respectively.
Some of the GPR sections from the small hamra soil box are shown in Figure 11, with zero (air-dried), 20.25 ± 0.51, and 54.00 ± 1.35 L of distilled water. The measured average EM velocity values were 0.156 ± 0.009 m/ns, 0.126 ± 0.013 m/ns, and 0.089 ± 0.003 m/ns, respectively. The bulk dielectric constants calculated from those velocities were 3.71 ± 0.42, 5.69 ± 1.20 and 11.27 ± 0.74, respectively, and the calculated moisture contents were 1.0 ± 0.5%, 3.5 ± 1.5% and 10.5 ± 0.9%, respectively.

4. Discussion

The changes in EM properties from the GPR measurement calibration experiments were correlated to the moisture content measurements (oven dry) in each tested soil. Increasing the moisture content directly affected the relative dielectric constant and the wave propagation velocity, influencing the GPR signals.
The calibration experiment presented the GPR subsurface imaging in loess and hamra soils for a wide range of moisture contents. GPR measurement of the grumusol soil type showed partial capabilities, with reliable results at lower moisture contents, starting from 14.4% up to about 31.6% (see error analysis above (in Section 3)), and less reliable results at the upper moisture contents. These results are meaningful for agricultural applications, especially in arid and semi-arid regions, where natural water is scarce, and farmers usually irrigate to field capacity (additional irrigation above field capacity will be lost because of gravitational forces).
The calibration experiment results showed that the GPR method slightly underestimated the moisture content in hamra soil by an average 1.9% with 2.3 SD (standard deviation) and in loess soil by 4.6% with 3.4 SD compared to the standard method. Grumusol calibration results showed an average underestimation of moisture content of 10.5% with 3.4 SD, which is approximately half that of the standard method.
The simulation experiment showed increasing VSMC values modelled in the small loess soil box, affecting the bulk dielectric constant and the average EM wave velocity. The soil box’s diffractors and bottom reflector were clearly observed, although the upper diffraction and reflection waves from the first and second diffractors interfered with those from the lowest diffractor and the bottom reflector. Furthermore, the deepest diffractor was quite close to the bottom, so we could expect some difficulty in identifying it, which disturbed the analysis at a greater depth (>50 cm). In the large metal soil box, the diffractions were separated, and the anomalies did not overlap significantly, so each diffraction had enough space to be well represented in the measured GPR section.
The irrigation experiments enabled measurements of field-size diffractors, which can be used for velocity analysis and VSMC calculation. The experiments began with small soil boxes that provided partial results in both loess and hamra soils. The two shallower diffractors were observed clearly and were fit for analysis in most of the VSMC range, but the deep diffractor and the soil box bottom were not. The clarity of their identification decreased considerably as a consequence of wave interference of the upper diffractors and reflections of the soil box edges and bottom. Since the wooden boxes were found to be too small for the irrigation experiment setting, the larger metal frame soil boxes were used.
In the large soil box, all three diffractors and the reflection from the bottom were identifiable thanks to the spacing effect, but, unfortunately, only the first measurement was successful. After the first irrigation, the GPR imaging was blurred, and no velocity analysis could be performed. This lack of clarity was probably due to the high clay fraction (50 ± 5%), with a strong effect on the electrical conductivity value and strong signal attenuation, which was also influenced by the soil structure, grains, and moisture content, resulting in a very low signal-to-noise ratio. Moreover, the soil was disturbed, as it was collected from the field, sifted, and ground and thus it was no longer in its unified natural condition, with respect to grain behavior and internal forces. Some of these problems can be fixed over time; for example, when the small loess soil box was filled, the GPR signal results were not clear, so the diffractors were replaced with larger ones, and some more time was given for the soil to restore itself. After these steps, the GPR data were clear.
For some of the GPR measurements, a lack of clarity resulted from anomalies related to soil moisture. The irrigation water in the soil produced strong reflections and diffraction waves, especially at the beginning of the irrigation experiment when there was a large difference in the effective dielectric constants between the upper layer with high moisture content and the lower layer with relatively low moisture content. Thereafter, the water spread more equally throughout the soil boxes, and the change in the effective dielectric constants was less dramatic, so the diffractors were no longer masked.
This technique could be applied directly on undisturbed soils in the field, even though the distance to the anomaly is unknown. The diffractions hyperbola shape in the GPR profiles can be analyzed by dividing it into many point-velocity hyperbolas or interpolated to 2D velocity profiles, converted to effective dielectric constant values, and used in a soil mixture model to obtain the subsurface spatial soil water content, relying on the calibration procedure described above (Section 2) combined with a few in situ soil samples measurements, mainly of soil texture.
The results of the study showed that the GPR can measure the spatial water content in various agricultural soils. The clarity of the profiles enabled the analysis and interpretation of the physical properties with high accuracy. The non-invasive measurements of a soil’s EM parameters can be relevant for other applications, such as subsurface targets detection, localization, and classification (e.g., explosive hazards) by the analysis of the hyperbola diffraction anomaly. GPR monitoring surveys of repeating measurements (as conducted in the calibration and irrigation experiments described above in Section 2) can also show the subsurface changes over time caused by various factors such as irrigation, cultivation, crop type, seasons, etc. These changes can be analyzed to increase the amount of information on the soil physical parameters in a field.

5. Conclusions

Laboratory GPR studies in soil boxes showed that it is possible to measure the moisture content from a known depth using diffraction hyperbola, at gradually increasing soil moisture content (from dry to saturation state), with high accuracy. The GPR profiles were analyzed by signal parameters (e.g., arrival time) to extract EM velocities, then converted to effective dielectric constant values and used in a soil mixture model that combined a few soil samples, mainly with different soil texture, to obtain the subsurface spatial soil moisture content. These results were correlated to those from soil samples subjected to the standard measurement method of oven-drying. Good agreement was found between the measurement methods (volumetric oven-drying and GPR).
High spatial resolution (lateral and vertical) is particularly important for this application, where the required penetration depth is low (root zone); small subsurface discontinuities must be clearly identifiable to obtain information on the EM properties of the soil, so that the soil moisture content can be calculated with high accuracy for spatial agricultural practices.
The GPR analysis showed high applicability for a large range of moisture contents in soils with various textures. Several elements can be addressed to achieve better results: (i) higher antenna bandwidth for higher signal resolution; (ii) higher measurement stack number to increase the signal-to-noise ratio; (iii) laboratory measurements in larger soil boxes to reduce the edge effects.
The average moisture content measurements were conducted with an accuracy of ±1.0–2.0%. GPR imaging in the field will be based on these laboratory results and others and will enable computing the average moisture content with high accuracy.
This method relied on a calibration study in a soil box (21 L volume) to achieve the presented accuracy and to improve our understanding of the physical relations between relative permittivity and moisture content. Some knowledge is required to build a simple and reliable volume mixture model that will describe the relationship between the EM characteristics of the medium and its moisture content. To achieve the accuracy levels presented in this article, a calibration process is required, preferably using a small soil box for a rapid and accurate technique. The intensive development in recent years of GPR instrumentation now allows close-to-real-time, rapid, and highly accurate geophysical imaging and extraction of a soil’s physical properties, among which soil moisture content.
The results of this study provide the basis for GPR measurements of moisture content in the field, using the calibration and irrigation experiment models from this research.
The GPR measurement results must be correlated to those from standard methods (e.g., oven-drying) as in the calibration experiment (with varying soil texture and moisture content) to increase result accuracy and reliability. Further studies are needed to obtain more information on the relations between the dielectric constant and the moisture content in various agricultural soils and on the correlation between the GPR method and other physical methods, such as FDEM, spectroscopy, synthetic aperture radar, and time-domain reflectometry, among others.
There are many regions (arid and semi-arid) with water-related problems, such as salinized soils and areas that use mostly retrieved water in unwashed soils. The farmers in those areas need to optimize irrigation and production. Because of the need to use water efficiently, this study is mostly significant for soils with a low moisture content, from air-dry soils to those at field capacity, as these face the presented challenges.

Author Contributions

O.S., N.G., and U.B., contributed to the design and implementation of the research, O.S., U.B., and M.R. contributed to the analysis, processing and interpretation of the data, O.S. and N.G. wrote the article with input from all authors.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Huisman, J.A.; Hubbard, S.S.; Redman, J.D.; Annan, A.P. Measuring Soil Water Content with Ground Penetrating Radar: A Review. Vadose Zone J. 2003, 2, 476–491. [Google Scholar] [CrossRef]
  2. Topp, G.C.; Parkin, G.W.; Ferre, T.P.A. Soil water content. In Soil Sampling and Methods of Analysis, 2nd ed.; Carter, M.R., Gregorich, E.G., Eds.; Canadian Society of Soil Science, CRC Press: Boca Raton, FL, USA, 2008; pp. 939–962. [Google Scholar]
  3. Ben-Dor, E.; Chabrillat, S.; Demattê, J.A.M.; Taylor, G.R.; Hill, J.; Whiting, M.L.; Sommer, S. Using Imaging Spectroscopy to study soil properties. Remote Sens. Environ. 2009, 113 (Suppl. 1), S38–S55. [Google Scholar] [CrossRef]
  4. Goldshleger, N.; Livne, I.; Chudnovsky, A.; Ben-Dor, E. New Results in Integrating Passive and Active Remote Sensing Methods to Assess Soil Salinity: A Case Study from Jezre’el Valley, Israel. Soil Sci. 2012, 177, 392–401. [Google Scholar] [CrossRef]
  5. Goldshleger, N.; Basson, U. Utilization of Ground-Penetrating Radar and Frequency Domain Electromagnetic for Investigation of Sewage Leaks. In Environmental Applications of Remote Sensing; Marghany, M., Ed.; InTech: London, UK, 2016. [Google Scholar] [Green Version]
  6. Goldshleger, N.; Ben-Dor, E.; Benyamini, Y.; Agassi, M. Soil Reflectance as a Tool for Assessing Physical Crust Arrangement of Four Typical Soils in Israel. Soil Sci. 2004, 169, 677–687. [Google Scholar] [CrossRef]
  7. Gerhards, H.; Wollschläger, U.; Yu, Q.; Schiwek, P.; Pan, X.; Roth, K. Continuous and simultaneous measurement of reflector depth and average soil-water content with multichannel ground-penetrating radar. Geophysics 2008, 73, J15–J23. [Google Scholar] [CrossRef]
  8. Mangel, A.R.; Moysey, S.M.; van der Kruk, J. Resolving precipitation induced water content profiles by inversion of dispersive GPR data: A numerical study. J. Hydrol. 2015, 525, 496–505. [Google Scholar] [CrossRef]
  9. Mangel, A.R.; Moysey, S.M.; van der Kruk, J. Resolving infiltration induced water content profiles by inversion of dispersive ground-penetrating radar data. Vadose Zone J. 2017, 16. [Google Scholar] [CrossRef]
  10. Serbin, G.; Or, D. Ground-penetrating radar measurement of crop and surface water content dynamics. Remote Sens. Environ. 2005, 96, 119–134. [Google Scholar] [CrossRef]
  11. Meles, G.A.; Van der Kruk, J.; Greenhalgh, S.A.; Ernst, J.R.; Maurer, H.; Green, A.G. A new vector waveform inversion algorithm for simultaneous updating of conductivity and permittivity parameters from combination crosshole/borehole-to-surface GPR data. IEEE Trans. Geosci. Remote Sens. 2010, 48, 3391–3407. [Google Scholar] [CrossRef]
  12. Dan, Y.; Raz, Z.; Yaalon, D.H.; Koyumdjisky, H. Soil Map of Israel; Ministry of Agriculture: Bet Dagan, Israel, 1975.
  13. Soil Survey Staff. Keys to Soil Taxonomy, 12th ed.; USDA-Natural Resources Conservation Service: Washington, DC, USA, 2014.
  14. Soil Conservation Department, Ministry of Agriculture. Israel, ArcGIS–Soil Survey Map of Israel (1:20,000), Esri Online Services. 1952. Available online: https://www.arcgis.com/apps/MapJournal/index.html?appid=ddb7490e8d214e22a7f1ab8196225521 (accessed on 7 August 2018).
  15. Annan, A.P. Ground-Penetrating Radar. In Near-Surface Geophysics; Society of Exploration Geophysicists: Tulsa, OK, USA, 2005; pp. 357–438. [Google Scholar]
  16. Davis, J.L.; Annan, A.P. High-Resolution Sounding Using Ground-Probing Radar. Geoscience 1986, 13, 205–208. [Google Scholar]
  17. Davis, J.L.; Annan, A.P. Ground-Penetrating Radar for High-Resolution Mapping of Soil and Rock Stratigraphy. Geophys. Prospect. 1989, 37, 531–551. [Google Scholar] [CrossRef]
  18. Cassidy, N.J. Electrical and Magnetic Properties of Rocks, Soils and Fluids. In Ground Penetrating Radar Theory and Applications; Jol, H.M., Ed.; Elsevier: Amsterdam, The Netherlands, 2009; pp. 41–72. [Google Scholar]
  19. Cassidy, N.J. Evaluating LNAPL contamination using GPR signal attenuation analysis and dielectric property measurements: Practical implications for hydrological studies. J. Contam. Hydrol. 2007, 94, 49–75. [Google Scholar] [CrossRef] [PubMed]
  20. Hargens, C.W. Electromagnetic fields and waves. J. Frankl. Inst. 1961, 272, 236–237. [Google Scholar] [CrossRef]
  21. Turner, G.; Siggins, A.F. Constant Q attenuation of subsurface radar pulses. Geophysics 1994, 59, 1192–1200. [Google Scholar] [CrossRef]
  22. Landa, E.; Shtivelman, V.; Gelchinsky, V. A method for Detection of Diffracted Waves on Common Offset Sections. Geophys. Prospect. 1987, 35, 359–374. [Google Scholar] [CrossRef]
  23. Landa, E.; Keydar, S. Seismic Monitoring of Diffraction Images for Detection of Local Heterogeneities. Geophysics 1998, 63, 1093–1100. [Google Scholar] [CrossRef]
  24. Topp, G.C.; Davis, J.L.; Annan, A.P. Electromagnetic Determination of Soil Water Content: Measurements in Coaxial Transmission Lines. Water Resour. Res. 1980, 16, 574–582. [Google Scholar] [CrossRef]
  25. Reynolds, J.M. An Introduction to Applied and Environmental Geophysics, 2nd ed.; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2011. [Google Scholar]
  26. Basson, U. Mapping of Moisture Content and Structure of Unsaturated Sand Layers with Ground Penetrating Radar. Master’s Thesis, Tel Aviv University, Tel Aviv, Israel, 1992. [Google Scholar]
  27. Topp, G.C.; Ferré, P.A. Water content. In Methods of Soil Analysis, Part—Physical Methods; Dane, J.H., Topp, G.C., Eds.; Soil Science Society of America: Madison, WI, USA, 2002; Chapter 3.1; pp. 417–545. [Google Scholar]
  28. Kelly, K.R.; Ward, R.W.; Treitel, S.; Alford, R.M. Synthetic seismograms: A finite difference approach. Geophysics 1976, 41, 2–27. [Google Scholar] [CrossRef]
  29. Yilmaz, Ö. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data; Doherty, S.M., Ed.; Society of Exploration Geophysicists: Tulsa, OK, USA, 1987. [Google Scholar]
  30. Everett, M.E. Near-Surface Applied Geophysics; Cambridge University Press: Cambridge, UK, 2013. [Google Scholar]
  31. Annan, A.P. Practical Processing of GPR Data. In Proceedings of the Second Government Workshop on Ground Penetrating Radar; Sensors & Software Inc.: Mississauga, ON, Canada, 1999. [Google Scholar]
Figure 1. Soil groups map of Israel [12] with the research sites marked with red circles (Gilat, Bnei-Dror, and Ginegar). Soil legend: E—hamra soils, N—loessial arid brown soils, H—grumusol soils.
Figure 1. Soil groups map of Israel [12] with the research sites marked with red circles (Gilat, Bnei-Dror, and Ginegar). Soil legend: E—hamra soils, N—loessial arid brown soils, H—grumusol soils.
Remotesensing 10 01667 g001
Figure 2. Ground-penetrating radar (GPR) laboratory calibration experiments in 21 L plastic soil boxes. (a) Metal plate at the bottom; (b) loess soil, mixed manually (15% volumetric soil moisture content (VSMC) added in steps of 5%); (c) 800 MHz antenna in static measurement position.
Figure 2. Ground-penetrating radar (GPR) laboratory calibration experiments in 21 L plastic soil boxes. (a) Metal plate at the bottom; (b) loess soil, mixed manually (15% volumetric soil moisture content (VSMC) added in steps of 5%); (c) 800 MHz antenna in static measurement position.
Remotesensing 10 01667 g002
Figure 3. Side sketch of the soil boxes in the irrigation experiments. (a) Small soil box; (b) large soil box. See text for details.
Figure 3. Side sketch of the soil boxes in the irrigation experiments. (a) Small soil box; (b) large soil box. See text for details.
Remotesensing 10 01667 g003
Figure 4. Small and large soil boxes for GPR laboratory irrigation experiments. (a) Empty small soil box; (b) small soil box filled with loess soil; (c) large soil box bottom, layer of netting; (d) layer of gravel; (e) layer of geotechnical fabric; (f) large soil box filled with grumusol soil.
Figure 4. Small and large soil boxes for GPR laboratory irrigation experiments. (a) Empty small soil box; (b) small soil box filled with loess soil; (c) large soil box bottom, layer of netting; (d) layer of gravel; (e) layer of geotechnical fabric; (f) large soil box filled with grumusol soil.
Remotesensing 10 01667 g004
Figure 5. GPR section of calibration of the soil box experiment; the amplitude scale is on the right. (a) Specific wiggle trace; (b) GPR profile of static measurement. The upper orange line is time zero (t0), and the lower orange line is the reflection from the metal bottom (t1).
Figure 5. GPR section of calibration of the soil box experiment; the amplitude scale is on the right. (a) Specific wiggle trace; (b) GPR profile of static measurement. The upper orange line is time zero (t0), and the lower orange line is the reflection from the metal bottom (t1).
Remotesensing 10 01667 g005
Figure 6. Grumusol soil GPR calibration static measurement in air-dried condition and in the presence of up to 4.2 L of distilled water added to the soil box. (a) Air-dried condition; (b) 1.05 L added water; (c) 2.1 L added water; (d) 3.15 L added water; (e) 4.2 L added water. The vertical axes are the GPR measurement times. The sections have different amplitude values, which decrease slightly with each step.
Figure 6. Grumusol soil GPR calibration static measurement in air-dried condition and in the presence of up to 4.2 L of distilled water added to the soil box. (a) Air-dried condition; (b) 1.05 L added water; (c) 2.1 L added water; (d) 3.15 L added water; (e) 4.2 L added water. The vertical axes are the GPR measurement times. The sections have different amplitude values, which decrease slightly with each step.
Remotesensing 10 01667 g006
Figure 7. Soil calibration GPR measurements compared to a standard (oven-dried) sample with linear trend lines (by regression). (a) Loess; (b) hamra; (c) grumusol.
Figure 7. Soil calibration GPR measurements compared to a standard (oven-dried) sample with linear trend lines (by regression). (a) Loess; (b) hamra; (c) grumusol.
Remotesensing 10 01667 g007
Figure 8. Simulation experiment in a small loess soil box with three diffractors and a bottom reflector. The modeled moisture content (MC) range was between the theoretical minimum and maximum conditions, from dry (0%) to saturation (100%) by 5% steps: (a) 0% MC; (b) 5% MC; (c) 10% MC; (d) 15% MC; (e) 20% MC; (f) 25% MC; (g) 30% MC. The effective scalar dielectric constants and velocity values of the sections are presented.
Figure 8. Simulation experiment in a small loess soil box with three diffractors and a bottom reflector. The modeled moisture content (MC) range was between the theoretical minimum and maximum conditions, from dry (0%) to saturation (100%) by 5% steps: (a) 0% MC; (b) 5% MC; (c) 10% MC; (d) 15% MC; (e) 20% MC; (f) 25% MC; (g) 30% MC. The effective scalar dielectric constants and velocity values of the sections are presented.
Remotesensing 10 01667 g008
Figure 9. GPR sections in air-dried soil boxes from the laboratory irrigation experiments with a central frequency of 800 MHz. (a) Grumusol; (b) hamra; (c) loess.
Figure 9. GPR sections in air-dried soil boxes from the laboratory irrigation experiments with a central frequency of 800 MHz. (a) Grumusol; (b) hamra; (c) loess.
Remotesensing 10 01667 g009
Figure 10. Representative GPR sections in the small loess soil box with different amounts of water added to the soil box. (a) No added water; (b) 25 L added water; (c) 55 L added water; (d) 80 L added water; the calculated moisture contents were 2.6 ± 0.6%, 9.8 ± 0.9%, 17.3 ± 1.1%, and 24.7 ± 1.3, respectively. The yellow markers indicate the clearly identified diffractors, and the green markers indicate the approximate location of the targets that are expected to have provided diffractions and reflections.
Figure 10. Representative GPR sections in the small loess soil box with different amounts of water added to the soil box. (a) No added water; (b) 25 L added water; (c) 55 L added water; (d) 80 L added water; the calculated moisture contents were 2.6 ± 0.6%, 9.8 ± 0.9%, 17.3 ± 1.1%, and 24.7 ± 1.3, respectively. The yellow markers indicate the clearly identified diffractors, and the green markers indicate the approximate location of the targets that are expected to have provided diffractions and reflections.
Remotesensing 10 01667 g010
Figure 11. GPR sections in the hamra small soil box. The moisture contents were: (a) 1.0 ± 0.5% (air-dried); (b) 3.5 ± 1.5%; (c) 10.5 ± 0.9% (saturation).
Figure 11. GPR sections in the hamra small soil box. The moisture contents were: (a) 1.0 ± 0.5% (air-dried); (b) 3.5 ± 1.5%; (c) 10.5 ± 0.9% (saturation).
Remotesensing 10 01667 g011
Table 1. Soil bulk density (kg/m3) measured in laboratory experiments for all soils.
Table 1. Soil bulk density (kg/m3) measured in laboratory experiments for all soils.
ExperimentLoessHamraGrumusol
Calibration1.331.411.13
Irrigation1.21.691.22
Table 2. Soil porosity (m3/m3) measured in laboratory experiments for all soils.
Table 2. Soil porosity (m3/m3) measured in laboratory experiments for all soils.
ExperimentLoessHamraGrumusol
Calibration0.4970.4670.575
Irrigation0.550.3630.539
Table 3. Grumusol soil calibration: static calculations of the GPR-derived VSMC parameters.
Table 3. Grumusol soil calibration: static calculations of the GPR-derived VSMC parameters.
Added Water (L)t0 (ns)t1 (ns)dt (ns)v (m/ns)kVSMC (%)
0.002.3405.8503.5100.1107.437.5
1.052.3406.1103.7700.1028.578.9
2.102.3406.1753.8350.1018.879.3
3.152.2106.5654.3550.08911.4412.6
4.202.2757.0854.8100.08013.9615.7

Share and Cite

MDPI and ACS Style

Shamir, O.; Goldshleger, N.; Basson, U.; Reshef, M. Laboratory Measurements of Subsurface Spatial Moisture Content by Ground-Penetrating Radar (GPR) Diffraction and Reflection Imaging of Agricultural Soils. Remote Sens. 2018, 10, 1667. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10101667

AMA Style

Shamir O, Goldshleger N, Basson U, Reshef M. Laboratory Measurements of Subsurface Spatial Moisture Content by Ground-Penetrating Radar (GPR) Diffraction and Reflection Imaging of Agricultural Soils. Remote Sensing. 2018; 10(10):1667. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10101667

Chicago/Turabian Style

Shamir, Omer, Naftaly Goldshleger, Uri Basson, and Moshe Reshef. 2018. "Laboratory Measurements of Subsurface Spatial Moisture Content by Ground-Penetrating Radar (GPR) Diffraction and Reflection Imaging of Agricultural Soils" Remote Sensing 10, no. 10: 1667. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10101667

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop