Electromagnetic induction mapping at varied soil moisture reveals field-scale soil textural patterns and gravel lenses

Abstract

Knowledge of the spatial distribution of soil textural properties is important for determining soil moisture storage and soil hydraulic transport properties. Capturing field heterogeneity without exhaustive sampling and costly sample analysis is difficult. Our objective was to employ electromagnetic induction (EMI) mapping in low apparent electrical conductivity (EC_a) soils at varying soil water contents to capture time invariant properties such as soil texture. Georeferenced EC_a measurements were taken using a ground conductivity meter on six different days where volumetric water content (θ_v) varied from 0.11 to 0.23. The 50 m × 50 m field included a subsurface gravelly patch in an otherwise homogeneous silt-loam alluvial soil. Ordinary block kriging predicted EC_a at unsampled areas to produce 1-m resolution maps. Temporal stability analysis was used to divide the field into three distinct EC_a regions. Subsequent ground-truthing confirmed the lowest conductivity region correlated with coarse textured soil parent materials associated with a former high-energy alluvial depositional area. Combining maps using temporal stability analysis gives the clearest image of the textural difference. These maps could be informative for modeling, experimental design, sensor placement and targeted zone management strategies in soil science, ecology, hydrology, and agricultural applications.

Keywords

soil electrical conductivity / soil texture mapping / temporal stability analysis

1 Introduction

Quantitative determination of the spatial properties of soils at the field scale remains a research challenge. The development of precision agriculture with the targeted application of fertilizer or irrigation could be further enhanced with improved, quantitative, soil textural data at the field scale. Geophysical methods are gaining more acceptance as a way of obtaining spatially distributed data that can be correlated with soil properties^[1^]. Although techniques such as ground-penetrating radar (GPR) and electromagnetic induction (EMI) are available to generate valuable spatial information, determining efficient ways to extract useful information remains a challenge^[2,3^]. Mapping protocols and calibration procedures have been developed for soil salinity surveys, which remains an active area of research^[4^]. However, in non-saline agricultural soils we are often interested in the spatial delineation of soils with different characteristics, usually textural differences, with the minimum amount of soil sampling and laboratory calibration, which adds expense to surveying.

Identifying time-invariant static properties such as texture is important as they have a direct influence on the quantity and distribution of dynamic soil properties such as soil water content. Soil water status is critical to plant growth, crop quality, chemical fate and transport, and microbial processes^[5^]. Soil structure and texture are important properties controlling the hydraulic conductivity and infiltration capacity of a soil system^[6^]. These two properties in turn determine how much water is stored in a soil system as well as how much water is diverted from the system due to runoff. Soil water content is an important factor for many agricultural practices such as the timing of tillage or fertilizer application and recent research effort has focused on this^[7,8^]. More importantly in dry climates, knowledge of the soil water content can reduce water waste and increase efficient use of this precious resource by targeting irrigation more effectively^[5^]. Thus, knowledge of soil textural properties in space can help to infer soil hydrological behavior.

Traditional soil mapping strategies tend to be subjective and rely on the expertise of the soil surveyor^[9^]. Discrete point locations in the form of a pedon are chosen to sample soil physical and chemical properties in order to determine soil properties for a given soil mapping unit. Thus, soil regions are dependent on the original scheme of sampling which in turn is dependent on vegetation coverage, topography, and geological features. While this might be sufficient to classify large areas, it is biased toward grouping and overlooks small scale soil variability, which is important for management and understanding hydrological and ecological patterns. Our understanding of the distribution of subsurface physical properties is also limited by the sparse sampling plan and consequent interpolation method. Thus defining subsurface physical property boundaries remains largely subjective using traditional sampling and survey methods.

Geophysical methods are frequently used in characterizing the subsurface by measuring a surrogate property that is related to underlying physical properties. Some popular geophysical methods used for characterizing the subsurface are: GPR, electrical resistivity tomography/imaging (ERT/I) and EMI^[10–12^]. A GPR system transmits a high frequency (MHz–GHz) electromagnetic radiation into the subsurface and receives a reflected signal that has been transformed by the electrical properties of the subsurface material. The reflected wave is dependent upon the dielectric permittivity of various subsurface layers. The application of GPR for spatial soil mapping is restricted by the attenuation of the transmission signal in clayey and highly conductive soils^[13^]. An ERT/I system measures the distribution of electrical resistivity in the subsurface by sending a direct current signal into the ground through one set of electrodes and measuring the induced voltage through another set of electrodes. While the ERT/I technique is a powerful technique for investigating different depths of a transect, it tends to be time consuming and most suited to static deployment. An EMI system transmits a low frequency signal (kHz) into the subsurface without the need to establish contact with the ground. The alternating current produces a magnetic field in the subsurface, which in turn induces secondary current loops related to the subsurface electrical conductivity. These create secondary magnetic field loops and the instrument measures the superposition of the combined primary and secondary fields^[14,15^]. This noninvasive technique is appropriate for field-scale measurement due to its rapid response, ease of integration into mobile vehicular measuring platforms and nondestructive/non-contact requirements.

Apparent electrical conductivity (EC_a) is a proxy for subsurface physical properties and provides a measure of charge mobility due to an application of an electric field, and is defined as the ratio between current density (J, A·m⁻²) and electrical field (E, V·m⁻¹) according to Ohm’s law. Several physical and chemical soil properties influence field-scale EC_a measurements. Friedman^[16^] conveniently partitions the major factors into three categories: bulk soil, solid particle and soil solution. The bulk soil category comprises the factors that are defined by the organization of a three phase soil system such as porosity (n) and water content (q). Factors in the solid particle category include particle shape and orientation, particle-size distribution, cation exchange capacity (CEC), and wettability. Ionic strength (s_w), cation composition and temperature are factors that are classified under the soil solution category.

