Introduction

An optimal range in soil pH is one of the fundamental prerequisites for successful agricultural crop production. However, under humid climates of Central Europe, soils tend to acidify as basic cations, in particular Ca2+ and Mg2+, are continuously replaced by H+ and Al3+ ion from the exchange sites of soil colloids (European Commission, 2012). This natural soil genetic process is accelerated by human activity, including the deposition of nitric and sulphuric acids through precipitation, the H+ released through bacterial nitrification of ammonium-based fertilizers, animal manures and residue/green manure decomposition and the removal of base nutrients by harvesting the cultivated crops (Brady & Weil, 2008; Kerschberger & Marks, 2007). Since soil pH affects many key determinants for plant growth, acidification needs to be regularly compensated by lime application in soils without carbonate buffering capacity in order to prevent yield loss and replenish soil fertility (Goulding, 2016).

To determine the lime requirement (LR) of an agricultural soil, its soil acidity and the quantity of bases that must be applied to replace reserve acidity on clay and organic matter, resulting in raising soil pH from acidic to neutral, needs to be determined. As soils are buffered systems that resist changes in pH, the measurement of pH alone, active acidity, is not sufficient to determine soil LR (Godsey et al., 2007). In contrast, the quantity of bases needed to replace the reserve acidity and increase the pH value to an optimum value to best support a particular crop rotation can be directly quantified by determining the base neutralizing capacity (BNC) of the soil (Blume et al., 2016; Vogel et al., 2020). The BNC is measured in the laboratory by discontinuous titration adding increasing concentrations of a basic solution (e.g. Ca(OH)2) to a soil sample. After a defined equilibration time, the pH increase of the solution is measured potentiometrically to determine the amount of acid neutralized. From the resulting titration curve, the quantity of lime required to obtain a desired target pH value can be determined.

When a site-specific liming strategy is to be developed, LR data at a very fine spatial resolution are needed to account for within-field spatial soil variability. However, the cost of laboratory analysis for a sufficiently dense soil sampling grid necessary to define site-specific lime rates is prohibitive, hence either insufficient samples are collected to identify lime rates to specific areas, or a whole-field approach is used, resulting in inadequate lime recommendations to optimize field pH and might result in a waste of resources. Consequently, alternative strategies for LR determination for precision agriculture application are necessary.

A few studies have demonstrated the potential of using proximal soil sensors to derive lime prescription maps. Viscarra Rossel and McBratney (1997) and Viscarra Rossel et al. (2005) reported the results of using a prototype on-the-go soil pH and lime requirement measurement system that consists of a soil sampling and sieving mechanism, a soil analytical component using a pH ISFET (ion sensitive field effect transistor) in order to measure soil pH changes in combination with data collection and measurement algorithms. Based on the work of Adamchuk et al. (1999), Lund et al. (2005) developed an automated soil sampling system for on-the-go measurement of soil pH integrated on the Veris Mobile Sensor Platform (MSP; Veris technologies; Salinas, KS, USA) and combined with an apparent electrical conductivity (ECa) sensor and near-infrared spectroscopy (NIRS) to generate soil pH and lime requirement maps. Kuang et al. (2014) utilized an on-the-go visible and near-infrared (visNIR) spectroscopy sensor to map the within-field variation of organic carbon, pH and clay content for LR determination. Von Cossel et al. (2019) deployed low-input sensor-based soil mapping of electromagnetic induction (EMI) (EM38 MK I; Geonics, Canada) in combination with in situ and ex situ pH measurements to determine the soil’s LR. Bönecke et al. (2020) applied an on-the-go multi-sensor approach using the Veris Mobile Sensor Platform and the Geophilus Electricus (Geophilus GmbH, Germany) proximal soil sensing system (Lück and Rühlmann, 2013). By combining the sensor data with reference analyses of soil characteristics that are well-correlated with soil acidity and soil pH buffer capacity, high-resolution soil maps of pH, texture and soil organic matter (SOM) were generated. These were utilized to develop site-specific LR maps based on the standard liming algorithm of the Association of German Agricultural Investigation and Research Institutions (VDLUFA; von Wulfen et al., 2008).

