Mapping Within-Field Soil Health Variations Using Apparent Electrical Conductivity, Topography, and Machine Learning

: High-resolution maps of soil health measurements could help farmers ﬁnetune input resources and management practices for proﬁt maximization. Within-ﬁeld soil heath variations can be mapped using local topography and apparent electrical conductivity (ECa) as predictors. To address these issues, a study was conducted in Texas Blackland Prairie soils with the following objectives: (i) to assess and map within-ﬁeld soil health variations using machine learning; (ii) to evaluate the usefulness of topography and ECa as soil health predictors; and (iii) to quantify the relationship between ECa and soil health index and use ECa to estimate soil health spatial distribution. We collected 218 topsoil (0–15 cm) samples following a 35 m × 35 m grid design and analyzed for one-day CO 2 , organic C, organic N, and soil health index (SHI) based on the Haney Soil Health Tool. A random forest model was applied to predict and map those properties on a 5 m × 5 m grid where ECa, and terrain attributes were used as predictors. Furthermore, the empirical relationship between SHI and ECa was established and mapped across the ﬁeld. Results showed that the study area was variable in terms of one-day CO 2 , organic C, organic N, SHI, and ECa distribution. The ECa, wetness index, multiresolution valley bottom ﬂatness, and topographic position index were among the top predictors of soil health measurements. The model was sufﬁciently robust to predict one day CO 2 , organic C, organic N (R 2 between 0.24–0.90), and SHI (R 2 between 0.47–0.90). Overall, we observed a moderate to strong spatial dependency of soil health measurements which could impact within-ﬁeld yield variability. The study conﬁrmed the applicability of easy to obtain ECa as a good predictor of SHI, and the predicted maps at high resolution which could be useful in site-speciﬁc management decisions within these types of soils.


Introduction
Soil health is defined as, "the capacity of soil to function as a vital living system, within an ecosystem and land-use boundaries, to sustain plant and animal productivity, maintain or enhance water and air quality, and promote plant and animal health" [1]. Concepts of soil health are available from Doran and Safley [2], Kibblewhite, et al. [3], and Maharjan, et al. [4], among others. Indeed, healthy soils are these soils with sufficient mineral nutrients, moisture, and growth-promoting microorganisms to increase productivity and promote human health [5,6]. Soil health indicators provide a quantitative approach to assessing soil quality and include a range of key soil properties [7]. Soil properties such as texture, bulk density, pH, soil organic matter, water holding capacity, drainage, and soil organism's functional relationship and diversity are the basic indicators of soil health [3]. Soil understand soil health and crop yield variations. The main objectives of this study were to: (i) assess and map within-field variability of selected soil health indicators using machine learning technique; (ii) evaluate the usefulness of local topography and ECa as predictors of soil health; and (iii) map the correlation between ECa and soil health and estimate the later with ECa collected on-the-go.

Site Description
The study site is located in the Blackland Prairies soils in Texas, USA (31 • 03 29.81 N, 97 • 21 15.08 W) ( Figure 1A). The site is a rectangular-shaped agriculture field covering an area of about 27 ha. The field typically grows corn and cotton in rotation and has been under cultivation for >50 year. The field applies conservation practices, and the two most common around the area are terracing to reduce soil erosion and minimum tillage. Soils in the study area are >2 m deep and are developed on weathered chalk or clayey residuum from calcareous mudstone of the Cretaceous Age. These soils are classified as fine-silty, carbonatic, thermic Udorthentic Haplustolls (Austin silty clay series) or fine, smectitic, thermic, Udic Haplusterts (Houston Black clay series) according to Soil Taxonomy [37]. Surface elevation changes from 185 m in the west up to 200 m to the east with a moderate slope up to 3.5%. The climate is generally hot with an average annual high and low temperatures of 25 • C and 13 • C, respectively, with an average precipitation of about 937 mm annually.
Agronomy 2022, 12, x FOR PEER REVIEW 3 of 16 about 27 ha. The field typically grows corn and cotton in rotation and has been under cultivation for >50 year. The field applies conservation practices, and the two most common around the area are terracing to reduce soil erosion and minimum tillage. Soils in the study area are >2 m deep and are developed on weathered chalk or clayey residuum from calcareous mudstone of the Cretaceous Age. These soils are classified as fine-silty, carbonatic, thermic Udorthentic Haplustolls (Austin silty clay series) or fine, smectitic, thermic, Udic Haplusterts (Houston Black clay series) according to Soil Taxonomy [37]. Surface elevation changes from 185 m in the west up to 200 m to the east with a moderate slope up to 3.5%. The climate is generally hot with an average annual high and low temperatures of 25 °C and 13 °C, respectively, with an average precipitation of about 937 mm annually.