EMI-based EC_a measurements have been used by researchers attempting to infer different properties and characterize a wide range of processes at the field scale for a host of different applications^[17^]. Some of the applications include: estimating claypan depth^[18^]; predicting herbicide application leaching potentials in specific areas^[19^]; petrocalcic horizon depth^[20^], inferring topsoil depth in claypan soils^[21^]; and identifying the locations and depths of septic-system failure^[22^]. Corwin and Lesch^[23,24^] outline standard operating procedure for EC_a surveys applied to precision agriculture; specifically surveys which calibrate EC_a to the electrical conductivity of saturation extract samples (EC_e) for use in salinity studies, and discuss several different applications of field-scale EC_a maps.

There are some studies that have used EMI mapping techniques coupled with soil sampling to delineate subsurface properties. Greve and Greve^[25^] have applied EMI mapping to better define soil map unit delineation widths. In classical soil mapping, map unit transition zones are represented by lines since it is time consuming and labor intensive to exactly quantify the width of the transition zone^[25^]. Their study used auger sampling and applied a spatial rate of change calculation on a kriged EMI map to better define the transition zones of the map units. EMI mapping was used to infer the subsurface morphology of an agricultural field to identify areas with off-site agrochemical migration^[26^]. The authors mapped the field for two consecutive days to see the effect of moisture change on EC_a. They were able to infer that the EC_a pattern similarity observed at two different water contents near field capacity was due to soil morphology. Areas with a faster change in conductivity as the field dried were associated with high unsaturated hydraulic conductivity and as prime candidates for chemical migration. Kitchen et al.^[27^] investigated the effectiveness of using EMI mapping for delineating productivity zones for agricultural management in claypan soils of Missouri. The study showed that the productivity zones delineated using EC_a and elevation data agreed up to 70% with those delineated from 10 years of combine monitored yield maps. Farahani and Buchleiter^[28^] conducted multi-year EC_a surveys to classify sandy and non-saline fields into low, medium and high EC_a zones. They measured the temporal variability of one mapping event from another by investigating how the measurements deviated from the 1:1 line.

Different field experiments have shown that the major properties that contribute to the apparent electrical conductivity of a soil include the electrical conductivity of the soil solution (EC_e), water content and texture. Our experiment was designed to test EMI performance at the low end of EC_a measurements. Differences in EC_a between sand and 2:1 clay mineral soils often can range as much as 60 to 100 mS·m⁻¹. However, in many agricultural soils, textural differences may be more subtle, but still of hydrological or economic importance. We chose a 50 × 50 m field site with a relatively uniform silt-loam soil, but on an alluvial fan with relict, subsurface gravel channels that cannot be readily observed from the surface; our objective was to determine if the EMI was sensitive enough to identify the location of these channels. Changes in soil water content could potentially mask differences in soil texture and complicate the interpretation of and electrical mapping of the subsurface in terms of defining boundaries. Our strategy was to map the field using EMI at a range of soil water contents to try and identify locations with consistently higher or lower EC_a than the global mean. In addition, analysis was conducted to determine the optimal water content for identifying textural differences. Therefore, the objective of this study was to use repeated EMI mapping and temporal stability analysis at different soil water contents to delineate soil textural patterns in a low EC_a agricultural field.

2 Theoretical background

2.1 Temporal stability analysis

Vachaud et al.^[29^] characterize the time-invariant association between spatial location and classical statistical parametric values as the concept of temporal or rank stability. The method depends on a spatial location keeping its rank in the cumulative probability function for different sampling times^[29^]. In the case of soil water content this has been shown to work relatively well for level ground, but less so for sloping ground^[30,31^]. Since our study area was nearly level and the main assumption was that EMI mapping can capture a time-invariant subsurface physical property through repeated mapping, a temporal stability analysis technique was a good way of quantifying and analyzing the data. We have modified Vachaud et al.’s^[29^] equations as used for moisture storage to apparent electrical conductivity (EC_a).

For each support block (i) in the field on a mapping event (j), the apparent electrical conductivity was defined as ${{\displaystyle \mathrm{E}\mathrm{C}\mathrm{a}}}_{ij}$. The difference ${{\displaystyle \mathrm{\Delta }}}_{ij}$, (Eq. (1)), was evaluated by subtracting the average EC_a for all n locations at mapping event, j, $\overline{{{\displaystyle \mathrm{E}\mathrm{C}\mathrm{a}}}_j}$, (Eq. (2)), from ${{\displaystyle \mathrm{E}\mathrm{C}\mathrm{a}}}_{ij}$.

(1)${{\displaystyle \mathrm{\Delta }}}_{ij}={{\displaystyle \mathrm{E}\mathrm{C}\mathrm{a}}}_{ij}-\overline{{{\displaystyle \mathrm{E}\mathrm{C}\mathrm{a}}}_j}$

where

(2)$\overline{{{\displaystyle \mathrm{E}\mathrm{C}\mathrm{a}}}_j}={\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \sum _{i=\mathrm{1}}^n{{\displaystyle \mathrm{E}\mathrm{C}\mathrm{a}}}_{ij}}$

The array ${{\displaystyle \mathrm{\Delta }}}_{ij}$ was normalized by dividing it by $\overline{{{\displaystyle \mathrm{E}\mathrm{C}\mathrm{a}}}_j}$ to produce a new variable, the relative difference (${{\displaystyle \delta }}_{ij}$).