Most of these methods rely on specifically calibrated transfer functions of the sensor data with laboratory measurements of soil properties that affect soil acidity, pH buffer capacity to determine soil LR. The main drawback of this approach is their indirect estimation of the LR via soil acidity affecting soil properties, necessitating their cost and time-consuming determination through laboratory analyses. To overcome this problem, a more straightforward method might be preferable which directly determines the LR through use of multiple sensor data. Determining the soil BNC is such a direct LR method that determines the effect of lime addition on the pH, individually for each soil sample in order to determine LR from a resulting titration curve (Vogel et al., 2020).

The objective of this study was to develop a site-specific variable rate LR procedure that is based on the soil’s BNC in combination with multi-sensor platform soil mapping. Specific objectives were: (i) to assess the quality of a multi-sensor-based prediction of LRBNC, (ii) to determine the sensor or sensor combination(s) most sensitive to describe the LRBNC, (iii) to quantify the within-field spatial variability of LRBNC and (iv) to generate high-resolution LRBNC maps of the test fields.

Materials and methods

Site description

Nine agricultural fields were selected on three farms in a quaternary landscape of northeast Germany. They show field sizes between 20 and 76 ha (Table 1). The study area is part of the Northeast German Plain which belongs to the broader geomorphological region of the North European Plain (Fig. 1). It was largely formed by the Pleistocene glaciations of the terrestrial Scandinavian ice sheets as well as by subsequent periglacial and interglacial Holocene geomorphic processes. In the study area, the present-day landforms and soils were particularly shaped by the advances of the Weichselian (115–12 ka) and the preceding Saalian glacial belt (150–130 ka; Krbetschek et al. 2008). Climatically, the study area is situated in a transitional zone between oceanic climate of Western and continental climate of Eastern Europe. Due to a relatively low altitudinal range of the land surface of ~ 0 to 200 m above sea level, regional climatic differences are small. Thus, following the Koeppen–Geiger Climate Classification System, the climate of the study region is classified as temperate oceanic with an increasing influence of continental circulations. The mean annual air temperature is ~ 9 °C. The coldest and warmest months in a year are January and July with mean temperatures of -1 and 18 °C, respectively. With a mean annual total precipitation of less than 550 mm, it is one of the driest regions in Germany.

Table 1 Descriptive statistics of standard laboratory results of reference soil samples
Full size table
Fig. 1
figure 1

Locations of BNC sampling sites on the studied fields (A: GW32; B: GW6; C: KL41 [left], 42 [right]; D: KL6; E: KL60; F: GW21; G: PP1401 [left], 1392 [right]). Note the two different scales, Projection: UTM ETRS89 33 N; Aerial Photographs: Google | DigitalGlobe; Bing Maps, Microsoft

Full size image

The three farms are located in the east (Komturei Lietzen, KL, Lat: 52.483766, Long: 14.333079); Landwirtschaft Petra Philipp, PP, Lat: 52.376035, Long: 14.461919) and in the north (Gut Wilmersdorf, GW, Lat: 53.110092, Long: 13.909461) of the federal state of Brandenburg (Northeast Germany). They are mainly located in the Pleistocene young morainic landscape of the Weichselian glaciation as well as in the Holocene river valley of the Oderbruch showing high within-field soil variability. In accordance with the German soil classification system KA5 (Eckelmann et al., 2005), soil textures range from pure sand (class: Ss) to loamy clay (class: Tl) showing a dominance of sand and loam (classes: Sl, Su, St, Ls). Even though the development of soil acidity in most of these soils require regular lime amendment, in some soils, the pH is greater than 7 due to the presence of surface soil carbonates embedded as part of the glacial till and landscape position. The crop rotation of all fields is cereal-dominated.

Proximal sensor mapping