Data Collection
Three sets of input data, namely, soil health indicators (one-day CO2, organic C, organic N, and SHI), apparent electrical conductivity, and topographic data, were collected from the field, and are described below.

Soil Sampling and Measurement of Soil Health Indicators
The soil sampling scheme in the field was based on a 35 m × 35 m grid layout in which one undisturbed soil sample at each grid node (218 nodes in total) was collected with a hydraulic soil probe containing an insert. Geographical coordinates of each sample location were recorded using a hand-held global positioning system with an accuracy of 0.76 m. Soil samples from the top 15-cm soil depth were then brought to the USDA-ARS, Grassland, Soil, and Water Research Laboratory to determine soil health indicators, and SHI based on Haney Soil Health Tool [11]. Those indicators include one-day CO2 evolution, water-extractable organic C (WEOC), and water-extractable organic N (WEON). One-day CO2 evolution is the amount of CO2 (mg/kg) released in the 24 h time period by soil microbes after the soil has been dried and rewetted, and it is a measure of the soil microbial activity related to soil fertility. The WEOC is the amount of organic C that is extracted from the soil with water, and it reflects the quantity of carbon in the soil that is readily available to microbes.

Data Collection
Three sets of input data, namely, soil health indicators (one-day CO 2 , organic C, organic N, and SHI), apparent electrical conductivity, and topographic data, were collected from the field, and are described below.

Soil Sampling and Measurement of Soil Health Indicators
The soil sampling scheme in the field was based on a 35 m × 35 m grid layout in which one undisturbed soil sample at each grid node (218 nodes in total) was collected with a hydraulic soil probe containing an insert. Geographical coordinates of each sample location were recorded using a hand-held global positioning system with an accuracy of 0.76 m. Soil samples from the top 15-cm soil depth were then brought to the USDA-ARS, Grassland, Soil, and Water Research Laboratory to determine soil health indicators, and SHI based on Haney Soil Health Tool [11]. Those indicators include one-day CO 2 evolution, water-extractable organic C (WEOC), and water-extractable organic N (WEON). One-day CO 2 evolution is the amount of CO 2 (mg/kg) released in the 24 h time period by soil microbes after the soil has been dried and rewetted, and it is a measure of the soil microbial activity related to soil fertility. The WEOC is the amount of organic C that is extracted from the soil with water, and it reflects the quantity of carbon in the soil that is readily available to microbes. Similarly, WEON is the amount of the total water-extractable nitrogen (WEN) minus the inorganic N that includes both nitrate and ammoniacal-nitrogen (NO 3 -N and NH 4 -N), and it can be easily broken down and released by soil microbes to the soil in readily plant-available inorganic N.
In the laboratory, soil samples were first oven-dried at 50 • C and crushed before sieving through a 2 mm mesh. From the sieved soil, a 4 g subsample was poured into a 50-mL centrifuge tube and was extracted with 40 mL de-ionized water. The extract was then analyzed for NO 3 -N and NH 4 -N using a Seal Analytical Rapid Flow Analyzer. The WEOC and WEN were determined using a Teledyne-Tekmar Apollo 9000 C: N analyzer and the WEON was derived as a difference between the total water-extractable N and nitrate-and ammoniacal-nitrogen (i.e., WEON = WEN-NH 4 -N-NO 3 -N). To determine one day CO 2 evolution, a 40 g sieved soil sample was put into a 50 mL perforated plastic beaker and was bought to a moisture level equivalent to field capacity solely through capillary action and the CO 2 released was determined by Licor 840 CO 2 detector [38]. All the laboratory procedures followed to determine soil health indicators can be found in Haney, Haney, Smith, Harmel, and White [11]. From these soil indicator values; SHI was derived as Equation (1). The SHI includes the weighted contribution of microbial activity, water-extractable organic C, and water-extractable organic N. The SHI score represents the overall health of the soil system and its value ranges between 0 and 50 or slightly over 50 in most agriculture systems, wherein scores > 7 are considered relatively good.
where one day CO 2 is the amount of CO 2 released in 24 h, WEOC as water-extractable organic C, and WEON as water-extractable organic N in mg/kg, the numbers 10 and 100 are the weighting factors.