(3)${{\displaystyle \delta }}_{ij}={\displaystyle \frac{{{\displaystyle \Delta }}_{ij}}{\overline{{{\displaystyle ECa}}_j}}}$

For each mapping event, Eq. (3) was then used to evaluate a column of relative differences for all locations.

Once the relative differences were determined, we then calculated the average relative difference ($\overline{{{\displaystyle \delta }}_i}$) for each location across m mapping events, given by:

(4)$\overline{{{\displaystyle \delta }}_i}={\displaystyle \frac{\mathrm{1}}{m}}{\displaystyle \sum _{j=\mathrm{1}}^m{{\displaystyle \delta }}_{ij}}$

The standard deviation of the relative differences of a location can also be evaluated similarly. A quantitative measure for testing the time stability between two mapping days is the Spearman’s rank correlation coefficient (r_s)^[32^], given by:

(5)${{\displaystyle r}}_s(j,j^{\prime })=\mathrm{1}-{\displaystyle \frac{\mathrm{6}{\displaystyle {\sum }_{i=\mathrm{1}}^n{({{\displaystyle R}}_{ij}-{{\displaystyle R}}_{i{j}^{\prime }})}^{\mathrm{2}}}}{n(n^{\mathrm{2}}-\mathrm{1})}}$

The test, specified in Eq. (5), is evaluated by comparing the rank of a support block, i, on a specific mapping event j (R_ij) to its rank on another mapping event j’ (R_ij’). A rank correlation matrix can be constructed for all the mapping days by evaluating Eq. (5) for j, j’ = 1 to m, where j ≠ j’. The matrix, r_s(j,j’) is an upper triangular matrix with values of one in the diagonal. The closer the values of r_s(j,j’) are to 1, the stronger the temporal stability of all the locations in the field.

2.2 Data transformation using the normal score procedure

The underlying assumption of kriging is that the data are normally distributed^[33^]. The normal score transform is useful in normalizing many environmental variables that have large outlying values (positively skewed) to provide a normal distribution. The normal score transform function is derived by matching the original skewed cumulative distribution function (cdf) to a standard normal cdf (e.g., Fig. 1).

Functions F(z) and G(y) are cdfs of the original random function (RF) Z(x) and the standard normal RF Y(x) respectively^[34^]:

(6)$F(z)=\mathrm{Pr}\mathrm{o}\mathrm{b}\{Z(x)\le z\},G(y)=\mathrm{Pr}\mathrm{o}\mathrm{b}\{Y(x)\le y\}$

The transform, j(.), that takes any cdf F(z) to a standard Gaussian cdf G(y) is given as:

(7)$Y(x)=\phi (Z(x))=G^{-\mathrm{1}}[F(Z(x))]$

where G⁻¹(.) is the inverse Gaussian cdf or quantile function of the RF Y(x).

Fig.1 Normal score transformation (j) of the skewed cumulative distribution function (CDF), F(z), from ECa mapping event of July 9, 2006 (a) into a normal CDF G(y) (b).

Full size|PPT slide

In algorithm form, the normal score transform procedure can be reduced to the following three steps^[35^]:

(1) The n original data values, z(x), are first ranked in ascending order and tied z values are separated according to the local averages of the data surrounding each tied value.

(2) The sample cumulative frequency of the datum z(x) with rank k is then computed as:

(8)${{\displaystyle p}}_k^{*}=k/n-\mathrm{0.5}/n$

where all data receive the same weight 1/n

(3) The normal score transform of the z datum with rank k is matched to the p_k*-quantile of the standard normal cdf:

(9)$y(x)=G^{-\mathrm{1}}[F^{*}(Z(x))]=G^{-\mathrm{1}}({{\displaystyle p}}_k^{*})$

2.3 Spatial prediction

Kriging relies on an underlying spatial structure of a measured variable in order to predict its value at unsampled locations^[33,36^]. Most prediction methods, including kriging, average the weighted values of the adjacent sampled values (z(x_i)) in order to predict the variable at the unsampled point z*(x₀).

(10)$z^{*}({{\displaystyle x}}_{\mathrm{0}})={\displaystyle {\sum }_{i=\mathrm{1}}^n{{\displaystyle \lambda }}_i·z}({{\displaystyle x}}_i)$

The kriging problem simplifies to solving for a vector of weights, l_i, that will minimize a generalized least squares (GLS) equation. The spatial dependence of the process, represented in the residuals of the GLS regression equation, is solved when:

(11)$\mathit{\lambda }={\displaystyle \frac{c(\mathit{x})}{C}}$

where C is the matrix of covariances, C(x_i,x_j), between all possible pairs of the n sample sites and c(x_i) is a column vector of covariances between the prediction point and each of the n sample sites.

To solve for l, we needed to evaluate the matrix of covariances C, which was done using a semivariogram function, written:

(12)$\gamma (h)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}N(h)}}{\displaystyle {\sum }_{i=\mathrm{1}}^{N(h)}{[z({{\displaystyle x}}_i)-z({{\displaystyle x}}_i+h)]}^{\mathrm{2}}}$

where the function computes the average squared differences of the values of the random variable z(.) at a vector of data pairs x and x+h, where N(h) is a number of data pairs within a given class of distance. A parametric model was used to describe the experimental semivariogram to provide a continuous, positive and smooth description of the covariance matrix, C.