In 2017 and 2018, nine fields were mapped using the Mobile Sensor Platform (MSP) developed and manufactured by Veris Technologies™ (Salinas, KS, USA). It is currently the only multi-sensor system commercially available for obtaining simultaneous potentiometric, geoelectric and optical measurements (Fig. 2, Bönecke et al., 2020):

  1. i.

    Potentiometry:

    On-the-go potentiometric measurements (Fig. 2b) are performed by two ion selective antimony electrodes on naturally moist soil samples. While driving across the field, a sampler shank is lowered into the soil and a soil core flows through the sampler trough. Next, the soil sampler is raised out of the soil and presses the sample against the two electrodes for two separate measurements. Then, the arithmetic mean of the two voltage measurements is recorded. After measurement, the sampler shank lowers back into the soil to replace the old soil core at the back end of the sampler trough when the new sample enters from its front end. Meanwhile, the electrodes are cleaned with water by two spray nozzles and the device is ready for the next measurement process. During field operation, the measurements are georeferenced by differential global navigation satellite system (GNSS) co-ordinates that are recorded when the sampler shank is raised out of the soil. The conversion of voltage into a pH value is realized by a preceding calibration with pH 4 and pH 7 standard solutions (Schirrmann et al., 2011). Depending on the measuring time required for the sensor, pH values are recorded every 10 to 12 s (Lund et al., 2005; Fig. 3a).

  2. ii.

    Optical reflectance:

    Soil reflectance has been studied since the 1970s as an effective means for estimating the SOM content of the soil (Sudduth & Hummel, 1993). In the present study, the Veris OpticMapper was used (Fig. 2c). It is a dual-wavelength on-the-go optical sensor measuring differences in the diffuse light reflectance of the soil. It consists of a single photodiode and two light sources (red LED, wavelength: 660 nm; NIR LED, wavelength: 940 nm). In its forward face to the direction of movement, the OpticMapper positions a coulter to cut crop residues while the optical module is mounted on the bottom of a furrow ‘shoe’ between two side wheels that set the sensing depth. The wear plate is pressed against the bottom of the furrow approximately 40 mm below the soil surface with a consistent pressure to provide a self-cleaning function. The modulated light is passed through a sapphire window onto the soil. The photodiode then receives the modulated reflected light and converts it into a voltage, which is further processed and logged. The optical data and GNSS co-ordinates are recorded at a rate of 1 Hz (Kweon & Maxton, 2013; Fig. 3C, D).

  3. iii.

    Geoelectrics:

    The apparent electrical resistivity (ERa) was measured at a rate of 1 Hz with a galvanic coupled resistivity instrument using six parallel rolling coulter electrodes (Fig. 2a, 3b). ERa values are internally converted into apparent electric conductivity (ECa) output. The electrode configuration provides readings over two depths with median depths of exploration of 0.12 (ECa shallow) and 0.37 m (ECa deep; Gebbers et al., 2009). This enables the identification of significant soil textural and/or soil moisture changes between soil horizons. Since pH and OpticMapper measurements are carried out in the topsoil, only ECa shallow readings were used in the present study.

Fig. 2
figure 2

