Feasibility assessment on use of proximal geophysical sensors to support precision management

A study was conducted at three sites in North Dakota to strengthen understanding of the usefulness of different proximal geophysical data types in agricultural contexts of varying pedology. This study hypothesizes that electromagnetic induction (EMI), gamma‐ray sensor (GRS), cosmic‐ray neutron sensor (CRNS), and elevation data layers are all useful in multiple linear regression (MLR) predictions of soil properties that meet expert criteria at three agricultural sites. In addition to geophysical data collection with vehicle‐mounted sensors, 15 soil samples were collected at each site and analyzed for nine soil properties of interest. A set of model training data was compiled by pairing the sampled soil property measurements with the nearest geophysical data. Eleven models passed expert‐defined uncertainty criteria at Site 1, 16 passed at Site 2, and 14 passed at Site 3. Electrical conductivity (EC), organic matter (OM), available water holding capacity, silt, and clay were predicted at Site 1 with an R‐squared of prediction (Rpred2)$(R_{pred}^2)$ > .50 and acceptable root mean square error of prediction (RMSEP). Bulk density (BD), OM, available water capacity, silt, and clay were predicted with Rpred2$R_{pred}^2$ > .50 and acceptable RMSEP at Site 2. At Site 3, no soil properties were predicted with acceptable RMSEP and an Rpred2$R_{pred}^2$ > .50. These results confirm feasibility of our method, and the authors recommend the prioritization of EMI data collection if geophysical data collection is limited to a single mapping effort and calibration soil samples are few.