Block kriging extends the above method from a point estimation of a spatially continuous variable to the average value over a small area or block^[34^]. This is useful when the support block of a physical measurement is beyond a point as in EMI measurements.

3 Materials and methods

3.1 Study site

The field study was conducted in the summer and autumn of 2006 and 2007 at Greenville Farm, Utah Agricultural Experiment Station, North Logan, Utah (41°46′1′′ N, 111°48′40′′ W). Mean annual precipitation is 422 mm and the mean annual temperature is 8.6°C, providing the site with a xeric soil moisture regime and a mesic soil temperature regime. The area of study was a 50 m × 50 m square area with smooth, nearly level topography located on the eastern end of a larger field. The field had been fallow for 2 years at the time of study.

The soil at the Greenville Farm is of the Millville Series. The soil parent material is medium-textured alluvium (silts and very fine sand) deposited as distal fan and overbank flood deposits on the Green Canyon alluvial fan (Fig. 2a). The soil pH was 8.2 due to the high concentration of primary CaCO₃. The soil is classified in the family of coarse-silty, carbonatic, mesic Typic Haploxerolls. The typical pedon that represents the soil in most of the study area contains<1% rock fragments and the texture (silt loam) is uniform with depth. However, gravel lenses exist in these fields and become most apparent during soil ripping (J. Slade, personal communication, 2006). Discussion with the farm manager and our visual observations were used to identify the general location of the gravel lens (Fig. 2b) in what is otherwise a relatively homogeneous soil in the study field. Representative pedons of the gravel-free soil and the soil with the gravel lens were described in backhoe-excavated pits using standard methods.

Fig.2 (a) Hillshade image showing the location of the study area (white box) on the distal portion of the alluvial fan (outlined with dotted line) deposited by streams draining Green Canyon (white line); (b) the 50 m × 50 m field showing the EMI survey route (dotted line), the general location of the gravel lens (crosshatch), and locations of the two soils (pits 1 and 2) described (white boxes).

Full size|PPT slide

3.2 Instruments

A set of 11 0.15-m TDR probes with thermocouples were set up in a plot near the mapping location in soil representing the gravel-free Millville series in order to measure water content and temperature. The sensors were buried at depths of 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.7, 1.0, 1.5 and 2.0 m below the surface. Volumetric water content, bulk electrical conductivity, and temperature data were collected every 30 min. The volumetric water content (q_v) at 0.5 m was taken to be the water content of the control plot, corresponding to the depth of maximum weighting for the EMI measurement. Hourly rainfall and evapotranspiration data were also collected from a weather station located on an adjacent grass field.

Georeferenced EC_a measurements were taken noninvasively using the DUALEM-1S ground conductivity instrument (Dualem, Milton, ON, Canada) along with a Trimble ProXT GPS unit (Sunnyvale, CA, USA). The DUALEM-1S (DUAL-geometry electromagnetic) is a geo-conductivity sensor with a transmitter operating at the frequency of 9 kHz and two receivers with different orientations. We used the horizontal coplanar geometry or V-V_DLM mode where both the transmitter and the receiver, with a 1-m separation, use vertical dipoles. The depth of exploration for the V-V_DLM setup is about 1.5 m^[15^], with the depth of maximum weighting at about 0.5 m. Data from the EMI instrument was transmitted serially through a 9-socket DB-9 connector and was acquired simultaneously with the GPS data using a hand-held geographic information system (HGIS, StarPal Inc., Fort Collins, CO, USA) program inside an Allegro CX hand-held field computer (Juniper Systems, Logan, UT, USA). In a previous study we tried using a Geonics EM-38, but the required calibration procedure to null out the magnetic susceptibility circuit resulted in inconsistent readings at low EC_a values (<20 mS·m⁻¹). The DUALEM-1S only measures EC_a and comes with factory set calibration, resulting in stable, consistent readings in low EC_a soils.

3.3 Mapping

The EMI instrument was turned on for 30 min before mapping, which was usually carried out early in the morning or in the evening, to avoid temperature drift effects on the instrument^[15^]. The EMI instrument was held at about 10 cm above ground while traversing the field by walking in rows spaced 3 m apart (Fig. 2b). The EMI mapping process required about 45 min with EC_a data being collected every second. The field was mapped a total of 30 times in summer/autumn 2006 and 2007.

Maps were selected from six different days to cover a range of field moisture regimes as recorded in the control plot. The volumetric water content (q_v) ranged from a low of 0.11 (September21, 2007) to a high of 0.23 (October 7, 2006). Intermediate water-content mapping days had the following q_v: 0.13 (October 15, 2007); 0.16 (July 9, 2006); 0.18 (September 28, 2006); and 0.20 (September 26, 2006).

To reduce the effect of diurnal temperature fluctuations, the EC_a data for each day were corrected to a standard temperature of 25°C. The temperature corrected EC_a data were then checked for continuity and anomalous values using a time-series view of the data. Anomalous values can be caused by buried metal fragments, wires, pipes and similar objects. These were identified and removed from the data set as a quality control measure.

3.4 Data analysis

We carried out exploratory data analysis to produce basic statistics and histograms for each mapping event. The quality controlled EC_a data were normal score transformed using S-GeMS^[37^] to prepare the data for a kriging process. The field geometry was then subdivided into 2500 blocks of 1m² area representing the support area of the EMI instrument. The data was kriged using VESPER^[38^] and then returned to S-GeMS to be back transformed. Temporal stability analysis was then applied to the six EC_a maps. Since our support area is a 1m² block, each map was divided into 2500 zones and the relative difference of each zone was calculated for each mapping event. The relative differences for each zone were then averaged across the six mapping events (mean relative difference) and ranked in ascending order.