A Mobile Sensor Platform (MSP) developed and manufactured by Veris Technologies™ (Salinas, KS, USA) with ECa instrument (1), B soil pH Manager (water tank (2), soil sampler (3) with sample (4), pH electrodes (5), and cleaning nozzles (6), and C OpticMapper with and optical module (7) in between the ECa coulter electrodes (photos: Torsten Schubert)

Full size image
Fig. 3
figure 3

Veris sensor data of test field GW6 showing (A) pH, (B) ECa [mS m−1], (C) OpticMapper Red [dimensionless], and (D) Optic Mapper Infrared [dimensionless]

Full size image

LR based on base neutralizing capacity (LRBNC)

The soil sampling sites for BNC laboratory analysis were taken in accordance with the procedure proposed by Bönecke et al. (2020), considering that the targeted samples: (i) cover the entire range of sensor data (feature space), (ii) are spatially representative, and (iii) are well distributed throughout the area of investigation. The sample size per field varies between 10 (KL41) and 35 (PP1392).

A total of 164 soil samples (Fig. 1) were analyzed for base neutralizing capacity (BNC). The BNC is defined as the amount of soil acidity that is neutralized by a base in a given time interval to a certain pH value (Meiwes et al., 1984). To directly determine the LR of the soils studied based on their base neutralizing capacity (LCBNC), the protocol of Meiwes et al (1984) was followed (Utermann et al., 2000).

In detail, the protocol included the following steps: The soil samples were air-dried and passed through a 2 mm sieve. Then, 150 g of each sample was divided into six subsamples of 25 g. One of these subsamples served as a control and was mixed with 50 ml deionized water, while the other subsamples were mixed with 25 ml of 2 N CaCl2 and 25 ml of 8 N NaOH solutions of five concentrations. This yielded six concentration levels of Ca(OH)2 added to the soil: 0, 0.25, 0.5, 1.25, 2.5 and 5 mmolc (25 g soil)−1. By adding Ca2+ and Na+ ions to the soil solution H+ and Al3+ ions are desorbed from the surface of soil colloids and neutralized by OH− ions (Meiwes et al., 1984). After 18 h of mechanical shaking, pH values were measured with a glass electrode (WTW SenTix® 81, Xylem Analytics, Weilheim, Germany) in the supernatant solution. For quantification of the buffering, the pH values and their corresponding concentrations of Ca(OH)2 added were displayed in a scatterplot and a titration curve was fitted to the six points. Based on the model, the amount of Ca(OH)2 in mmolc (25 g soil)−1 to achieve a target pH of 6.5 was derived and converted to kg CaCO3 (ha*dm)−1 by multiplying by 2,000 (Meiwes et al, 1984; Utermann et al., 2000). Because fertilization guidelines for the United Kingdom (Defra, 2010) and most other countries advise farmers to maintain a soil pH of 6.5 for cropped land (Goulding, 2016), this was chosen as the target pH value. Of course, choosing a pH of 6.5 does not reflect the fact that arable crops differ in their sensitivity to soil acidity.

Standard laboratory analyses of soils studied

To provide a general field-wise characterization of the soils studied, the following laboratory analyses where carried out on oven-dried (75 °C) and sieved (< 2 mm) soil samples:

  1. i.

    The soil pH value was measured in 10 g of soil and 25 ml of 0.01 M CaCl2 solution following DIN ISO 10390. The pH was measured with a glass electrode after a reaction time of 60 min.

  2. ii.

    The particle size distribution of the fraction < 2 mm was determined according to the German standard in soil science (DIN ISO 11277) by wet sieving and sedimentation after removal of organic matter with hydrogen peroxide (H2O2) and dispersal with 0.2 N sodium pyrophosphate (Na4P2O7).

  3. iii.

    Soil organic carbon (SOC) was analyzed by elementary analysis using the dry combustion method (DIN ISO 10694) after removing inorganic carbon with hydrochloric acid. Finally, the amount of SOC was converted into SOM following Eq. 1 (Peverill et al., 1999) assuming that SOM contains approximately 58% of organic carbon:

    $$ SOM\left[ {\text{\% }} \right] = SOC\left[ {\text{\% }} \right] \cdot 1.72 $$
    (1)

Titration curve fitting and sensor-based prediction of LRBNC

All data were processed with the free software environment for statistical computing and graphics R (version 3.6.1) (R Core Team, 2018). To fit a BNC curve to the six titration points, non-linear regression modeling was conducted using the nls function. The sensor-based prediction of LRBNC was done using a stepwise multi-variate linear regression (MLR) analysis with forward selection (R package ‘caret’; Kuhn, 2020). It iteratively adds the most contributive independent variables to the predictive model until the model improvement is no longer statistically significant. This aims at finding the combination of variables, which achieve the best model performance minimizing the prediction error (Bruce & Bruce, 2017; James et al., 2014). The MLR models are of the type:

$$ z = b_{0} + b_{1} x_{1} + b_{2} x_{2} + \ldots + b_{n} x_{n} + \varepsilon $$
(2)

where \({\text{z}}\) represents the dependent variable, x1, x2, …, xn the ancillary data measured at the same site, b0, b1, b2, …, bn the n + 1 regression coefficients, and ε the random error. In order to assess the explanatory power of each regression model, the adjusted coefficient of determination of the linear regression between predicted and measured values (adj-R2) was determined, considering the number of covariates in the model. Moreover, the average prediction error (RMSE) was calculated in a tenfold cross-validation dividing the dataset into k folds, using k – 1 folds for training and one fold for validation and repeating that procedure k times, each having a different fold for validation.

Since some of the sensor data may be correlated, prior to the modeling, the independent variables were tested for interdependencies. If two variables received a Pearson’s R greater than 0.5, one was predicted from the other using a univariate linear regression (ULR) model. Then the residuals (e) of that model were calculated following Eq. 3:

$$ e = y - \hat{y} $$
(3)

where y is the observed value and \(\hat{y}\) the predicted value. The residuals were then utilized as an uncorrelated substitute for one of the correlated independent variables. By that procedure, only the information content that is unique for each independent variable is included in the MLR analysis.

To gain an increased understanding of the relationships between the independent variables and the dependent variable as well as to identify the features of the multivariate sensor data having the greatest effect on the model performance of LRBNC, a sensitivity test was carried out by computing standardized regression coefficients (SRC). Before conducting the stepwise MLR analysis, the sensor data were standardized by subtracting the sample mean from the original values and dividing by the sample standard deviation in order to remove the influence of different units and to place all covariates on the same scale. By that procedure, the SRC of the best performing MLR model are a direct measure of sensitivity, i.e. indicative of the magnitude of influence of a single sensor datum on the LRBNC model as a whole (Hamby, 1994; Saltelli et al., 1993).

Generation of LR maps and analysis of spatial variability of LRBNC

The point-based sensor data were interpolated using the geostatistical method of ordinary block kriging (R package ‘gstat’; Pebesma, 2004) with robust variogram estimation, outlier elimination and weighted least squares approximation. For a more detailed description of the applied interpolation method, the reader is referred to Boenecke et al. (2020). For regionalization of the stepwise MLR analysis, the best performing MLR models were finally applied to the raster-based sensor data using a GIS-based raster calculator in order to generate LRBNC maps of the nine study fields.

The spatial variability of BNC-based lime requirement (LRBNC) for the nine test fields was quantified by semivariance analysis (Deutsch & Journel, 1998; Goovaerts, 1997; Webster & Oliver, 2007). The semivariogram can provide information about the maximum of semivariance (sill parameter) as well as the range of spatial autocorrelation (range parameter; Webster & Oliver, 2007). Additionally, the nugget parameter summarizes the measurement error and sample micro-variability. For semivariogram modeling, firstly, the method of moments (Webster & Oliver, 2007) was used to obtain the empirical semivariogram, which relates the average squared differences between observed values to their respective distance class (lag interval). Secondly, a theoretical semivariogram model was fitted to the empirical semivariogram using robust estimates to prevent effects from extreme outliers (Cressie, 1993). For variogram fitting, weighted least squares approximation (fit method 7 in gstat) as well as localized cut-offs were used at distances when a first local maximum was reached or the model first flattened out. From the parameters of the empirical semivariogram, insights into the spatial variability of LRBNC can be gained. The sill parameter refers to the magnitude of variability. The range parameter, on the other hand, defines the spatial context in which variability is expressed, where smaller ranges indicate small-scale distribution patterns. Taking into consideration the interdependence between sill and range, a high spatial variability is characterized by a high sill and low range value. Hence, the sill-to-range ratio (SRR) can be used in order to quantify spatial variability of LRBNC at the scale of soil sensing utilized in the operation.

Results and discussion

Field-wise soil characterization regarding acidity and lime requirement (LR)

The spatial statistics regarding the most relevant LR affecting soil properties as well as of the BNC parameters of the nine fields are shown in Tables 1, 2 and 3, respectively. The pH of the soils have median values between 5.9 and 6.6 indicating only little to no LR. However, minimum and maximum pH values of 3.8 to 5.3 and 6.7 to 7.3 show that a high within-field soil variability exists and, thus, demonstrate that the median pH value of a field alone is rather insufficient to serve as an indicator for LR determination. SOM contents are rather low throughout the study region having minima of 0.8 to 1.1%, maxima of 1.7 to 5.6% and median values of 1.2 to 2.8%. This situation is typical for the geologically young, sandy, non-stagnic soils in Brandenburg. The dominating soil texture classes (following the German Soil Texture Classification System KA5; Eckelmann et al., 2005) reveal high soil heterogeneities on the investigated fields even though sandy textures prevail showing loamy sands (Sl), silty sands (Su), clayey sands (St) and clayey sandy loams (Lts). However, also sections of pure sand (Ss) and loamy clay (Tl) exist which results in a very differentiated soil acidity and LR.

Table 2 Descriptive statistics of BNC analysis of reference soil samples I
Full size table
Table 3 Descriptive statistics of BNC analysis of reference soil samples II
Full size table

From the BNC analysis, it can be seen that the target pH increase after the addition of increasing amounts of Ca(OH)2 can be described as an exponential growth curve where the pH value reaches a threshold value when the quantity of lime tends to infinity. It has the form:

$$ pH_{target} = \alpha - \beta \cdot \gamma^{{Ca\left( {OH} \right)_{2} { }}} $$
(4)

where α, β and γ are the regression coefficients of the exponential function. For a general description of the BNC data and the soil’s pH buffer capacity (pHBC), the reader is referred to Vogel et al. (2020). The field-wise characterization of the BNC reveals that the total pH increase over all base additions (δpHtotal) strongly varies between the fields having minima of 2.5 (KL60) to 5.1 (KL41) and maxima of 6.2 (PP1392) to 7.6 (KL60) pH units. This is caused by different acidity and pHBC characteristics of the investigated soils. The range of δpHtotal per field, on the other hand, is a function of within-field variability of pHBC, showing a minimum of 0.9 pH units for KL41 and a maximum of 5.1 pH units for KL60. That means that KL41 is much more homogeneous in terms of pHBC and soil acidity than KL60, as corresponding to the pronounced variability of soil texture and SOM of KL60 (StdDev in Tab. 1). The BNC-based lime requirement (LRBNC) to reach the target pH value of 6.5 ranges between -1,117 (KL60) and 1,484 kg CaCO3 ha−1 (KL60) which illustrates that large sections of the investigated fields, showing negative LRBNC values, do not require any lime fertilization. However, negative and positive LRBNC values occur on the same field underpinning the necessity of site-specific pH management for yield optimization. In accordance with the pHBC, the spatial variability of LRBNC is lowest for KL41 showing a within-field range of only 436 kg CaCO3 ha−1 and highest for KL60 having a range of 2601 kg CaCO3 ha−1. Comparing these findings with the standard soil characterization regarding the most relevant LR affecting soil properties, it is noticeable that KL60, which is situated in the valley of the River Oder, also shows the highest ranges of pH (2.4 units), sand (67%) and SOM content (4%) of all investigated fields. In contrast, for KL41 the ranges of values are lowest with only 1.6 units for pH and 19% for sand and among the lowest with 1.1% for SOM (Table 1).

Since LRBNC also defines negative LR values, it provides useful information on the magnitude of soil acidification necessary when a soil has a too high pH value. On that basis, the farmer is enabled to evaluate if these sections of the fields are simply left out of lime treatment or if a treatment with physiologically or chemically acidic fertilizers may be reasonable to increase soil productivity. However, apart from PP1392, all investigated fields rather require lime fertilization than active acidification.

Sensor-based prediction of LRBNC

Prior to the MLR modeling, the independent variables were tested for interdependencies (Table 4). Highest correlations were detected for OpticMapper Infrared (OM-IR) and OpticMapper Red (OM-Red). Of minor importance was the correlation between shallow apparent electrical conductivity (EC-sh) and OM-IR as well as between pH and OM-IR.

Table 4 Intercorrelation between soil sensor readings illustrated by Pearsson’s R
Full size table

The best performing multi-variate linear regression (MLR) models for LRBNC prediction and their figures of merit are described in Table 5 and Fig. 4. All models received very high to high adjusted R2 (adj-R2) values between 0.70 and 0.91 and low RMSE values of 65 to 234 kg CaCO3 ha−1. This demonstrates that LRBNC can be successfully predicted with the present approach and the proximal sensor technique used. By comparison, an ‘on-the-go’ soil pH and lime requirement measurement system based on a pH ISFET (ion sensitive field effect transistor) sensor tested by Viscarra Rossel et al. (2005), achieved an accuracy of estimated LRs of about 600 kg ha−1. Lund et al. (2005) received an RMSE of 643 kg ha−1 for a 20 ha field in Kansas (USA) using sensor pH, ECa and NIRS data and a locally weighted partial least squares regression analysis.

Table 5 Results of the stepwise regression analysis with forward selection in order to predict LRBNC
Full size table
Fig. 4
figure 4

Correlations between measured and MLR-predicted LRBNC for A GW6, B GW21, C GW32, D KL6, E KL41, F KL42, G KL60, H PP1392, and I PP1401

Full size image

It is striking that the best performing models for the nine study sites show a rather different combination of independent variables. The only concordance of the MLR models is the premier importance of the pH value for determining LRBNC. For two fields (GW6, KL42), the best predictive model was based solely on the pH value. In two cases, two covariates were used, i.e. pH plus EC-sh (KL41), pH plus OM-IR (PP1401), and pH plus the ratio between OM-IR and OM-Red (ratio_OM-IR_OM-Red; KL60). A model of three covariates performed best at three fields using pH plus EC-sh and ratio_OM-IR_OM-Red (PP1392), pH plus the residuals of OM-IR and OM-Red (res_OM-IR_OM-Red) and the residuals of EC-sh and OM-IR (res_EC-sh_OM-IR; GW21) as well as OM-IR plus pH and EC-sh (GW32). One MLR model contained four covariates, i.e. pH, EC-sh, res_OM-IR_OM-Red and the residuals of pH and OM-IR (res_pH_OM-IR; GW32). Finally, a model of five independent variables performed best at KL6 including res_OM-IR_OM-Red, pH, ratio_OM-IR_OM-Red, OM-IR and EC-sh.

The high sensitivity and partly exclusivity of pH in the determination of LRBNC is in contrast with findings of Viscarra Rossel and McBratney (2001) who predicted LR as a function of pH, SOM, clay content and exchangeable Al in south-eastern Australia. They state that soil pH and exchangeable Al explained only moderate proportions of the variation in LRs. Hence, these properties alone do not provide accurate estimates of a soil’s lime requirement.

Considering the generally large effect of the pH value to successfully model LRBNC, it could be argued that pH mapping alone could be sufficient in predicting LRBNC of a field. To evaluate that, univariate linear regression (ULR) models were also developed between LRBNC and the sensor data. Their performance is illustrated in Fig. 5 demonstrating that in eight out of nine fields LRBNC correlates best with the sensor pH value receiving an R2 between 0.59 and 0.87 and a mean of 0.72. By contrast, the univariate correlations with the ECa and OM data were much lower (EC-sh: R2, 0[Min]…0.31[Mean]…0.68[Max]; OM-IR: R2, 0…0.27…0.76; OM-Red: R2, 0…0.16…0.53).

Fig. 5
figure 5

R2 of univariate linear regression models between sensor data and LRBNC (EC-sh: shallow apparent electric conductivity; OM-IR: OpticMapper Infrared; OM-Red: OpticMapper Red)

Full size image

There are two important reasons for the good performance of the sensor pH compared to the other sensor data. One reason is that, even though the sensor output is voltage, the pH sensor very directly assesses the soil pH value due to the selectivity of the antimony electrode (Lund et al., 2005). The relationship between potentiometric reading and pH value is well established and thus the sensor-based pH measurement can be regarded as reliable (Subirats et al., 2015). The second reason is that apparent electric conductivity and optical reflectance, on the other hand, can rather be considered integrative soil parameters, which are affected by a variety of soil characteristics, in particular soil moisture (Corwin & Lesch, 2003; Lück et al., 2009). Compared to the pH electrode, EC readings and optical reflection data are less selective. Thus, the correlation with soil texture and SOM, respectively, which in this study are among the key determinants for lime requirement, can sometimes be low. This is in line with findings by Corwin and Lesch (2003) and Sudduth et al. (2005). Moreover, another possible reason for the poor score of the OpticMapper data is the above mentioned low SOM contents of the fields studied showing a maximum range of 4.8% and standard deviation of 1.4%. Notwithstanding, the potential of the OpticMapper in determining LRBNC was demonstrated at GW21 where OM-IR and OM-Red obtained an R2 of 0.76 and 0.53, respectively and where res_OM-IR_OM-Red ranked highest in the best performing MLR model.

Comparing R2 of the ULR and adj-R2 of the MLR models, it can be seen that the application of MLR increased the performance of the LRBNC models by 3 to 27% with a mean value at 9%. As a consequence, multi-variate sensor data fusion considerably increased model performance compared to the univariate predictions. This is in accordance with Lund et al. (2005) stating for one test field that the LR model improved by 31% when a multi-sensor (pH, ECa, NIRS) instead of a single-sensor (pH) approach is applied.

The different influence of the covariates in the prediction models shows that model performance is field dependent. This emphasizes the necessity of a field-wise calibration of the sensor data and complicates the development of cross-field calibration models. This corresponds to findings of Schirrmann et al. (2012) from the same geographical region.

Regionalization of sensor-based LRBNC maps

The best performing multi-variate calibration models were applied to generate high resolution LRBNC maps of the nine fields of investigation (Fig. 6). The LRBNC can be categorized into three different soil acidity or LR domains: (i) sections with too low pH values that need to be treated with lime (blue colors), (ii) sections that show LRs near zero that are characterized by a pH at the optimum (grey), and (iii) sections where the pH is too high and needs to be lowered to reach the optimum pH of 6.5 and, thus, received negative LRs (red colors).

Fig. 6
figure 6

Field-wise regionalized LRBNC maps of A GW6, B GW21, C GW32, D KL6, E KL41 (left), KL42 (right), F KL60, G PP1401 (left), PP1392 (right). Note the two different scales, Projection: UTM ETRS89 33 N; Aerial Photographs: Google | DigitalGlobe; Bing Maps, Microsoft (Color figure online)

Full size image

From the within-field spatial patterns of LRBNC (Fig. 6) and the results of the semivariance analysis (Table 6), it can be seen that the test fields show a more or less high spatial variability in LR. The magnitude of variability (sill parameter), i.e. the range of values of lime and/or acid requirement, is lowest for KL42 and highest for KL60. In contrast, the spatial context of autocorrelation (range parameter), i.e. the distance in which the variability occurs, is smallest for GW6 (70 m) and largest for 382 m (PP1401). Since a high spatial variability depends on both high sill and low range parameters, the sill-to-range ratio (SRR) can be used in order to quantify spatial variability of LRBNC. SRR is highest at KL60 showing the highest spatial variability, which corresponds to the earlier stated high variability of LR-affecting soil properties. For KL60, high LR variability mainly derives from the highest sill value of all test sites whereas the range is rather high (rank 8). Intermediate SRRs were obtained for GW6, GW21, GW32 and KL6. Low values show PP1401 and PP1392 mainly due to high range values whereas KL41 and KL42 have the lowest SRRs and spatial LR variability of the nine test sites. As a consequence, site-specific variable rate lime fertilization is recommended for most of the test fields in order to increase soil fertility and optimize agricultural productivity. Especially at KL41 and KL42, where LRBNC is relatively homogeneously low, the designation of rather large sub-plots for single-rate lime applications can be justified.

Table 6 Semivariogramm parameters of LRBNC of the test fields and ranking of within-field LRBNC variability
Full size table

The raster histograms in Fig. 7 show the field-wise quantification of spatial distribution of the three above mentioned LR domains. Three fields (GW6, GW21, GW32) have their areal maximum in domain one showing high soil acidity that needs to be managed by lime fertilization. Four fields (KL6, KL41, KL42, PP1401) are dominated by LR domain two that have an optimal pH and need no lime application and two fields (KL60, PP1392) have the strongest areal share in LR domain three showing too high pH values where liming should be omitted. In contrast, for these areas, it could be thought about using fertilizers that, as a result of physiological or chemical reactions, acidify the soil.

Fig. 7
figure 7

Field-wise raster histograms of regionalized LRBNC for A GW6, B GW21, C GW32, D KL6, E KL41, F KL42, G KL60, H PP1392, and I PP1401

Full size image

Conclusions

Spatially varying lime requirements (LR) within nine agricultural fields in northeast Germany were successfully predicted using base neutralizing capacity (BNC) data from laboratory and proximal multi-sensor mappings. Compared to the current best management practices in LR determination, this direct approach has potential to reduce the time and cost of laboratory analyses with a simultaneous increase in the spatial resolution of LR data. The best performing models from stepwise multi-variate linear regression (MLR) analysis received high adjusted R2 between 0.70 and 0.91 and low RMSE values ranging from 65 to 204 kg CaCO3 ha−1. Sensor data fusion increased the model performance by 3 to 27% with a mean at 9%. High resolution LRBNC maps of the nine fields were produced. LRBNC could be categorized into three different soil acidity or LR domains: (i) areas of lower than optimal pH values that need lime treatment, (ii) areas that have a pH at the optimum at which no lime is necessary, and (iii) areas with pH values greater than 7 where liming should be omitted and an estimate of the quantity of chemically acidifying fertilizers to reduce pH are provided. Within-field variability in LR was quantified using the sill-to-range ratio from semivariance analysis for the sensing density imposed on the fields. In seven out of nine prediction models, the sensor pH value was the most important predictor variable. Thus, it might be cost-efficient just to use a pH sensor for determining LR if soil characteristics were similar within a region of fields. However, results and conclusion apply only for this soil-scape. In order to validate these findings for other regions, additional BNC studies should be carried out in different soil-scapes.