INTRODUCTION
Predictive soil mapping with proximal geophysical data has potential to benefit precision agriculture because proximal sensors such as the gamma-ray sensor (GRS) (IAEA, 2003;van der Veeke et al., 2021), the cosmic-ray neutron sensor (CRNS) (Desilets et al., 2010;Zreda et al., 2008), and electromagnetic induction (EMI) (Abdu et al., 2008;Gibson & Franz, 2018) have footprint sizes that can characterize soil on the subfield scale. In the United States, field sizes vary widely depending on region and crop type, but the median field size is 27.8 ha (Yan & Roy, 2016), and management on the subfield scale (around 0.4 ha) is possible because of recent advancements in fertilizer, planter, sprayer, and irrigation equipment (Hamrita et al., 2000;O'Shaughnessy et al., 2019). Precision agriculture manages inputs such as water, fertilizer, and seeding rate and variety on a subfield scale to maximize profit, which often means maximizing yield while optimizing the timing and placement of input resources. Soil texture, pH, available water capacity, cation exchange capacity (CEC), electrical conductivity (EC), and organic matter (OM) content are all related to setting and obtaining yield goals and considered valuable information for precision management decisions (Shearer & Ward, 1999). Possible variable management responses to soil maps that portray subfield variation include irrigation, seeding rate, tillage settings, liming, and application of compost and fertilizer (van Egmond et al., 2010). This paper aims to determine if predictive soil maps useful for variable management can be created with a combination of proximal geophysical data sources and in-field soil sampling.
Data fusion, or using multiple data sources as predictive data, is a common approach for predicting soil properties. Combinations explored in the past have included EMI, GRS, elevation, and visible and near-infrared data. Rodrigues et al. (2015) used EMI and GRS and found that, using principal components, regression models of clay and cation exchange capacity were significant (p < .05) at five out of eight study sites. Castrignanò et al. (2012) found that different soils in Western Australia that produced similar responses in a single sensor (sandy, sandy gravelly, sandy salt-affected, and clayey soils) could be discriminated when using combined EMI, GRS, and elevation. Additionally, Castrignanò et al. found correlation (r ≥ .46) between GRS and soil organic carbon (SOC), plant-available K, and P, and found weaker correlation (r ≤ .31) between EMI and P and pH. Elevation was correlated with SOC, plant-available K, and P with correlation coefficients of .28 to .39. In another study, Ji et al. (2019) could not predict extractable K and P with combined information from elevation, GRS, EMI, and visible and nearinfrared data. However partial least-squares regressions of soil OM, pH, lime buffering capacity, Ca, Mg, and Al were usu-

Core Ideas
• Framework given to assess relationships between proximal sensing and soil properties. • Correlations between sensing data and soil properties varied among three study sites. • Electromagnetic induction was most consistently useful in soil property prediction.
ally improved by substituting the data fusion approach for a single sensor, obtaining R 2 > .5.
The present study combines EMI, GRS, CRNS, and elevation data. Each of these measurements is theoretically related to a variety of soil characteristics. The EMI measures apparent bulk electrical conductivity (ECa), which is affected by soil water content, soil temperature, clay content, mineralogy, bulk density (BD), and salinity McBratney et al., 2005). Gamma-ray sensors detect naturally emitted gamma radiation from K-40 and the gamma-rays emitted by the U-238 and Th-232 decay series. Detected gamma radiation is influenced by soil water content, parent material mineralogy, OM, and texture (Carroll, 1981;Dierke & Werban, 2013). The CRNS measures low-energy neutron counts (NC) (∼0.25-1,000 eV), which are an established method for soil water estimation (Zreda et al., 2008). Low-energy NC may also serve as a proxy for overall soil variability related to properties such as OM content or available water capacity (Andreasen et al., 2017;Finkenbiner et al., 2019). Elevation is connected to soil formation, soil water, and organic carbon content (Florinsky et al., 2002). The predictive data suite comprised of EMI, GRS, CRNS, and elevation data allows this study to freshly examine any possible correlations between these proximal sensing data and soil properties of interest.
Soil mapping with data fusion is intriguing not only because of the wide range of potential data combinations, but also because the relative performances of different predictive data types vary across settings. Wong and Harper (1999) concluded that the usefulness of gamma-ray spectroscopy alone is limited because relationships between K-40 counts and soil properties did not hold everywhere for sites in Western Australia. This suggests that site specific calibrations are required for soil property predictions with gamma-ray data and that inclusion of additional sensors would be informative. For instance, in Queensland and South Australia, Rodrigues et al. (2015) found that predictions of clay and cation exchange capacity were improved by using principal components from EMI and gamma-ray data as predictors vs. predicting with EMI or gamma-ray data alone. However, the geophysical information most strongly correlated with a given soil property differed between field sites of varying pedology and geographic location. Rodrigues et al. also explored universal calibration for Australian soils by combining sensor and soil sample data from all their sites and found adjusted R 2 values of .27 and .22 for predictions of CEC and clay, respectively. The varying results of data fusion in different settings mean that each analysis of a new site adds valuable information to our understanding of which geophysical sensor is most crucial in given situations. Because the GRS, CRNS, and EMI explore different wavelengths on the electromagnetic spectrum, we expect each sensor to obtain novel information. Similar to use of visible and near-infrared bands to calculate the normalized difference vegetation index (NDVI), this paper pursues integration of GRS, CRNS, EMI, and elevation into new information that characterizes the field.
Understanding the proper situations for different sensors will inform producers and researchers as they navigate the numerous commercial soil mapping technologies available. At least one soil mapping company, SoilOptix (Canada), has arisen that provides gamma-ray mapping technology and support, using sensors produced by Medusa Radiometrics (Netherlands). The EMI and direct current resistivity are standard soil mapping capabilities offered by numerous companies. The CRNS is still an emerging technology, but the sensor is commercially available through several companies such as Hydroinnova LLC. Given the current accessibility of commercial EMI, GRS and CRNS surveys, determining the predictive ability of these tools in new agricultural contexts is extremely timely.
Predictive soil mapping methods in the literature include support vector machine, random forest, classification and regression trees, partial least squares regression (PLSR), bagging-PLSR, multivariate adaptive regression splines, K nearest neighbor, and co-kriging (Castrignanò et al., 2012;Ji et al., 2019;Piikki et al., 2013;Rossel et al., 2007;Söderstrom et al., 2016). In addition to the more complex modeling approaches, multiple linear regression (MLR) also has extensive precedent due to its simplicity (Mahmood et al., 2013;van der Klooster et al., 2011;van Egmond et al., 2010) and high interpretability. In our analysis we use MLR because it is pragmatic given the expected and desired small soil sample sizes usually attainable by producers and crop consultants.
The cost and time required for soil core sampling limits methods in both precision agriculture and other aspects of the agriculture industry. Another sector needing maximal information return on few soil samples is the monitoring, verification, and reporting of soil organic carbon (SOC). The SOC market has the potential to be a viable source of income for producers, but the system is limited by poor information on producers' actual SOC storage. Third party companies verify the C credits that farmers sell, and this verification service comprises roughly 75% of the total cost of producing C credits (Plume, 2021). This study addresses how well SOC (or OM here) can be predicted from limited samples with the help of geophysical surveys, and which geophysical data types are preferred.
The hypothesis of this study is that EMI, GRS, CRNS, and elevation data layers are all useful in MLR predictions of soil properties that meet expert criteria at three agricultural sites in North Dakota. The objectives are to: (a) determine soil property predictions for BD, texture (percentage sand, silt, and clay), available water capacity, and OM that meet validation criteria at each site; (b) recommend which predictive geophysical data type among EMI, GRS, CRNS, and elevation is expected to produce successful MLR predictions most often; and (c) evaluate feasibility of using data fusion and MLR with small sample size for soil property prediction in precision agriculture.