Fig.3 Block kriged ECa maps of the field at varying water content with each shade gradation representing 12.5 percentile of the overall data, with the lightest color representing low ECa and the darkest color representing high ECa values. The general location of the gravel area is inside the dashed line enclosure.

Full size|PPT slide

4 Results

The six selected mapping days representing different volumetric water contents were used to produce field-scale EC_a maps presented in Fig. 3. The maps show the EC_a of the field at varying water content with each color gradation on the map representing a 12.5 percentile of the data range, with the lightest color representing low EC_a and the darkest color representing high EC_a values. The average EC_a of the field for each volumetric water content is shown in Fig. 3 and ranged from 7.60 mS·m⁻¹ at a q_v of 0.11 m³·m⁻³ to 14.7 mS·m⁻¹ for the highest water content of 0.23 m³·m⁻³.

Visual inspection of the figures shows the low EC_a area corresponding with the general location of the gravel lens and a higher EC_a region to the north. It is also noticeable that the pattern around the low EC_a area changes with water content. The map in Fig. 3e shows the least consistent pattern. The EC_a was measured at the end of a rainfall event and we interpret the pattern to indicate the redistribution of water changing the soil EC_a as the field was mapped. This indicates that mapping should be avoided during or immediately following rainfall events when rapid changes in EC_a may occur with wetting.

Further analysis of the data was conducted using the coefficient of variation (the ratio of the standard deviation to the average EC_a of the field). This tended to be lowest at the low and high water contents (Fig. 3a, Fig. 3b, Fig. 3f) and highest at the somewhat intermediate water contents (Fig. 3c, Fig.3d, Fig. 3e). This perhaps indicates that the greatest contrast for a single mapping occurs at an intermediate water content, discussed later on.

Box plots are presented in Fig. 4 and show the EC_a distribution for each mapping event as the volumetric water content increased. A clear increase in the field average EC_a is observed as a function of water content. The lowest interquartile range (IQR, Q3-Q1) of 1.23 mS·m⁻¹ was observed at the lowest q_v of 0.11 and the highest IQR of 2.75 mS·m⁻¹ was observed at the highest q_v of 0.23. The standard deviation increased at a rate lower than the rate of increase of the mean because the CV did not increase with increasing average EC_a, likely an underlying reflection of the quadratic relationship between EC_a and soil water content.

Fig.4 Box plots showing the ECa distribution of each mapping event as a function of volumetric water content. The box plots give the 25 (Q1), 50 (Q2, median) and 75 (Q3) percentile of the data in the box as well as the 5 and 95 percentile of the data at the whiskers. Minimum and maximum values are shown by dots.

Full size|PPT slide

Spearman’s rank correlation coefficient (r_s) was used to get a quantitative measure of the time stability of spatial locations between different mapping days (Table 1). The three wetter days of mapping (q_v = 0.18, 0.20 and 0.23) had relatively high Spearman’s rank correlation coefficients. The highest occurs between mapping events with q_v of 0.18 and 0.20 with r_s = 0.92. The two driest days of mapping (water contents of 0.11 and 0.13) have the second highest r_s in Table 1 at 0.89. The wet and dry mapping days that have the highest correlation coefficients are those mapped a few days apart. The mapping events at water contents just below field capacity (e.g. q_v = 0.16 and q_v = 0.20) are the least correlated with the other mapping events, indicating the strongest electrical contrast.

Tab.1 Spearman’s Rank Correlation coefficient between the six mapping events at varying volumetric water content, qv

Volumetric water content	qv = 0.11	qv = 0.13	qv = 0.16	qv = 0.18	qv = 0.20	qv = 0.23
qv = 0.11	1
qv = 0.13	0.89	1
qv = 0.16	0.79	0.57	1
qv = 0.18	0.69	0.79	0.54	1
qv = 0.20	0.54	0.68	0.43	0.92	1
qv = 0.23	0.78	0.72	0.84	0.86	0.77	1

Fig.5 Temporal stability analysis. (a) The mean relative difference (%) of the spatial blocks ranked in ascending order; (b) a map of the mean relative difference delineated into three regions according to the inflection points of the mean relative difference graph.

Full size|PPT slide

The relative differences for each zone were averaged across the six mapping events (mean relative difference) and ranked in ascending order to produce Fig. 5a. The lowest ranked spatial location (zone) is 27% lower than the averaged mean EC_a of the six mapping events, while the highest ranked zone is 18% higher. We interpreted the inflection points as indicating transition zones between soil units, representing boundary delineations to classify the field into three regions (Fig. 5b). The regions were classified using the following ranges of mean relative difference percentages: Region 1 (lowest conductivity area), -27% to -11%; Region 2, -11% to 5%; and Region 3 (highest conductivity area), 5% to 18%.

Fig.6 Selected morphological properties of soils described at the Millville control soil representative of Region 3 (a) and the soil with gravelly horizons representative of Region 1 (b), pits 1 and 2 respectively on Fig. 2. Colors are for moist soil; gravel volume was estimated visually in the field; and clay concentration was estimated by feel.

Full size|PPT slide