Field Apparent Electrical Conductivity Survey and Mapping
A field survey was conducted in 2020 to collect ECa data across the field using a DualEM-1S sensor (DUALEM Products. Available online: https://dualem.com/products/ (accessed on 16 April 2022)). The DualEM-1S sensor works with the principles of electromagnetic induction. It has one transmitter (vertical dipole) and two receivers (vertical and horizontal dipole) separated by a distance of about 1 m. About 70% of the cumulative response of the horizontal configuration of the receiver comes from 0.5 m soil depth, whereas it corresponds to a depth of 1.5 m for the vertical configuration [39]. During the field operation, the sensor was tightly placed in a hard-core plastic sled and was dragged in the field with an all-terrain vehicle. A global positioning system (GPS Pathfinder Pro XRS) was connected to the sensor to record geographic coordinates in real-time. The sensor received power from the all-terrain vehicle and was connected to a field computer Yuma 2 Rugged Tablet. The communication between the sensor, GPS unit, and the field computer was acquired by running a Termite Software Program (Termite: a simple RS232 terminal. Available online: https://www.compuphase.com/software_termite.html (accessed on 16 April 2022))and its Geographic Information System (GIS) display by using Handheld GIS Mapping Software (HGIS GPS mapping. Available online: https://www.atsgps.com/hgis-gps-mapping.html (accessed on 16 April 2022)). The field survey was done in a single scan, and it took about 4.5 h to complete. There was a total of 85 drive-paths each placed roughly 5 to 8 m apart along the northeast-southwest direction of the field. ECa data were recorded every second that corresponded to a measuring distance of about 2.5 to 3 m along the drive path, making a total of 14,171 observations covering the entire study area. Each ECa measurement represented a GPS-recorded electrical conductivity in dS/m acquired for horizontal as well as the vertical configuration of the receiver; however, the present study only used the ECa data from its horizontal configuration, as it is more responsive to surface soil depth (~0.5 m) and properties.
Once the field scan was complete, recorded data in txt format were downloaded and imported to Microsoft Excel and to ArcGIS [40] for data cleaning and subsequent analysis. Data cleaning included removing data points falling outside the field boundary, overlapping points, out-of-drive-track points, negative points, and possible outliers determined as values falling outside this range [ECa value < (mean − 3 × standard deviations (SD)); > (mean + 3 × SD)]. A continuous ECa map at a grid size of 5 m × 5 m was generated by ordinary kriging as it outperformed inverse distance weighting, radial bias functions, and empirical Bayesian kriging methods in terms of prediction error (root mean squared error) based on leave-one-out cross-validation. A grid size of 5 m was selected for mapping to make the dimensions of the smallest spatial unit (the grid) compatible with the swath length of most farm equipment used in the field. Before the kriging operation, spatial autocorrelation of ECa observations (n = 14,171) was modelled with a variogram (nugget: 0.0, partial sill: 0.07, range: 320 m, model: spherical) considering 12 lags with a lag distance of 30 m in VESPER program [41].

Topographic Information
A digital elevation model (DEM) was compiled from light detection and ranging technology and the terrain attributes derived from the DEM were used as a source of topographic information in the study area. The DEM was downloaded from the USDA-Natural Resources Conservation Services, Geospatial Data Gateway website (Geospatial data gateway. Available online: ( https://gdg.sc.egov.usda.gov/GDGHome.aspx (accessed on 5 December 2021)). The DEM was originally compiled at a grid resolution of 1 m × 1 m, but it was resampled to a 5 m × 5 m grid for this study by applying the bilinear interpolation technique. Before deriving terrain attributes, the DEM was pre-processed by filling unnecessary sinks which would otherwise disrupt regular water flow and distribution on the soil surface. A total of 15 terrain attributes were derived from the DEM in SAGA GIS [42] platform and were used as predictors of soil health indicators, and SHI in the study area. The terrain attributes used were aspect (Aspect), elevation (Elevation), flow accumulation (FlowAccu), landscape position (Landforms), mid-slope position (Midslp), minimum curvature (Curvmin), multiresolution ridge top flatness (MRRTF), multiresolution valley bottom flatness (MRVBF), overland flow distance (Overlflw), wetness index (SAGAWI), slope gradient (Slope), slope height (Slopeht), slope-length factor (LSFactor), topographic position index (TPI), and valley depth (Valdep). More detail on these terrain attributes and their derivation can be found in Guo, et al. [43]. Table 1 lists all the terrain attributes and ECa and their correlation with one-day CO 2 , organic C, organic N, and SHI from the study area. Table 1. Pearson's correlation coefficient between measured soil properties, apparent electrical conductivity, and terrain attributes.

Assessment of Autocorrelation and Spatial Dependency
Spatial autocorrelation of one-day CO 2 , organic C, organic N, SHI, and ECa was investigated using variograms that plot semivariance among observation pairs against their separation distances (Equation (2)). It measures the average degree of dissimilarity between unsampled values and a nearby data value [44] and thus can describe autocorrelation among observations at various distances. We calculated omnidirectional experimental variograms of measured properties following Equation (2) and modeled its spatial characteristics using three different variogram models (Spherical model, Exponential model, and Gaussian model). The variogram parameters nugget (C 0 ), partial sill (C 1 ), and range of the variogram (a) from the best-fitted model among the three tested that provided the lowest value of Akaike Information Criteria were recorded and were used to characterize spatial autocorrelation or spatial dependence of the measured properties using nugget: sill ratio (NSR). The NSR represents a contribution of the nugget to the overall spatial structure of the variogram and it can be calculated as NSR = C 0 /(C 0 + C 1 ). The NSR tells how well the measured property is spatially dependent or autocorrelated (NSR < 0.25, strong; 0.25 < NSR > 0.75, moderate; and NSR > 0.75, weak) [45]. The experimental variograms were calculated considering 30 lags with a lag distance of 20 m for all measured properties. Lags are the distances between observation pairs at which the variogram is calculated. Variogram calculation and model fitting was performed in the VESPER program [41].
where, γ (h) is the variogram for a distance h(lag) between observations z(x a ) and z(x a + h) with N (h) the number of observation pairs separated by h.

Spatial Prediction, Mapping, and Model Validation
Prediction of one-day CO 2 , organic C, organic N, SHI, and ECa and its mapping across the study area was performed with a machine learning technique-Random Forest. The RF model was employed using a package 'randomForest' in the R environment [46], with three different parameters defined in the model: the number of trees (ntree) set to be 500, the minimum number of data points be used in each terminal node (nodesize) set to be 1, and the number of variables tried at each node (mtry) was 4 determined with 'tuneRF' function. The model also quantified the relative importance (RI) of predicting variables based on the ranking of percentage increase in prediction error (mean squared error) (%IncMSE) on its removal from the model. Before executing the RF model, point observations were randomly split into training (75% of the data) for model calibration and test datasets (25% of the data) for model performance evaluation. The mean and SD values of the measured properties in training and test datasets were comparable. A correlation matrix of one day CO 2 , organic C, organic N, and SHI and predictor variables was generated and the variables that were at least significant at p value of <0.01 for each property were used in the RF models for subsequent predictions. Once the predictions were made, the maps were exported to ArcGIS [40] for necessary map edits and layout.
The performance of the prediction models was evaluated on both the calibration as well as the validation data sets using the statistical indices of R 2 , Lin's Concordance Correlation Coefficient (LCCC), root mean squared error (RMSE), and mean error (ME). The LCCC was calculated as Equation (3), the calculation of the remaining indices could be found elsewhere, e.g., [47].
where X and Y are the mean of the observed (X) and predicted values (Y), ρ is the correlation coefficient between observations and predictions, σ X and σ Y as corresponding standard deviations, and σ 2 X and σ 2 Y corresponding variances.

Correlation between ECa and Soil Health Index and Its Mapping
Spatial correlation between ECa and SHI maps was quantified using a local correlation between two maps [48], and the values were mapped across the study area at 5 m spatial resolution. Around each pixel of the raster, a 3 × 3 pixel neighborhood was defined to record local correlation using the 'cor' function, and the operation was repeated throughout the raster using the 'focal' function in the R environment [46]. Further, an empirical relationship between ECa and SHI was established using a linear model where SHI was used as a dependent variable on ECa for each predicted pixel from the map. The linear equation could then be used to predict SHI from easy to obtain ECa data.

Descriptive Summary of Observation Data
Summary statistics of one-day CO 2 , organic C, organic N, SHI, and ECa for the measurement locations are listed in Table 2. One day CO 2 measurements ranged from 9.9 to 123 mg/kg with a mean of 42.7 (±17.5 SD). Average organic C and organic N were 154 and 11.7 mg/kg, respectively, with their corresponding coefficient of variation (CV) of 18.9, and 24%. Similarly, SHI ranged from 2.6 to 16.6 with a mean of 8.5 (±2.2 SD), and that of ECa from 0.5 to 1.5 dS/m with a mean of 1.0 ((±0.2 SD). Among all properties, one-day CO 2 had the highest CV (41%) followed by SHI, and organic C had the smallest CV (18.9%) of all. All properties data had positive skewness coefficients except for organic N and ECa which were slightly negatively skewed with skewness coefficients of −0.4, and −0.2, respectively.

Spatial Autocorrelation
The experimental variogram of the measured properties was best modeled with a spherical variogram as the model showed the lowest value of Akaike Information Criteria among the tested models-Spherical, Exponential, and Gaussian models (Figure 2). The range of the variogram of one-day CO 2 , organic C, and organic N was between 360-370 m, and NSR between 0.25 and 0.41 showing a strong to moderate spatial dependence of these properties. Similarly, a strong spatial dependence was observed for SHI, which had a range of 396 m, and a partial sill of 4.5. Among these four properties, SHI had the lowest nugget variance and organic C had the highest nugget of all. ECa was strongly spatially dependent with a nugget variance of zero, and a range of 383 m.

Spatial Autocorrelation
The experimental variogram of the measured properties was best modeled with a spherical variogram as the model showed the lowest value of Akaike Information Criteria among the tested models -Spherical, Exponential, and Gaussian models (Figure 2). The range of the variogram of one-day CO2, organic C, and organic N was between 360-370 m, and NSR between 0.25 and 0. 41 showing a strong to moderate spatial dependence of these properties. Similarly, a strong spatial dependence was observed for SHI, which had a range of 396 m, and a partial sill of 4.5. Among these four properties, SHI had the lowest nugget variance and organic C had the highest nugget of all. ECa was strongly spatially dependent with a nugget variance of zero, and a range of 383 m.

Correlation of Soil Properties with Predictors
A Pearson's correlation coefficient of one-day CO 2 , organic C, organic N, and SHI with predicting variables that included terrain attributes and ECa is listed in Table 1 with their corresponding level of significance. Among the 16 predictors, 14 predictors were significantly correlated (p-value < 0.01) with one-day CO 2 , 10 predictors with organic C, 9 predictors with organic N, and 13 predictors with SHI, respectively. The predictors that were not significantly correlated with one-day CO 2 were MRRTF and Overlflw. Similarly, predictors that were not significant with organic C were Aspect, ECa, FlowAccu, Midslp, Curvmin, and Overlflw. The same predictors that were not significant with organic C except FlowAccu were also non-significant for organic N. Two additional predictors namely, Slope and LSFactor were found non-significant to predict organic N as well; however, both were significantly correlated with organic C. In the case of SHI, all predictors except Aspect, MRRTF, and Overlflw were significantly correlated.

Model Performance and Predicted Maps
The RF model used to predict and map soil health indicators and SHI was evaluated on 25% set-aside data that were not used in model calibration, and Table 3 lists model evaluation results for both training and test datasets. The predicted R 2 value ranged from 0.85 to 0.90 for training data, and from 0.24 to 0.48, respectively, for the test data; one-day CO 2 had the highest R 2 , and organic C had the lowest R 2 of all. The model produced the highest LCCC for one-day CO 2 and the lowest RMSE value for SHI, and the lowest LCCC and the highest RMSE for the prediction of organic C in the test data set. The ME statistics showed that the model was negatively biased to predict one-day CO 2 . Overall, the RF model showed that the SHI was best predicted with the higher value of R 2 and LCCC and the lower value of RMSE, and organic C was least predicted with the lowest value of R 2 and LCCC and the highest value of RMSE (Table 3). The predicted map of one-day CO 2 , organic C, organic N, and SHI are shown in Figure 3, and their corresponding important predictors are shown in Figure 4, respectively. These maps were generated by using the RF model in which only the significant predictors of each property (Table 1) were used in their subsequent predictions. Based on the predicted maps, the investigated field seems variable in terms of soil health indicators and SHI distribution; CVs of these properties ranged between 8% and 20%. Among all properties, organic C distribution in the field was least variable (mean = 157 mg/kg; CV = 8.2%) and organic N was moderately variable (mean = 11.7 mg/kg; CV = 13%). The maps showed that a major portion of the central part of the field had a lower level of one-day CO 2 and organic C content resulting in a lower SHI value in that area compared to the rest of the field. On the other hand, most of the southern part of the field had higher values of those properties including organic N. Specific to organic N which was higher in the lower 3 quarters of the field, its maximum value was mapped from the southernmost part, and in the western corner of the field. These maps also showed the influence of field terracing on the spatial distribution of measured properties across the field (Figure 3). Agronomy 2022, 12, x FOR PEER REVIEW 9 of 15

Important Predictors for Soil Health
Among the 16 predictors (Table 2), SAGAWI was ranked in the list of top five predictors of all properties indicating its higher influence on the spatial distribution of these properties in the study area ( Figure 4). The RI of SAGAWI ranged between 13% for one-day CO2 to 29% for organic C and organic N prediction. Similarly, ECa was among the top three predictors of one-day CO2 and SHI for which the highest RI of nearly 21% was observed. However, ECa was not considered a good predictor of organic C and organic N distribution, even though it had some correlation with those properties. Predictors such as Curvmin, slope, FlowAccu, and MRRTF were less important compared to multiresolution valley bottom flatness, TPI, Landforms, and slope height to predict soil health indicators and SHI. Predictors like multiresolution valley bottom flatness, SAGAWI, and TPI were the top terrain attributes to predict one-day CO2, organic C, organic N, and SHI in the study area. Overall, the spatial distribution of soil health indicators and SHI in the study area was greatly influenced by ECa and local topography and their level of influence depended on the properties considered ( Figure 4).

Important Predictors for Soil Health
Among the 16 predictors (Table 2), SAGAWI was ranked in the list of top five predictors of all properties indicating its higher influence on the spatial distribution of these properties in the study area (Figure 4). The RI of SAGAWI ranged between 13% for one-day CO 2 to 29% for organic C and organic N prediction. Similarly, ECa was among the top three predictors of one-day CO 2 and SHI for which the highest RI of nearly 21% was observed. However, ECa was not considered a good predictor of organic C and organic N distribution, even though it had some correlation with those properties. Predictors such as Curvmin, slope, FlowAccu, and MRRTF were less important compared to multiresolution valley bottom flatness, TPI, Landforms, and slope height to predict soil health indicators and SHI. Predictors like multiresolution valley bottom flatness, SAGAWI, and TPI were the top terrain attributes to predict one-day CO 2 , organic C, organic N, and SHI in the study area. Overall, the spatial distribution of soil health indicators and SHI in the study area was greatly influenced by ECa and local topography and their level of influence depended on the properties considered ( Figure 4).

Relationship between ECa and soil Health Index
Spatial correlation between SHI and ECa maps based on each predicted pixel quantified mapped as in Figure 5A and the linear relationship between them is shown in 5B, respectively. T was a correlation of about 0.55 between ECa and SHI based on all pixels and its distribution sh some patchy patterns across the field. The western part of the field showed a lower or neg correlation between ECa and SHI compared to the rest of the field; the central and the eastern of the field showed a mixed low and high correlation between ECa and SHI. Based on the l relationship between SHI and ECa, 30% of SHI variability was explained by ECa only, and th former could be predicted by ECa with an equation: y = 6.21 + 2.71*x.

Relationship between ECa and soil Health Index
Spatial correlation between SHI and ECa maps based on each predicted pixel quantified and mapped as in Figure 5A and the linear relationship between them is shown in 5B, respectively. There was a correlation of about 0.55 between ECa and SHI based on all pixels and its distribution showed some patchy patterns across the field. The western part of the field showed a lower or negative correlation between ECa and SHI compared to the rest of the field; the central and the eastern part of the field showed a mixed low and high correlation between ECa and SHI. Based on the linear relationship between SHI and ECa, 30% of SHI variability was explained by ECa only, and that the former could be predicted by ECa with an equation: y = 6.21 + 2.71x. Agronomy 2022, 12, x FOR PEER REVIEW 11 of 15 Figure 5. (A) Local correlation map between apparent electrical conductivity (ECa) and soil health index, and (B) a linear relationship between the two. Data points in 5B correspond to all the predicted pixels at 5m grid resolution, continuous line with the shaded area are a linear fit and the standard error, and the green eclipse shows a 95% prediction eclipse.

Discussion
In this study, a 35 m × 35 m grid design was applied to identify sampling locations for mapping within-field soil health status across a field. As reported in Adhikari, Smith, Collins, Haney, and Wolfe [35], the selected grid size of 35 m × 35 m and the number of samples collected were representative of the study area to capture field variability, as the grid size was smaller than the range of the modeled variograms of measured soil properties [49]. Furthermore, grid sampling is a widely used soil sampling scheme in field-scale soil attribute mapping, e.g., [50,51] because of its simplicity and potential for greater efficiency [52]. Grid samplings are useful to access within-field variations in plant-available nutrients and facilitate site-specific management [53,54]. Spatial autocorrelation of soil health indicators and SHI including ECa was investigated with geostatistical methods (i.e., variogram) which are common in spatial studies in soil and agronomic research [55]. As suggested by [56], the spherical model is a common variogram model to describe the autocorrelation of soil properties, and our study reported the same as the experimental variogram of all properties were best modeled by a spherical model. We observed that one-day CO2, SHI, and ECa had a strong spatial dependence, and organic C and organic N had a moderate dependence in the study area. A strong spatial dependence of SHI at a field scale was previously reported by Amirinejad, et al. [57] in which the variogram range of organic C and SHI varied between 300-380 m, similar to the variogram ranges observed in this study. The range values of one-day CO2, organic C, and organic N were between 360-370 m, and it is the average distance to which those properties are correlated in the field [58].
Spatial variability in soil properties affects crop performance, yield, and adoption of site-specific management practices. The study site showed a moderate to a high degree of variability (i.e., CV 15-35%: moderate; CV > 35%: high [59] as the CV of measured soil properties including SHI, and ECa ranged between 18.9 to 41% ( Table 2). The spatial variability of those properties was predicted and mapped with a RF model which is becoming popular in field-scale soil properties mapping because of its higher prediction performance, among other benefits. Studies by Schmidt, et al. [60] and Zhang, Ji, Saurette, Easher, Li, Shi, Adamchuk, and Biswas [32] found a higher and more consistent performance of the RF model compared to regression tree and generalized linear model in predicting soil properties at field level; our results are consistent with these findings as we observed a similar model performance (R 2 as high as 0.90). However, a recent study from Denmark found the performance of kriging (R 2 = 0.90) as high as the RF model (R 2 = 0.91) in predicting fieldscale variability of soil organic matter [50]. Authors argued that the higher performance of kriging might be due to a high sampling density (20 m × 20 m grid sampling scheme) in their field.
The present study used a range of terrain attributes and ECa which were collected on-the-go as predictors of one-day CO2, organic C, organic N, and SHI. Research showed that collecting ECa data on the go is an efficient way of capturing field variability as ECa is mostly related to soil physicochemical properties (e.g., porosity, soil moisture, salinity, cation composition, cation exchange

Discussion
In this study, a 35 m × 35 m grid design was applied to identify sampling locations for mapping within-field soil health status across a field. As reported in Adhikari, Smith, Collins, Haney, and Wolfe [35], the selected grid size of 35 m × 35 m and the number of samples collected were representative of the study area to capture field variability, as the grid size was smaller than the range of the modeled variograms of measured soil properties [49]. Furthermore, grid sampling is a widely used soil sampling scheme in fieldscale soil attribute mapping, e.g., [50,51] because of its simplicity and potential for greater efficiency [52]. Grid samplings are useful to access within-field variations in plant-available nutrients and facilitate site-specific management [53,54]. Spatial autocorrelation of soil health indicators and SHI including ECa was investigated with geostatistical methods (i.e., variogram) which are common in spatial studies in soil and agronomic research [55]. As suggested by [56], the spherical model is a common variogram model to describe the autocorrelation of soil properties, and our study reported the same as the experimental variogram of all properties were best modeled by a spherical model. We observed that one-day CO 2 , SHI, and ECa had a strong spatial dependence, and organic C and organic N had a moderate dependence in the study area. A strong spatial dependence of SHI at a field scale was previously reported by Amirinejad, et al. [57] in which the variogram range of organic C and SHI varied between 300-380 m, similar to the variogram ranges observed in this study. The range values of one-day CO 2 , organic C, and organic N were between 360-370 m, and it is the average distance to which those properties are correlated in the field [58].
Spatial variability in soil properties affects crop performance, yield, and adoption of site-specific management practices. The study site showed a moderate to a high degree of variability (i.e., CV 15-35%: moderate; CV > 35%: high [59] as the CV of measured soil properties including SHI, and ECa ranged between 18.9 to 41% ( Table 2). The spatial variability of those properties was predicted and mapped with a RF model which is becoming popular in field-scale soil properties mapping because of its higher prediction performance, among other benefits. Studies by Schmidt, et al. [60] and Zhang, Ji, Saurette, Easher, Li, Shi, Adamchuk, and Biswas [32] found a higher and more consistent performance of the RF model compared to regression tree and generalized linear model in predicting soil properties at field level; our results are consistent with these findings as we observed a similar model performance (R 2 as high as 0.90). However, a recent study from Denmark found the performance of kriging (R 2 = 0.90) as high as the RF model (R 2 = 0.91) in predicting field-scale variability of soil organic matter [50]. Authors argued that the higher performance of kriging might be due to a high sampling density (20 m × 20 m grid sampling scheme) in their field.
The present study used a range of terrain attributes and ECa which were collected on-the-go as predictors of one-day CO 2 , organic C, organic N, and SHI. Research showed that collecting ECa data on the go is an efficient way of capturing field variability as ECa is mostly related to soil physicochemical properties (e.g., porosity, soil moisture, salinity, cation composition, cation exchange capacity, ionic strength, particle size distribution, particle shape and orientation, wettability) [61], and it can be collected intensively in an easy and inexpensive way [62,63]. Due to their correlation with soil properties, ECa and terrain attributes are widely used in spatial prediction and mapping of soils properties at the field scale, e.g., [51,[64][65][66][67]. Among the predictors used, we found a weak correlation of ECa with organic C, and organic N, whereas it was highly correlated with SHI and one-day CO 2 . Our study highlighted the importance of terrain attributes such as SAGAWI, MRVBF, and TPI in predicting soil health indicators in the study area. SAGAWI had the highest influence on Organic C, Organic N, and SHI (RI > 25%) suggesting that the field areas that tend to remain moist or low-lying areas had a strong and positive (Table 1) influence on these properties. A positive influence of wetness index on soil carbon distribution was also reported in other studies, e.g., [50,68]. As low-lying or concave areas in the field tend to accumulate organic matter due to erosion from upslope areas and favor slow decomposition of organic matter [69,70], terrain attributes like SAGAWI, MRVBF, and Valdep have a high potential to predict soil carbon and nitrogen distribution in those areas. We also observed the influence of field terracing on soil properties distribution (Figure 3). Those terraces were built to reduce soil erosion, and increase infiltration 30-40 year ago, which might have influenced field hydrology and soil and nutrient transport [71]. The terrain attributes we used for prediction were able to identify those terracing signatures and perhaps captured their influence on soil properties variations across the field.
The correlation map between ECa and SHI shows a pixel-by-pixel relationship between these properties and it could provide a lot more information in the spatial context. This could be very useful in within-field soil variability assessment and for site-specific management decisions as it shows a direct relationship between these properties from every paddock of the field. Establishing an empirical relationship between ECa and SHI is helpful, as it highlights the potential of on-the-go sensing to estimate SHI in the field easily and efficiently with a minimum cost.

Conclusions
This study was conducted to quantify and map within-field variability of the soil health status of Blackland prairies soils in Texas. Two hundred and eighteen topsoil samples were collected and analyzed for one-day CO 2 , organic C, organic N, and soil health index. A High-resolution (5 m × 5 m grid) map of each property was generated using a random forest model in which terrain attributes and ECa data were used as predictors. Furthermore, the spatial correlation between ECa and SHI was mapped and a linear relationship between them was established. Results showed that one-day CO 2 , organic C, and organic N, SHI, and ECa data were moderate to strongly spatial dependent in the study area, and the range of this dependency was between 259 to 396 m. The ECa was strongly positively correlated with one-day CO 2 , and SHI, whereas it was weakly correlated with organic C and organic N. There was a strong positive correlation between one-day CO 2 , organic C, and organic N, and SHI and SAGAWI, and a strong negative correlation with TPI. Based on the predicted maps, average values of one-day CO 2 , and SHI were 47.1 mg/kg (±7.9 SD), and 8.9 (±1.2 SD), respectively. The field shows high variability in terms of SHI. The ECa, SAGAWI, MRVBF, and TPI were among the top predictors of soil health indicators, and SHI in the study area, while ECa was a good predictor of one-day CO 2 , and SHI, but not for organic C and organic N. Validation of the RF model showed a higher prediction performance for one-day CO 2 and SHI but a lower performance for organic C. The ECa and SHI maps were correlated (r = 0.55), and the latter could be fairly estimated with ECa with a linear function y = 6.21 + 2.71x. There was a noticeable variation in the selected soil health indicators, and SHI in the study area at within-field scale, and such variation could play a role in yield variations across the field. We recommend using topography and ECa to predict soil health indicators and soil health indexes. Further, SHI could be fairly estimated with ECa data that could be collected on-the-go efficiently with minimum cost; however, it needs to be tested in other fields as well. We believe the high-resolution predicted maps of soil health indicators, and SHI including the correlation map between ECa and SHI could be useful in fine-tuning SSCM applications.