Study sites
Each of the three sites considered in this study is a roughly 53 ha agricultural field located in Southeast North Dakota. The sites were selected based on the following criteria: (a) they had an existing or prior USDA-NRCS Environmental Quality Incentives Program contract on variable rate irrigation, (b) the water table was below the crop rooting zone, and (c) the sites were relatively close to Fargo, ND. Average annual rainfall for the region during the time period 2007-2020 is 448 mm and annual potential evapotranspiration is 1,262 mm (North Dakota Agricultural Weather Network, 2021; Lisbon station). For southeastern North Dakota, the normal monthly low temperature ranges from -17.9˚C (January) to 14.7˚C (July), and the normal monthly high temperature ranges from -7.5˚C (January) to 27.8˚C (July). Normal temperatures reported here are from the period 1991-2020 as reported by the National Weather Service. All fields are center pivot irrigated. The crops grown at Sites 1, 2, and 3 in 2020 were soybean [Glycine max (L.) Merr.], potato (Solanum tuberosum L.), and maize (Zea mays L.), respectively, with growing seasons from May through September. Southeastern North Dakota generally experiences its first killing frost around 1 October and the soil is free of frost again around 1 April. At Site 1, wetlands fill depressions on the East (7.5 ha) and South (6.7 ha) sides of the field. Site 2 has a moraine feature in the Southeast corner of the field with a maximum height of 11 m above the rest of the field. A shallow depression is oriented West-Northwest through the middle of Site 3. In Southeast North Dakota, the surface geology is a patchwork of till, glacial outwash, deltaic deposits, glacial lacustrine sediment, and aeolian sand (Bluemle, 1975). Locations of the field sites amid the variable surface geology are F I G U R E 1 Surficial geology of Southeast North Dakota (units grouped by primary lithology type). Locations of Sites 1, 2, and 3 are plotted with black triangles shown in Figure 1. Site 1 lies on a Holocene aeolian sand deposit about 1 km South of the Sheyenne River. Site 2 is on a Holocene glacial outwash deposit of cross-bedded sand and plane-bedded gravel. Site 3 also sits on a Holocene deposit of bedded sand and gravel, about a 0.5 km Northeast of the modern James River (State of North Dakota, NDGISHUB Surface Geology). The unconsolidated sediments at all three sites are underlain by Cretaceous calcareous shale: the Greenhorn Formation at Site 1 and the Niobrara Formation at Sites 2 and 3 (State of North Dakota, NDGISHUB Bedrock Geology). Generally, the soil types at all the sites are loams or sandy loams. Figure 2 depicts the soil series present in each field in greater detail.

Geophysical data collection and processing
Geophysical surveys were performed on 15 Oct. 2020 at Sites 1 and 2 and on 16 Oct. 2020 at Site 3. The EMI, CRNS, and GRS data were simultaneously collected from a vehicle traveling approximately 10-15 km h -1 in transects spaced roughly 10 m apart.
The EMI was performed with a DUALEM-21 sensor, pulled in a plastic sled behind the vehicle, to obtain ECa data in mS m -1 . Shallow apparent bulk electrical conductivity (ECaS), deep apparent bulk electrical conductivity (ECaD), and ratio of shallow to deep electrical conductivity (ECaSDR) were recorded. Only the 2-m coil spacing array was used. The horizontal co-planar coil orientation penetrates the surface to about 1-m depth (ECaS), and the perpendicular coil orientation penetrates to roughly 2.5-to-3-m depth (ECaD) (Dualem Inc., 2013). Soundings were recorded every second, and the location of each measurement was recorded with a Hemisphere GPS XF101 DGPS (Juniper Systems, Inc.) unit.
Outliers and redundant data were removed from the raw ECa data to assure basic quality.
A passive, vehicle-mounted cosmic-ray neutron detector (eight ∼1.8 m CRS 2000/B tube capsules from Hydroinnova LLC) recorded accumulated NC in 1-min intervals (units of counts per minute, cpm). The measurement volume was a disk with a diameter ∼400 m and depth ranging from 0.12 to 0.76 m (Köhli, 2015;Zreda et al., 2008). Neutron moderation power of the soil is controlled by hydrogen, so the flux of epithermal or fast neutrons detected at the soil surface is inversely proportional to soil water content (Desilets et al., 2010;Zreda et al., 2012). From the NC, volumetric soil water content (SWC) was estimated in cm 3 cm -3 with a nonlinear calibration function following Franz et al. (2015).
Gamma-ray spectra were collected with a 2.5 L NaI(Tl) scintillation crystal with 512 channels, made by Hydroinnova. The detector was mounted on the vehicle and collection time for each spectrum was 10 s. Detector position was recorded via GPS at the beginning of each 10-s measurement period. The midpoint between the detector location at the beginning of the measurement period and end of the measurement period was used as the location for the corresponding gamma-ray spectrum. Gamma-rays that are detected by the spectrometer are emitted from the top 30-60 cm of the soil. The stationary 65% footprint of the gamma-ray spectrometer, when mounted at 1.5-m height, is described by a circle that has a radius of 3.8 m. Sixty-five percent of the radiation detected by the gamma-ray spectrometer is emitted by a volume that lies within this circle. The 95% footprint has a radius of 24 m, and therefore it can be advised for the interpretation that the static spectrometer collects gamma-ray from an area that has a radius <24 m (van der Veeke et al., 2021). Using a generic calibration based on detector specifications, Gamman software (Medusa Radiometrics) was used to analyze gamma-ray spectra and determine activity concentrations of K-40, U-238, and Vadose Zone Journal F I G U R E 2 Drone-surveyed red, green, blue images, soil types, and sample locations at Sites 1, 2, and 3 Th-232. Gamman performs energy stabilization and then uses the non-negative least squares full-spectrum analysis (NNLS-FSA) approach to find radioelement concentrations (Caciolli et al., 2012;Hendriks, 2001). A digital surface model (DSM) from each field, after harvest, was created from images collected with a DJI Phantom 4 RTK (DJI) unmanned aircraft system (UAS). The aircraft is equipped with a 20 megapixel RGB camera (5,472 × 3,648 pixels), and it was flown at 61 m (200 ft) aboveground level, with front and side overlap of 75%. To assure high spatial accuracy, the UAS was connected during the flights to an internet based virtual base network (VBN) provided by DigiFarm (Monticello, IA), which resulted in images geotagged with real time kinematics precision (0.02 m accuracy). For redundancy, we used eight ground control points (GCPs) spread across each field. Lids of 22.7 L (5 gallon) pails (area of 0.07 m 2 ) were used for that purpose, and a Trimble Geo7x GPS unit (Trimble), connected to the same VBN mentioned above, was used to survey the center of each GCP (0.02-m accuracy). The images were stored in a SD card during flights, and later they were transferred to a desktop computer to be processed. The images were first processed (stitched) with Pix4Dmapper from Pix4D (Pix4D SA), resulting in a DSM with average ground sample distance across sites of 0.018 cm per pixel. Because the field level analysis did not require such high resolution as of the DSMs generated from the stitching process, ArcGIS Pro (ESRI) software was used to resample those to a 1 m per pixel resolution prior to further analysis.
All covariate data measurements were translated to a 10 by 10 m grid, where the value of each grid node was the average of all surrounding data points within a specified radius. A 70-m radius was used for all data types except for gamma-ray data at Site 1, where the radius was decreased to 31 m to avoid an artificial spatial pattern that arose when a search radius of 70 m was used. Geophysical data smoothed to the 10 by 10 m grid was then interpolated in Surfer mapping software (Golden Software LLC) with ordinary kriging to full field extent and grid cell size of 2 m to create a complete covariate table for model prediction. The spatial continuity and stationarity assumptions of kriging are believed to be reasonable in this field soil mapping scenario. Multiple variogram models (linear, spherical, Gaussian, exponential) were constructed for each geophysical data type, and the resulting interpolation with lowest median absolute deviation of residuals (from 100 randomly selected points) was chosen for model training and prediction.