Ground-truthing was conducted to determine observable differences in these regions. The two soils described initially are consistent with Regions 1 and 3. Selected morphological properties of the soils are illustrated in Fig. 6. The Millville pedon (coarse-silty, carbonatic, mesic Typic Haploxerolls) of the control plot (Fig. 6a) represents a typical pedon of the gravel-free, silt-loam dominated Millville series found in Region 3. The pedon with the gravel lens (coarse-loamy, carbonatic, mesic Entic Haploxerolls) typical of Region 1 (Fig. 6b) has a substantial amount of gravel throughout the upper 65 cm, with a gravel lens (80% gravel by volume) from 45 to 65 cm.

Core samples representative of the three zones were taken in 30 cm increments. In Region 1, where we expected the gravelly soil, the coring device was only able to penetrate to a depth of about 40 cm until impeded by the gravel, so only one core was sampled. In Region 2, the coring device was able to penetrate deeper, going to depths of between 70 and 85 cm, with the depth increasing with distance from Region 1. Three cores were sampled in Region 2. In Region 3, five cores were taken to the extent of the device (90 cm), indicating essentially no gravel in this zone.

Each 30-cm core sample was homogenized and analyzed for particle-size distribution using the hydrometer method^[39^]. In addition to the clay percentage, the percentage of fine particles (silt and clay) was determined due to its importance in the water holding capacity of the soil. Region 2 had an average clay percentage of 11.9%±1.09%, whereas Region 3 had an average clay percentage of 12.9%±1.82%. We did not find a statistically significant difference between the average clay percentages of Regions 2 and 3. However, the average fine particle percentages in Regions 2 and 3 were 50.5%±4.39% and 54.6%±4.09%, respectively. We were able to reject the null hypothesis of equal means of fine particle percentage for both regions at the 5% significance level (P = 0.023) using a two sample t-test. The test was repeated with the nonparametric Wilcoxon rank sum test, and we were able to ascertain at the 5% significance level (P = 0.036) that the medians of fine particle percentage were different between Regions 2 and 3.

We interpret Region 1 as the gravel lens, which tapers off in Region 2, giving way to the homogeneous Millville silt-loam in Region 3. Region 1 is consistent with a relict high energy channel of the alluvial fan. As the energy of the channel subsided, finer textured material was deposited forming the present Millville Series soil. We can also see from our textural delineation map (Fig. 5b) that the coarser texture material also spread at the edge of the high energy channel depositing a limited amount of gravel around the surrounding area until tapering off. The spread of the gravel was predisposed toward the south to Region 2, where our depth of coring was limited and the fine particle percentage is lower compared to Region 3. We conclude that EMI imaging using repeated surveys at different water contents can be a useful tool in identifying subsurface morphological features such as gravel lenses.

Fig.7 Histograms of ECa at qv = 0.11 (solid line); qv = 0.23 (dashed line); and, qv = 0.16 (dotted line) field water content conditions

Full size|PPT slide

5 Optimal water content for mapping heterogeneity

Clearly, a multi-mapping strategy provides important information for delineating soil textural boundaries without the extra work required for physical calibration. However, repeated mapping cannot always be achieved and so a pertinent question to ask is, for soils with low electrical contrasts, like non-saline soils, is there an optimum water content for determining underlying spatial patterns of soil texture. To attempt to answer this question, we examined the histograms of the kriged maps (Fig. 7). We focused on the driest mapping event (q_v = 0.11), wettest mapping event (q_v = 0.23), and a mapping event in between with a water content of 0.16 (medium). The histogram for q_v = 0.16, shows three strong peaks, which appear to collapse and merge at the high and low water contents. This may indicate that the strongest contrast in EC_a occurs at water contents just below field capacity (q_v ~ 0.21).

If physically or economically limited to one mapping event, the strongest texturally induced EC_a differences in the Millville soil were observed below field capacity. This point (q_v = 0.16) has the highest coefficient of variation, 16% (Fig. 3), it is the least correlated with the other maps (Table 1), and it exhibits a multimodal histogram with distinct peaks (Fig. 7). Why does this occur? Theory and research supports the idea that field variability in bulk EC likely increases with the variability of water content^[40^]. The greatest electrical contrast will occur when one textural component is wetter and one drier. So when sand and clay areas are subjected to gravitational drainage we expect the sand to drain while the clay will remain wetter. This increases the electrical contrast between textural areas, and appears to occur close to field capacity, which is what we are exploiting to map the texture. However, there is a further issue that arises, which is how changes in the soil moisture due to drainage alter the vertical electrical conductivity. When the nonlinear nature of the spatial weighting function of the EM instrument comes into play^[41^], and vertical heterogeneity in soil texture and water content is added, it becomes harder to predict how the spatial variability of the bulk EC will respond to spatial and vertical variability in soil water content. This raises the question as to what extent the high variability in EC_a at q_v = = 0.16 is a result of the variability in the textural component of the EC_a or a result of the interaction in horizontal and vertical variations in water content and the spatial weighting function of the EM instrument. Given that higher EC_a tends to compress the EM depth of exploration^[41^], this perhaps adds to the contrast. This paper does not address this intriguing question which will set the limits on how well textural contrasts can be imaged using the EMI approach, especially using inversion. Combining techniques such as EMI and gamma ray may well prove to be the best approach^[42^].

The repeated EMI mapping of low EC_e soils at varying water content reveals the textural patterns of the subsurface as demonstrated by this study. The fact that the range of EC_a in this research was only about 12 mS·m⁻¹ makes the adaptation of this methodology into areas with a larger EC_a range more informative. Thus repeated EMI mapping can be useful in soil surveys intending to delineate areas with heterogeneous soils as well as in better defining transition zones between soil units. The methodology should be considered a tremendous benefit in the arsenal of tools used by the soil surveyor, especially for site-specific soil maps. All geophysical techniques exploit contrasts in target properties and as such require careful interpretation. The methodology could be useful in precision agriculture in demarcating productivity and management zones for improved utilization of resources and better yield without the need for extensive calibration. The methodology could also be used to improve sampling schemes, especially in pristine environments, by providing an extra layer of information on soil variability and determining locations where maximum or minimum change occurs. This information could be very helpful for the potential placement of monitoring equipment, sensors and observation nodes for monitoring soil hydrological processes in situ.