Soil sampling and laboratory analyses
Sites 1, 2, and 3 were sampled on 20, 21, and 26 Oct. 2020, respectively. Typical soil sampling is done every hectare, but cost and labor required for hydraulic property analysis limited this study to 15 samples per 53 ha. Moreover, in most agricultural applications, the number of soil samples will be limited to one sample per hectare or even fewer following university extension guidelines. Optimal placement of such limited soil samples has been discussed elsewhere; see Lesch et al. (2000) for USDA soil salinity sampling based on EMI and see Gibson and Franz (2018) for soil hydraulic property sampling based on EMI and CRNS. Here, sampling locations were primarily selected based on uniform spacing. Slight position modifications were made to capture soil types based on visual examination of ECaD data, SSURGO soil zones, and elevation data. Samples were also a minimum distance of 50 m from one another. Locations of the 15 soil samples collected at each site are given in Figure  Eight soil properties were estimated in the laboratory: CEC, pH, EC, OM, BD, texture (percentage sand, silt, and clay), field capacity (FC), and permanent wilting point (PWP). Available water capacity (AWC) was determined as FC-PWP, for a total of nine soil properties. Only relatively static soil properties were considered in this analysis. Cation exchange capacity, pH, EC, and OM were determined by the Soil Testing Lab (STL) at North Dakota State University, Fargo. Organic matter content was measured by weight loss on ignition (Combs & Nathan, 2015). Bulk density was determined in the laboratory based on oven drying of field samples. Tex-ture was determined with hydrometer mechanical analyses. Field capacity and PWP on a gravimetric bases were estimated with 1/3 bar and 15 bar water contents, respectively, from pressure plate analyses of 100 g of sample. Volumetric water content values of FC and PWP were determined as the product of ρ b and gravimetric water content assuming a density of 1.0 g cm -3 for water, and adjustment for stones was made following Gardner (1986). In 14 cases where soil sampling depth was limited by gravel (particularly at depths below 0.30 m at Site 3), missing data were populated with values from NRCS Soil Survey Geographic database (Soil Survey Staff, 2019). The SSURGO "CEC-7" values were considered equivalent to the NDSU STL's CEC values (L. Cihacek, personal communication, 2021).

Soil property statistical models
Simple statistical models of the soil properties measured in the lab were built using all geophysical data types as possible predictor variables. Each sampling depth interval was modeled separately. The training set was constructed by extracting all geophysical data at the grid node closest to each of the sample locations and joining it to the sampled soil property data.
Modeling was limited to multiple linear regression using ordinary least squares because only 15 soil samples were collected for calibration at each site. All modeling and prediction was carried out using the caret package in R (Version 4.0.2).
The following components of multiple linear regression were addressed: normality of error, multicollinearity, and homoscedasticity. A log base 10 transformation was applied to all predictive data to improve normality. Because we are interested in the significance of individual predictors, multicollinearity was handled by calculating the variance inflation factor (VIF) of the model and iteratively removing the variable with the highest VIF until all VIF scores of the remaining variables were <5. The VIF is given by where 2 is the coefficient of determination of the ith predictor variable regressed against all other variables. Model residuals were plotted against fitted values to evaluate homoscedasticity (see supplemental R code). To avoid overfitting, only models with three parameters (two predictors and an intercept) or fewer were evaluated. Multiple linear regression was performed on all possible two-and three-parameter combinations using the entire training set. The set of regression models for each soil property made up of all two-parameter models and the three-parameter models with the 10 highest coefficient of determination (R 2 ) statistics was further evaluated with leave-one-out crossvalidation (LOOCV) using the "caret" package in R. P values of the final model parameters given by LOOCV describe the significance of each predictor.