6 Conclusions

Spatial variability of soil properties presents difficulty for capturing the heterogeneity without sampling exhaustively and becoming overwhelmed with data analysis, and losing site of the dominant processes. The development and testing of a new procedure for mapping non-saline soils to differentiate static soil characteristics from dynamic ones was presented. We considered soil texture as a dominant static property, whereas soil water content is a dynamic property that changes rapidly over time. To separate these properties, we developed a multi-mapping methodology and analysis procedure. This allowed us to identify locations exhibiting consistent behavior over a range of soil water content. The procedure allowed us to noninvasively map field-scale soil textural patterns by separating the EMI response due to water content variation from static textural properties using temporal stability analysis. We collected and presented six georeferenced EC_a surveys at volumetric water contents ranging from 0.11 to 0.23 in a 50 × 50 m field containing silt loam soil. Block kriging facilitated prediction of unsampled areas to produce EC_a maps at meter resolution. Temporal stability analysis was then applied on the six EC_a maps and the field was divided into three regions. The lowest conductivity region was found to be associated with a relict high energy channel that deposited coarser materials (gravel) as the soil parent material. This noninvasive mapping approach has the potential to reveal the spatial distribution of time-invariant subsurface properties using repeated EMI surveys, especially when taken over a range of field soil moisture levels. Given the composite maps produced, decisions can be made as to whether calibration is required, or if the obtained information is sufficient for the purposes of identifying management zones. These maps provide valuable insights into soil textural properties to support modeling and experimental design for a broad array of disciplines.

References

[1] Corwin D L, Plant R E. Applications of apparent soil electrical conductivity in precision agriculture. Computers and Electronics in Agriculture, 2005, 46(1–3): 1–10

[2] Robinson D A, Binley A, Crook N, Day-Lewis F D, Ferre T P A, Grauch V J S, Knight R, Knoll M, Lakshmi V, Miller R, Nyquist J, Pellerin L, Singha K, Slater L. Advancing process-based watershed hydrological research using near-surface geophysics: a vision for, and review of, electrical and magnetic geophysical methods. Hydrological Processes, 2008, 22(18): 3604–3635

[3] Binley A, Hubbard S S, Huisman J A, Revil A, Robinson D A, Singha K, Slater L D. The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales. Water Resources Research, 2015, 51(6): 3837–3866

[4] Huang J, Barrett-Lennard E G, Kilminster T, Sinnott A, Triantafilis J. An error budget for mapping field-scale soil salinity at various depths using different sources of ancillary data. Soil Science Society of America Journal, 2015, 79(6): 1717–1728

[5] Hillel D. Environmental soil physics. San Diego: Academic Press, 1998

[6] Davie T. Fundamentals of hydrology. London: Routledge, 2003

[7] Huang J, McBratney A B, Minasny B, Triantafilis J. Monitoring and modelling soil water dynamics using electromagnetic conductivity imaging and the ensemble Kalman filter. Geoderma, 2017, 285: 76–93

[8] Huang J, Monteiro Santos F A, Triantafilis J. Mapping soil water dynamics and a moving wetting front by spatiotemporal inversion of electromagnetic induction data. Water Resources Research, 2016, 52(11): 9131–9145

[9] Marquez O, Julia GdB. Soil taxonomy: a basic system of soil classification for making and interpreting soil surveys. Soil Science, 1977, 123(4): 270

[10] Reynolds J M. An introduction to applied and environmental geophysics. New York: John Wiley & Sons, 1997

[11] Telford W M, Geldart L P. R.E.Sheriff R E. Applied geophysics. 2nd ed. New York: Cambridge University Press, 1990

[12] Rubin Y, Hubbard S S. Hydrogeophysics, New York: Springer Verlag, 2006

[13] Weihermuller L, Huisman J A, Lambot S, Herbst M, Vereecken H. Mapping the spatial variation of soil water content at the field scale with different ground penetrating radar techniques. Journal of Hydrology, 2007, 340(3–4): 202–216

[14] McNeill J D. Electromagnetic terrain conductivity measurement at low induction numbers.  Geonics Limited, 1980

[15] Abdu H, Robinson D A, Jones S B. Comparing bulk soil electrical conductivity determination using the DUALEM-1S and EM38-DD electromagnetic induction instruments. Soil Science Society of America Journal, 2007, 71(1): 189–196

[16] Friedman S P. Soil properties influencing apparent electrical conductivity: a review. Computers and Electronics in Agriculture, 2005, 46(1–3): 45–70

[17] Hendrickx J M H, Kachanoski R G. Nonintrusive electromagnetic induction, In: J.H. Dane, and G.C. Topp, ed. Methods of Soil Analysis: Part 4 Physical Methods. Madison: Soil Science Society of America, 2002, 1297–1306

[18] Doolittle J A, Sudduth K A, Kitchen N R, Indorante S J. Estimating depths to claypans using electromagnetic induction methods. Journal of Soil and Water Conservation, 1994, 49(6): 572–575