Map predictions of soil properties
For all predicted soil properties and depth intervals, the model with lowest root mean square error of prediction (RMSEP; Equation 2) was chosen and model predictions were calculated using the covariate table of interpolated geophysical data. Predicted values were truncated according to physical constraints, which were set as the minimum lower and maximum upper expected value within the field area according to the SSURGO data base (Table 1). If observed values were more extreme than those expected by the SSURGO data base, the minimum and maximum expected values were substituted as constraints. Predictions were then summarized with the fol-lowing measures: RMSEP, R-squared of prediction ( 2 ), minimum, maximum, standard deviation, and mean. RMSEP and 2 are given by wherêis the ith predicted value from the final model built by LOOCV, y i is the ith observed value, N is the sample size and̄is the sample mean. The RMSEP given by LOOCV was considered a reasonable estimate of the overall uncertainty in the models and was compared to generic uncertainty levels considered useful in agricultural management. Proposed uncertainty thresholds are defined in Table 1 based on the authors' expert knowledge.

Geophysical data
Maps of kriged geophysical data are given in Figure 3. Mean and uncertainty of ECa, radionuclide concentrations, elevation, and NC are given in Table 2. Among all the sites, Site 3 has the smallest range in ECa and elevation. Elevation and ECaD were moderately correlated at Sites 1 and 2, and elevation and all EMI variables were strongly correlated at Site 3. The ECaS and ECaD were negatively correlated with cosmicray NC with correlation coefficients of r = -.7 and r = -.53 at Site 1, and correlation coefficients of r = -.4 and r = -.39 at Site 2. At Site 2, elevation and cosmic-ray NC were correlated with a coefficient of r = .55. K-40 was correlated with U-238 (r = .33) and Th-232 (r = .40) at Site 1, and negligibly correlated with the other radioelements at Sites 2 and 3 (|r | < .26). Th-232 and U-238 were positively correlated at Site 1 (r = .41), negatively correlated at Site 2 (r = -.41), and negligibly correlated at Site 3.

Soil sampling
Descriptive statistics of sampled soil properties for the 0-to-0.30-m depth interval are given in Table 3 ( sample length less than the desired minimum of 15 cm, or in one case a length recording uncertainty. Correlations between sampled soil properties and geophysical layers are shown in Figure 4. Linear correlations between soil properties and geophysical data varied among sites and soil depths. The EMI data (ECaS, ECaD, and ECaSDR) was most often at least moderately correlated (|r | > .5) with soil properties compared with the other data types. Across all three sites, ECaS was most consistently correlated with percentage sand, silt, and clay for the 0-to-0.30-m sampling interval, with correlation coefficients ranging from r = .47 to r = .73. A strong correlation coefficient of r = .93 existed at Site 1 between ECaSDR and AWC. Beyond the EMI data, correlations with other sensor data were more sporadic, such as the strong correlation between elevation and OM at Site 2 (r = -.76) and moderate correlation between elevation and OM at Site 3 (r = -.55). Elevation was also moderately correlated with BD, available water capacity, sand, and silt (r = .64, -.59, .63, -.68, respectively). K-40 was moderately correlated with EC and CEC in T A B L E 2 Mean and uncertainties are reported for measured shallow apparent bulk electrical conductivity (ECaS), deep apparent bulk electrical conductivity (ECaD), elevation (Elev.), neutron counts (NC), K-40, U-238, and    the 0-to-0.30-m depth interval at Site 1 (r = .42, .56, respectively). At Site 1 in the 0.61-to-0.91-m depth interval, U-238 and Th-232 were correlated with AWC, sand, silt, and clay with absolute value of correlation coefficients between r = .45 and r = .58. U-238 and Th-232 also had a moderate to strong correlation in the 0.61-0.91-m interval at Site 1 with EC, OM, and BD, with absolute value of correlation coefficients between .53 and .83. U-238 had a moderate to strong correlation with texture, CEC, and AWC at Site 2. The only noteworthy radionuclide correlations at Site 3 were r = .58 and r = .59 for U-238 with CEC and pH, respectively. The SWC was moderately correlated with EC, CEC, and OM at Site 1, and was also moderately correlated with CEC at Site 3.

Multiple linear regression results
The number of soil properties and depth interval pairs (27 possible) that could be modeled by multiple linear regression with 2 close to .5 (>.4) was 13, 14, and 6 at Sites 1, 2, and 3, respectively (Table 4). The depth interval most often modeled successfully was 0-0.30 m. ECaS, U-238, NC, and elevation were eliminated as possible predictor variables at Site 1 to reduce multicollinearity. At Site 2, ECaS, Th-232, and NC were eliminated as possible predictor variables. At Site 3, ECaS, ECaD, ThUR, and NC were excluded from prediction. Plots of residuals vs. fitted values showed minor heteroscedasticity in some CEC, EC, and pH models, but no correction was attempted. Only models with the lowest RMSEP for each soil property and depth pair are reported in Table 4. An exhaustive model summary is available in supplemental R code document. At Site 1, pH models for 0.30-0.61 and 0.61-0.91 m with K-40, Th-232, or ThUR as predictors had 2 > .5 but p values >.65. Significant (p value < .01) predictor variables for the EC models in the 0-0.30-m interval were ECaSDR, ThUR, and SWC. A significant (p value < .01) model for OM in the 0.61-to-0.91-m interval with 2 was predicted with Th-232 and ThUR. All possible predictive variable combinations involving ECaSDR at Site 1 were able to predict AWC from 0 to 0.30 m with 2 >.8. The p values of the secondary predictor variable coefficients for AWC models range from .26 to .95 while the p values of the ECaSDR coefficients were between 3.61 × 10 -6 and 5.48 × 10 -7 . Models for 0-0.30 m AWC at Site 1 that did not contain ECaSDR had large overall p values (p value > .15). The ECaD was a significant predictor for AWC in the 0.30-to-0.61 and 0.61-to-0.91-m depth intervals. Models of 0-0.30 m sand, silt, and clay at Site 1 predicted with ECaSDR in any combination all achieved an 2 >.65 and overall p values < .01.
At Site 2, significant models with 2 > .5 predicted BD with ECaD as the primary predictor and K-40 or ThUR as the secondary predictors. Elevation with SWC and elevation with K-40 predicted OM with 2 > .5 in the 0-to-0.30-m interval. For the OM in the 0.61-to-0.91-m interval, all models including ECaD were significant and achieved 2 > .5. All models at Site 2 that included ECaD as a predictor, excluding the ECaD and elevation model, were significant and achieved 2 >.48 for 0-0.30 m AWC. The 0.30-0.61 m AWC was also well-predicted with all models containing ECaD; the lowest 2 was .53. Behind ECaD, the next best predictors for the first and second depth intervals of AWC were U-238 ( 2 = .42) and K-40 ( 2 = .37), respectively. Sand models for 0-0.30 m with 2 > .5 were predicted with ECaD and U-238, elevation, or ECaSDR. All clay models for 0.30-0.61 and 0.61-0.91 m that included ECaD achieved 2 T A B L E 4 Multiple linear regression models with lowest root mean square error of prediction (RMSEP), where models with 2 > 0.4 are underlined Note. Min., minimum; Max., maximum; RMSEP, root mean square error of prediction; 2 , R-squared of prediction; EC, electrical conductivity (mmhos cm -1 ); CEC, cation exchange capacity (cmol c kg -1 ); BD, bulk density (g cm -3 ); OM, organic matter; AWC, available water capacity (cm cm -1 ); ECaS, shallow apparent bulk electrical conductivity; ECaD, deep apparent bulk electrical conductivity; ECaSDR, ratio of shallow to deep apparent bulk electrical conductivity; K40, potassium; U238, uranium; Th232, thorium; ThUR, ratio of thorium to uranium; SWC, soil water content from cosmic-ray neutron probe. Model descriptions assume that an intercept is also included.

F I G U R E 5
The underlined models in Table 4 are mapped in space at Site 1. Soil properties predicted are pH, electrical conductivity (EC; mmhos cm -1 ), cation exchange capacity (CEC; cmol c kg -1 ), bulk density (BD; g cm -3 ), percentage organic matter (OM), available water capacity (AWC; cm cm -1 ), percentage sand, percentage silt, and percentage clay values between .45 and .66. K-40 alone predicted clay with 2 of .47 and .61 for the 0.30-to-0.61-and 0.61-to-0.91-m intervals, respectively.
One soil property was predicted at Site 3 with 2 > .50. EC was predicted by elevation in the 0.30-to-0.61-m interval with 2 = 0.84, but the model was not significant (p value = .94). A model for AWC in the 0.0-to-0.30-m depth interval with Th-232 and ECaSDR as predictors achieved an 2 of .46. Organic matter was predicted with K-40 and U-238 from 0.61 to 0.91 m with an 2 of .45. For all the soil properties combined, U-238 and elevation were the predictors with p values most frequently <.01. Overall p values were <.01 for 24 models at Site 3 compared with 82 and 147 models at Sites 1 and 2, respectively.

Spatial predictions of soil properties
Underlined models in Table 4 were mapped spatially in Figures 5-7. When comparing predictions to uncertainty thresholds (Table 1), 11 models pass at Site 1, 16 at Site 2, and 14 at Site 3. Despite meeting uncertainty thresholds, most of the models at Site 3 still have a low 2 , showing that the variation within the field may be too low for a spatial prediction to add any useful information to the mean. There were also model cases where 2 was >.5 but the uncertainty was too high (Table 4). Soil properties that were predicted at Site 1 with both an 2 > .5 and low enough uncertainty were EC, AWC, silt, and clay. At Site 2, BD, OM, AWC, silt, and clay were predicted with 2 > .5 and acceptable uncertainty. The AWC was the only soil property predicted at Site 3 with acceptable uncertainty and 2 close to .5 ( 2 = .46).

DISCUSSION
Because a number of the statistical models produced favorable 2 values, coupling of geophysical surveys with a small number of soil core samples with MLR proved to be a reasonable strategy that invites further development. The small number of soil samples (15) collected at each site in this study was on par with the reality that measuring soil core properties will always be time-consuming and expensive. Incorporating Vadose Zone Journal F I G U R E 6 The underlined models in Table 4 are mapped in space at Site 2 F I G U R E 7 The underlined models in Table 4 are mapped in space at Site 3 geophysical data allowed us to infer how the measured soil properties vary across a field with much smaller collection and analysis time than would be required for extensive grid soil sampling.
Our results agree with the inconsistent predictor-response relationships presented by Rodrigues et al. (2015) and Wong and Harper (1999). The models trained at each of the three sites were also applied to the other two sites to explore potential for a universal calibration among our sites. However, 74% of the resulting predictions had more variance between the predicted soil properties and the sample mean than between the observed soil properties and the sample mean. The results indicate that local calibration is still needed, and that further work needs to be done to approach a universal calibration for these sites similar to Rodrigues et al. (2015). At the sites in this study, EMI (ECaS, ECaD, ECaSDR) appeared to be the most useful geophysical layer. While GRS ThUR) and CRNS (NC, SWC) data added some information, the correlations with soil properties were not uniform across sites or soil depth intervals, and the model coefficients for GRS and CRNS variables were usually less significant than those of EMI when used as joint predictors.
These differences in the sensors' predictive performance raise a short discussion of the sources of error and bias introduced by their differences in measurement frequency and sample volume. Because all sensors were traveling at the same speed and the EMI sensor had the highest sampling frequency at one measurement per second, the distance between consecutive EMI measurements was smaller (2.78 m) than the distance between consecutive GRS and CRNS measurements (27.78 and 166.78 m, respectively). The fact that more EMI data points were collected than GRS and CRNS data points and that EMI soundings are most similar in scale to soil samples may explain why EMI generally performed as the best predictor in our models. However, despite the better fit in scale between EMI data and soil samples, Gibson and Franz (2011) found that CRNS was more strongly correlated with soil hydraulic properties than EMI. The sample volumes for the GRS and CRNS (circles with radii roughly 24 and 200 m, respectively) are large enough to capture the information between the larger measurement spacings so that the same ground is covered even though fewer measurements are taken. At the same time, both the GRS and CRNS are more sensitive to the ground volume closer to the detector. In addition to differences in horizontal spatial scale, the authors acknowledge errors introduced by differences between vertical scales of the soil sampling depth intervals and the sensing depths of the geophysical sensors. Each sensor has a slightly different sensitivity function with depth below the soil surface, and inversion of geophysical data was not attempted for the sake of simplicity in future practical use. We refer the reader to other papers that have studied inversion of EMI data (Callegary et al., 2007;Piikki et al., 2014) and the vertical and horizontal support volume of the CRNS and GRNS more in-depth (Köhli et al., 2015;van der Veeke et al., 2021).
For soil cores in this study, uniform sampling combined with ECaD, SSURGO, and elevation information was used to intuitively select informative locations for a few soil samples. The authors acknowledge that consideration of ECaD data in selection of soil sample locations (and not GRS and CRNS data) introduces bias toward EMI data as a successful predictor in soil property models. However, our sample selection method was constrained by performing the geophysical data collection and the soil sampling very close to one another temporally because of the need to collect the soil samples before an impending snowstorm. Thus, the selection of sample locations was done in rapid time, primarily directed by the need to uniformly space our allotted 15 samples. We did not have enough time to process the gamma-ray and neutron data to include in the sample selection process. Future work could examine the relative success of EMI, GRS, and CRNS data when all three sensor data types are formally considered in soil sample selection. Despite introduction of bias toward certain predictive data types, development of intelligent sampling strategies based on elevation, SSURGO zones, and geophysical surveys is essential for maximizing the potential of MLR for situations with sampling constraints. One situation with limited soil sample size is the SOC storage industry. For verification of SOC storage, practical recommendation by the Food and Agriculture Organization of the United Nations (FAO) is currently one composite sample per 10 ha (FAO, 2020). This would be equivalent to collecting only five soil samples per site in this study. The first basic approach the FAO recommends is stratified simple random sampling, where at least three strata -or zones -are determined by dividing the area of interest equally. The second approach is direct stratified sampling, where at least three strata are identified from previous information such as ECa maps. In both methods at least three composite soil samples are randomly collected in each stratum. Clearly, selecting representative sample locations of SOC is highly important for accurately understanding C storage and enabling the C credit industry to become lucrative for producers. A promising approach to intelligent selection of sample locations is to use k-means clustering of geophysical and other available data. For example, van Arkel and Keleita (2014) used k-means clustering of EMI and topography data to select critical soil moisture sampling locations for estimation of mean field-scale soil moisture and concluded that the approach was a good alternative to other selection methods because it performed nearly as well or better without the need for extensive pre-sampling. In a hypothetical method, k-means clustering of ECa and elevation data in feature space would create map zones within which centroidal voronoi tessellations could be calculated with Lloyd's algorithm (Lloyd, 1982). Sample locations would then be the centers of the voronoi tessellations within the different zones characterized by ECa and elevation. Extra constraints are required to ensure that k-means clusters in feature space are mapped to zones that are concave, connected, and large enough in real space.
Although the MLR predictions were moderately successful, we also acknowledge limitations in our predicting capabilities. At Site 3, we were surprised to find that none of the significant MLR models met 2 > .5. Probable cause of failure was low variability within Site 3 to the best of our knowledge. In addition, soils at Site 3 are highly disturbed because the site has been previously used for gravel production due to its proximity to the James River. The variance in ECaD at Site 3 was 10.55 mS m -1 compared with variances of 58.32 and 47.55 mS m -1 at Sites 1 and 2, respectively. Furthermore, sampled soil properties at Site 3 also tended to have a lower standard deviation than the measured soil properties at Sites 1 and 2 (Table 3). Due to the low variability, the best statistical inferences for soil properties at Site 3 were the sample means instead of regression model predictions. In addition to within-field variability, transient soil temperature and soil moisture may have also influenced EMI performance at the different sites. Brevik et al. (2006) found that the difference in ECa readings between different soils decreased as soil moisture decreased. Gibson and Franz (2018) found that the use of multiple mapping times combined with Empirical Orthogonal Functions could reduce the transient effects of soil temperature and soil moisture on EMI and CRNS data. The authors are not confident that type of parent material is a cause of statistical modeling failure. Site 3 material does appear to be less weathered because gravelly sands are dominant below about 0.5 m, but the parent materials at Sites 2 and 3 are both river-deposited cross-bedded sand and plane-bedded gravel. However, collection of soil mineralogy data in the future could distinguish between quartz sands and K-feldspar sands to allow better interpretation of the gamma-ray data (Heggemann et al., 2017;Priori et al., 2014). Inclusion of soil mineralogy data in the models would likely improve the predictive ability of gamma-ray data at all three sites. Small sample size and nonlinear relationships between predictor and response variables may have also limited MLR model performance, which typically improves with increased sample size (Khaledian & Miller, 2020).
In the future, the limitation of not knowing whether a site is a good candidate for linear regression may be overcome by examining readily available online data such as elevation and SSURGO data. For instance, Lo et al. (2016) used SSURGO data to calculate the field-averaged amount of undepleted available soil water in the root zone under conventional irrigation by examining differences in rootzone water holding capacity among soil units. Soil water at the end of the growing season above a critical management threshold was considered undepleted, and Lo et al. used the field average amount of undepleted available soil water to estimate benefits of implementing variable rate irrigation for mining undepleted soil water through planned depletion. A similar approach based on variability in SURRGO data may be possible to estimate feasibility of multiple linear regression modeling at a given site. Another practical development required for implementation in precision agriculture is to compile evidence-based uncertainty targets for MLR predictions. Work needs to be done to determine the actual uncertainty level required in each soil property to make a variable management decision. It would also be beneficial to determine decision threshold values for certain treatments where there is a yes-no decision. These thresholds would be used to evaluate a statistical model's ability to differentiate a field into zones that are either below or above a given decision threshold. Finally, both GRS and CRNS require investment in data processing, whether through obtaining expert support in software use, purchasing sensors with embedded software, or spending extended time learning the details of the method. This time or monetary cost motivates a recommendation of whether GRS and CRNS data is worthwhile at future sites. Although further study is needed, GRS and CRNS data may not be essential at sites with similar pedology and variation to those in this study.