[19] Jaynes D B, Novak J M, Moorman T B, Cambardella C A. Estimating herbicide partition coefficients from electromagnetic induction measurements. Journal of Environmental Quality, 1995, 24(1): 36–41

[20] Boettinger J L, Doolittle J A, West N E, Bork E W, Schupp E W. Non-destructive assessment of rangeland soil depth to petrocalcic horizon using electromagnetic induction. Arid Soil Research and Rehabilitation, 1997, 11(4): 375–390

[21] Sudduth K A, Drummond S T, Kitchen N R. Accuracy issues in electromagnetic induction sensing of soil electrical conductivity for precision agriculture. Computers and Electronics in Agriculture, 2001, 31(3): 239–264

[22] Taylor R, Lee B, Doolittle J. Ground conductivity survey of a septic system during and after failure. Leading Edge, 2003, 22(6): 555–566

[23] Corwin D L, Lesch S M. Application of soil electrical conductivity to precision agriculture: theory, principle and guidelines. Agronomy Journal, 2003, 95(3): 455–471

[24] Corwin D L, Lesch S M. Apparent soil electrical conductivity measurements in agriculture. Computers and Electronics in Agriculture, 2005, 46(1–3): 11–43

[25] Greve H G, Greve M B. Determining and representing width of soil boundaries using electrical conductivity and MultiGrid. Computers & Geosciences, 2004, 30(6): 569–578

[26] Wilson R C, Freeland R S, Wilkerson J B, Yoder R E. Inferring subsurface morphology from transient soil moisture patterns using electrical conductivity. Transactions of the American Society of Agricultural Engineers, 2004, 6(5): 1435–1441

[27] Kitchen N R, Sudduth K A, Myers B, Drummond S T, Hong S Y. Site-specific productivity zones delineated using bulk soil electrical conductivity. In:ASAW Conference Paper: 2003. Las Vegas: American Society of Agricultural Engineers, 2003

[28] Farahani H J, Buchleiter G W. Temporal stability of soil electrical conductivity in irrigated sandy fields in Colorado. Transactions of the American Society of Agricultural Engineers, 2004, 47(1): 79–90

[29] Vachaud G, Passerat Silans A, Balabanis P, Vauclin M. Temporal stability of spatially measured soil water probability density function. Soil Science Society of America Journal, 1985, 49(4): 822–828

[30] Kachanoski R G, de Jong E. Scale dependence and the temporal persistence of spatial patterns of soil water storage. Water Resources Research, 1988, 24(1): 85–91

[31] Grayson R B, Western A W. Towards areal estimation of soil water content from point measurements: time and space stability of mean response. Journal of Hydrology, 1998, 207(1–2): 68–82

[32] Kottegoda N T, Rosso R. Probability, statistics, and reliability for civil and environmental engineers. Boston: McGraw-Hill, 1997

[33] Webster R, Oliver M. Geostatistics for environmental scientists. New York: John Wiley and Sons, 2000

[34] Deutsch C V, Journel A G. GSLIB: Geostatistical Software Library and User’s Guide: 2nd ed.  New York: Oxford University Press, 1998

[35] Goovaerts P. Geostatistics for natural resources evaluation. New York: Oxford University Press, 1997

[36] Goovaerts P. Geostatistics in soil sciences: State-of-the-art and perspectives. Geoderma, 1999, 89(1–2): 1–45

[37] Remy N. S-GeMS: The Stanford geostatistical modeling software: a tool for new algorithms development, In: Leuangthong O and Deutsch C V. Eds. 7th Int. Geostatistics Conference, Quantitative Geology and geostatistics, New York: Springer, 2005

[38] Walter C, McBratney A B, Douaoui A, Minasny B. Spatial prediction of topsoil salinity in the Chelif valley, Algeria, using local kriging with local variograms versus local kriging with whole-area variogram. Australian Journal of Soil Research, 2001, 39(2): 259–272

[39] Gee G W, Bauder J W. Particle size analysis. In: Klute A, ed. Methods of Soil Analysis: Part 1. Physical and Mineralogical Methods: Second edition. Madison: Agronomy Society of America, 1986, 383–411

[40] Robinson D A, Abdu H, Lebron I, Jones S B. Imaging of hill-slope soil moisture wetting patterns in a semi-arid oak savanna catchment using time-lapse electromagnetic induction. Journal of Hydrology, 2012, 416–417: 39–49

[41] Callegary J B, Ferré T P A, Groom R W. Vertical spatial sensitivity and exploration depth of low-induction-number electromagnetic-induction instruments. Vadose Zone Journal, 2007, 6(1): 158–167

[42] Huang J, Lark R M, Robinson D A, Lebron I, Keith A M, Rawlins B, Tye A, Kuras O, Raines M, Triantafilis J. Scope to predict soil properties at within-field scale from small samples using proximally sensed g-ray spectrometer and EM induction data.Geoderma, 2014, 232–234: 69–80

Acknowledgements

This research was supported by an INRA SSGP Doctoral Fellowship from the Inland North-west Research Alliance (INRA) and by the Utah Agricultural Experiment Station, Utah State University, Logan, Utah, approved as journal paper 8952. The authors express gratitude to Ian Riley for technical and grammatical corrections and to two anonymous reviewers who provided questions and suggestions for improving the manuscript.

Compliance with ethics guidelines

Hiruy Abdu, David A. Robinson, Janis Boettinger, and Scott B. Jones declare that they have no conflict of interest or financial conflicts to disclose.

This article does not contain any studies with human or animal subjects performed by any of the authors.

RIGHTS & PERMISSIONS

The Author(s) 2017. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0