CONCLUSIONS
This paper contributes to understanding the usefulness of geophysical data types by introducing results from three new pedological and geographic settings in North Dakota. Over half of the best soil property predictions were based on multiple data sources instead of data from a single sensor (Table 4). Statistical models at two of the three sites met expert opinion for uncertainty targets in variable management decision-making. It is unclear which geophysical data type is expected to be the best predictor a priori at a given location. Our predictions were site-specific, and models trained at one site performed very poorly at the other sites. Our understanding of sensor performance and its relationship to field conditions and sensor support volumes could be further refined by incorporation of more information such as soil mineralogy and spatial variability of soil property values. However, based on results from these sites with our simplistic method, the authors recommend prioritizing EMI surveys if geophysical data collection is limited to a single mapping effort and calibration soil samples are few.

A C K N O W L E D G M E N T S
We acknowledge the support from USDA National Institute of Food and Agriculture Foundational Program Cyber-physical systems (2019-67021-29312). T.E.F. and D.D.S acknowledge the financial support of the USDA National Institute of Food and Agriculture, Hatch project nos. 1009760 and ND0140, respectively. Financial support was provided by the Joint FAO/IAEA Programme of Nuclear Techniques in Food and Agriculture through the Coordinated Research Project (CRP) D1.20.14 Enhancing agricultural resilience and water security using Cosmic-Ray Neutron Sensor (2019-2024). The research presented is also a result of a joint study by USDA-NRCS North Dakota (under agreement no. NR206633XXXXC001), North Dakota State University, and the University of Nebraska-Lincoln on irrigation water management and variable rate irrigation technologies; and was partially supported by the North Dakota Agricultural Experiment Station. We appreciate the equipment, time, and expertise for soil sampling provided by Jordaan Thompson-Larson, Erica Althoff, and others of the USDA-NRCS staff of North Dakota. We thank Sheldon Tuscherer, Dongqing Lin, and Mathew Blum for field and laboratory support. Mention of equipment and trade names is for information only and is not intended to constitute endorsement by the authors, their respective organizations, or the research sponsors, of one product to the exclusion of others available for similar purposes.