Introduction

In recent years, visible and near-infrared (visNIR) spectroscopy has been increasingly applied to quantify a broad variety of chemical and physical soil constituents (e.g. Ben-Dor & Banin, 1995; Christy, 2008; Dalmolin et al., 2005; Kodaira & Shibusawa, 2013; Liu et al., 2017; Reeves et al., 2002; Viscarra Rossel et al., 2009). This technique has several methodological benefits for soil analysis such as time and cost saving, little sample preparation without environmentally harmful chemicals and a non-destructive sample analysis. Furthermore, it can address several soil characteristics with a single spectrum (Reeves, 2010; Viscarra Rossel et al., 2006a). However, the latter advantage entails a major challenge of visNIR soil spectroscopy, i.e. the complex nature of spectral data and the manifold and mostly overlapping interactions of molecules with radiation that are not fully understood. Disentangling these interactions requires sophisticated data pre-processing and chemometric analysis methods to correlate specific spectral features with a single soil parameter (Dotto et al., 2017; Viscarra Rossel et al., 2006b; Wight et al., 2016) as well as an appropriate amount of reference soil samples for calibration and validation (Pinheiro et al., 2017).

Lime fertilization is one of the most fundamental management strategies in agronomy to increase soil fertility and optimize crop yields. A precise lime application requires detailed knowledge about the field’s soil acidity and its pH buffer capacity to deduce its lime requirement (LR), i.e. the amount of lime needed to raise the soil pH to an optimum value. This information can be derived directly from conventional field or pot experiments (e.g. soil–lime incubations) or from laboratory analysis (e.g. soil-base titrations). The base neutralizing capacity (BNC) is a soil parameter, which is determined by a soil-base titration where the response of the soil to the addition of increasing concentrations of a basic solution (e.g. Ca(OH)2) is recorded. Subsequently, the soil’s LR can be derived by calculating the amount of Ca(OH)2 used to reach the optimum soil pH (target pH). However, this method is too laborious, time-consuming and expensive to be extensively applied in routine soil testing. In contrast, pedotransfer functions could be used to estimate the lime requirements using the statistical relationship with other soil properties which can be measured more effectively (Merry & Janik, 1999), such as visNIR spectra.

During the last decades, the potential of visNIR and mid-infrared (MIR) spectroscopy to predict soil pH and LR for digital soil mapping was discussed in several publications. Therein, diverse prediction methods such as Multivariate Regression (MVR), Partial Least Squares Regression (PLSR), or Support Vector Machines (SMV) were tested. Due to the complexity of soils and soil spectra, each of these methods had their advantages in certain environments and certain tasks and probably there is no ultimate best method. Viscarra Rossel and McBratney (2008) reviewed studies conducted between 1986 and 2006 that used visNIR spectroscopy to predict soil pH and LR. They averaged the coefficients of determination (R2) for pH and LR and concluded that soil spectra in the NIR region result in higher R2s (pH: Ø of R2 = 0.68; LR: Ø of R2 = 0.62) compared to soil spectra in the visible (vis) region (pH: Ø of R2 = 0.36; LR: Ø of R2 = 0.25). Pinheiro et al. (2017) evaluated the performance of visNIR spectral data for the prediction of soil properties in the Central Amazon (Brazil) using PLSR and a set of 41 pre-processing methods. For soil pH, they obtained a rather low R2 of 0.4, which, they state, corresponds with findings of Terra et al. (2015) where Brazilian soil reactive properties like pH could not be well predicted from the visNIR spectra using support vector machines (SVM). Merry and Janik (1999) analyzed the relationships between pH buffering capacity (pHBC) and related Australian soil acidity properties, such as pH, carbon content, clay content and LR using Fourier transform (FT)-NIR and FT-MIR spectrometers. They found MIR to be superior to NIR in most of the PLSR models. For soil pH and LR, they obtained R2 values of 0.45 and 0.77 for NIR spectra and 0.80 and 0.85 for MIR spectra, respectively. Viscarra Rossel et al. (2001) predicted pH and LR for Australian soil using MIR diffuse reflectance spectroscopy combined with PLSR and obtained an R2 of 0.85 for pH and 0.75 for LR. D’Acqui et al. (2010) used a FT infrared spectrometer in the near-to-mid-infrared spectral range to characterize soils typical for Italian Mediterranean offshore environments. Their best prediction of the soil pH using PLSR showed an R2 after validation of 0.72. Viscarra Rossel et al. (2006a) reviewed a large number of previous studies focusing on the prediction of various soil properties including pH and LR using vis, NIR and MIR spectroscopy. They reported R2s for soil pH ranging from 0.54 (Wijaya et al., 2001; using visNIR and stepwise multiple linear regression (SMLR)) to 0.74 (Reeves & McCarty, 2001; using NIR and PLSR). R2s for LR ranged from 0.73 (Janik et al., 1998; using FT-NIR and PLSR) to 0.86 (Janik et al., 1998; using MIR and PLSR). Their own predictions of LR were moderate using NIR (R2 = 0.50; PLSR) and improved using MIR (R2 = 0.75; PLSR). Leenen et al. (2019) predicted pH and LR using MIR and PLSR for six soil locations (pH: R2 = 0.63–0.93; LR: R2 = 0.47–0.92). Finally, Metzger et al. (2020) published R2s in the range of 0.41 to 0.86 for LR using MIR and PLSR. Additionally, they moved from laboratory to field analysis resulting in a promising R2 of 0.85 for LR using a handheld MIR spectrometer (Metzger et al., 2021). These results indicate that a higher precision and accuracy in predictions of soil pH and LR can be achieved with MIR spectroscopy and PLSR. This is because visNIR spectra show no distinct features that can be unambiguously related to molecules in the sample. Instead, absorption bands strongly overlap, especially in heterogeneous media like soils. Such absorption bands belong to molecular overtone and combination vibrations and occur in the NIR, while MIR spectra capture fundamental molecular vibrations, which are more intensive than overtone and combination vibrations (Sjaunja, 2005). However, MIR spectrometers are still more expensive rendering visNIR spectroscopy a more feasible and affordable alternative.

The main hypothesis (H1) of the present study states that visNIR spectroscopy has the capability to predict LR, for the first time based on lab-analysed BNC, for a target pH value of 6.5 (LRBNC) with particular regard to agricultural fields of a quaternary landscape of Central Europe (North-east Germany). In order to confirm the hypothesis, several data pre-processing techniques and multivariate calibration procedures including new chemometric techniques were applied searching for the best correlation with spectral reflectance signatures. In the few publications using optical spectroscopy to predict LR, except in Metzger et al., (2020, 2021), prediction models were not validated with a new and independent validation set, where R2 values are generally lower compared to the calibration set. Thus, the performances of the prediction models will be validated by an independent data set. Furthermore, it is hypothesized (H2) that important wavelength regions can be identified, which are most relevant for LRBNC prediction. In a novel approach, it is finally tested if the hypothesis (H3) is true that LRBNC can be determined more precisely by predicting LRBNC directly from the soil spectra than by predicting LRBNC indirectly using the predicted BNC parameters.

Materials and methods

Research area

Nine agricultural fields on three farms were studied in a quaternary landscape of the North-east German Plain, which is part of the broader geomorphological region of the North European Plain (Fig. 1). The landscape is the result of repeated Pleistocene glaciations by the continental Scandinavian ice sheet as well as by subsequent periglacial and interglacial Holocene geomorphic processes. In the study area, the landforms and soils were particularly shaped by the advances of the Weichselian (115–12 ka BP) and the preceding Saalian glacial belt (150–130 ka BP; Krbetschek et al., 2008). Climatically, it is situated in a transitional zone between oceanic climate of Western and continental climate of Eastern Europe. Regional climatic differences are rather low due to a relatively low altitudinal range of the land surface of ~ 0 to 200 m a.s.l. Following the Koeppen–Geiger Climate Classification System, the climate of the study region can be classified as temperate oceanic with an increasing influence of continental circulations. The mean annual air temperature is around 9 °C. The coldest and warmest months 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.

Fig. 1
figure 1

Map of Central Europe with the location of the study sites in the federal state of Brandenburg (Germany) [Projection: UTM ETRS89 33 N]. The inlay map indicates the detailed location of the study sites in Brandenburg

Full size image

The three farms are situated in the east and in the north of the federal state of Brandenburg (North-east 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 a high within-field soil variability. Soil acidity of these soils often necessitates liming. However, at some places, the pH is naturally high due to the occurrence of carbonates from glacial till.

Soil samples were collected from the upper 30 cm of the soils’ Ap horizons. In accordance with the German soil classifcation 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) (Vogel et al., 2022). The pH value of the soils range from 3.8 (extremely acidic) to 7.4 (slightly alkaline) with median values of 6.2 (slightly acidic). The soil organic matter content is rather low throughout the study region, having minima of 0.8%, maxima of 5.6% and median values around 1.4% (Vogel et al., 2020).

Analysis of the base neutralizing capacity (BNC)

A total of 170 soil samples were analyzed for the base neutralizing capacity (BNC), which are a partial set of BNC data previously published in Vogel et al. (2020). The BNC is defined as the amount of soil acidity that is neutralized by a base in a given time interval targeting a certain pH value (Meiwes et al., 1984; Vogel et al., 2020). To directly determine the LR of the studied soils from their base neutralizing capacity (LRBNC), the procedure of Meiwes et al. (1984) and Utermann et al. (2000) was followed. It is a laboratory analysis where a base of increasing concentration is added to a soil sample, the resulting pH change is recorded and a titration curve between pH development and the added base concentrations is generated. Then, the BNC of the analyzed soil sample for receiving a target pH value can be deduced from the parameters of the fitted titration curve and finally, the respective LR can be calculated.

In more detail, an air-dried and 2 mm-sieved soil sample is divided into 6 subsamples of 25 g. Afterwards, a control sample is added with 50 ml of deionized water. Furthermore, 25 ml of 2 N CaCl2 and 25 ml of 8 N NaOH solution are added in five concentrations to obtain Ca(OH)2 in the following 6 concentration levels: 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− (Meiwes et al., 1984). After 18 h of overhead shaking, the pH value is measured with a glass electrode (SenTix® 81; WTW) in the supernatant solution. For the quantification of the BNC, the pH values and concentrations of Ca(OH)2 added are displayed in a scatterplot and a titration curve is fitted to the six points by means of a non-linear regression model using the nls function implemented in the free software environment for statistical computing and graphics R (R Core Team, 2018). For the soils studied, Vogel et al. (2020) demonstrated that the non-linear reaction of the pH value on the application of increasing quantities of Ca(OH)2 is best described by an exponential model (Eq. 1):

$${pH}_{target}=\alpha -\beta \cdot {\gamma }^{{Ca(OH)}_{2}}$$
(1)
$${BNC}_{{Ca(OH)}_{2}}=\frac{\mathrm{log}(6.5 - \alpha )+ \mathrm{log}(\beta )}{log(\gamma )}$$
(2)

where α, β and γ are the regression coefficients of the exponential function. The amount of Ca(OH)2 in mmolc (25 g soil)−1 needed to achieve a target pH of 6.5 was derived based on this model (Eq. 2) and the LR, expressed in kg CaCO3 (ha × dm)−1 (dm: decimeter), was calculated by multiplying BNC with the factor 2 000 (Meiwes et al., 1984; Utermann et al., 2000). A pH of 6.5 was chosen as the target pH value because fertilization guidelines, e.g. for the UK (Devra, 2010) and many other countries advice to maintain a soil pH of 6.5 for arable lands to bring most of the nutrients to their optimal availability for plants (Goulding, 2016). Of course, this does not reflect the fact that arable crops differ in their sensitivity to soil acidity. For more detailed information on BNC parameters and how to determine lime requirement by using the BNC, the reader is referred to Vogel et al. (2020).

Spectral analysis

Spectral measurements were conducted with an ultra-broadband UV–Vis-NIR spectrometer system by ArcOptix (Arcspectro UV–Vis-NIR fibered, ArcOptix S.A., Neuchatel, Switzerland). The device is equipped with a dispersive spectrometer using a silicon array detector for the ultra violet (UV) and visual (vis) range (200 to 1000 nm) and a Fourier-Transform spectrometer for the near-infrared region (FT-NIR; 900 to 2500 nm) with an extended range InGaAs photo diode. The spectral resolution is < 5 nm in the UV and Vis and < 8 nm in the NIR region with a reporting interval of 1 nm. Samples were irradiated by four halogen lamps at 45° inside a box which excluded ambient light. The diffuse reflected light was transmitted to the spectrometer via two glass fibres: one high OH fibre for the vis range and another one with low OH bindings specialized for the NIR region. For sample preparation, each soil sample was sieved to < 2 mm, air dried, filled in a petri dish and flattened with a spatula. Measurements were repeated three times in different positions by rotating the sample by about 90 degrees. Each spectrum reported by the spectrometer was the average of 10 internal replicates in the vis range with an integration time of 1100 to 1700 ms and 16 replicates in the NIR using the high gain factor setting. After about every 60 min, the spectrum of a 20% reflection standard (Lake Photonics, Uhldingen-Mühlhofen, Germany) with certified reflection grades was recorded. Furthermore, a dark current measurement was carried out twice a day.

Data modelling

Raw reflection spectra (I) of the soil samples, the reflection standard (I0) with its reference correction values for each wavelength (zλ) and the dark current (Id) were converted into reflectance spectra (R) using Eq. 3:

$${R}_{\uplambda }=\frac{I-{I}_{d}}{{I}_{0}-{I}_{d}}\times {z}_{\uplambda }$$
(3)

Reflection values from 250 to 359 nm were removed due to detector noise. The spectra were then smoothed by a Savitzky–Golay filter with a polynomial first order combined with 11 or 36 windows (Savitzky & Golay, 1964). Afterwards, several pre-processing techniques were tested for optimizing the prediction of the BNC properties. These included:

(i)Pseudo-absorbance transformation (pA; Eq. 4; analogous to Lambert-Beer’s law for transmitted light from 1852),

(ii)Kubelka-Munk transformation (KM; Eq. 5; Kubelka & Munk, 1931),

(iii)Standard normal variate transformation (SNV; Eq. 6; Barnes et al., 1989),

(iv)Multiple scatter correction (MSC; Eq. 7; Martens et al., 1983; Geladi et al., 1985),

(v-vi)Savitzky-Golay first and second derivative (SG1, SG2; Savitzky & Golay, 1964),

(vii)Orthogonal signal correction (OSC; Wold et al., 1998) and

(viii)A technique based on normalized differences (ND) using dual wavelengths indices (DWI; Eq. 8; Schirrmann et al., 2013), where all possible normalized differences of each wavelength are calculated requiring considerable processing power.

Details for SNV, MSC and the Savitzky-Golay derivatives can be found in (Rinnan et al., 2009).

$$pA=\mathit{log}\left(\frac{1}{R}\right)$$
(4)
$$KM=\frac{{(1-R)}^{2}}{2R}$$
(5)
$${SNV}_{w}=\frac{({y}_{w}-\overline{y})}{\sigma }=({y}_{w}-\overline{y})/\sqrt{\frac{\sum ({y}_{w}-{\overline{y })}^{2}}{N-1}}$$
(6)
$$MSC=\frac{{y}_{w}-{a}_{i}}{{b}_{i}}$$
(7)
$${ND}_{ij}=\frac{{R}_{i}-{R}_{j}}{{R}_{i}+{R}_{j }};i<j$$
(8)

where R is the reflectance spectrum, yw the reflectance at a specific wavelength, ӯ the mean value of all measured values of one spectrum, ai the additive effect calculated by ordinary-least-square-regression (OLS) for each spectrum on the mean spectrum of all spectra, bi the multiplicative effect calculated by ordinary-least-square-regression (OLS) for each spectrum on the mean spectrum of all spectra, and Ri and Rj the reflectance values at the ith and jth wavelength of a spectrum. The LRBNC (BNC-based LR determined for a target pH value of 6.5) and the following five parameters were used as dependent variables:

(i) pH0: initial pH value of the soil sample measured in deionized water before base addition,

(ii) δpHtotal: total pH increase over all five base additions, and

(iii-v) α, β and γ: regression coefficients of the exponential model.

To find the best correlation between spectral reflectance signatures and the six dependent variables, six regression methods were applied and the results compared:

  1. (i)

    partial least squares regression (PLSR; Wold, 1975) with a tenfold cross validation using the non-linear iterative partial least squares (NIPALS) algorithm,

  2. (ii)

    canonical powered PLSR (CPPLSR; Indahl et al., 2009),

  3. (iii)

    least absolute shrinkage and selection operator regression (LASSO; Tibshirani, 1996),

  4. (iv)

    least angle regression (LARS; Efron et al., 2004),

  5. (v)

    random forest (RF; Breiman, 2001; Ho, 1995), and

  6. (vi)

    forward stagewise subset selection combined with PLSR (FS-PLSR).

The FS-PLSR preselects wavelengths by an algorithm that belongs to the family of forward stagewise regressions (Sen & Srivastava, 1990) before using PLSR. The algorithm is inspired by sure independence screening (SIS) from Fan and Lv (2014) and it involves the following data processing steps: Firstly, the predictor variable (wavelength) with the absolute maximum of Kendall τ correlation to the response variable is selected (Kendalls τ is supposed to be a robust variant for correlations of nonparametric variables; Fan & Lv, 2014). Secondly, a robust linear model is built between the selected predictor variable and the response variable. Thirdly, the residues of this model are taken as new response variable and the complete procedure is started again until the same predictor is selected twice consecutively. Finally, the set of selected predictor variables is used to build a PLSR model.

All 170 spectra were randomly separated into a training dataset for calibration including a tenfold cross-validation (repeated for 20 times) and an independent test set for validating the model algorithms in a ratio of 3:1. This resulted in 127 samples for calibration and 43 samples for validation. The prediction models were evaluated by the following diagnostic variables: coefficient of determination (R2; Eq. 9), root mean square error of prediction (RMSEP; Eq. 10), (RPIQ; Bellon-Maurel et al., 2010; Eq. 11) and the number of used components.

$${R}^{2}= \frac{\sum_{i}^{N}({\widehat{y}}_{i}-{\overline{\widehat{y} })}^{2}}{\sum_{i}^{N}({y}_{i}-{\overline{y })}^{2}}$$
(9)
$$RMSEP= \sqrt{\sum_{i=1}^{N}\frac{({\widehat{y}}_{i}-{{y}_{i})}^{2}}{N}}$$
(10)
$$RPIQ= \frac{Q3-Q1}{RMSEP}$$
(11)

where N is the number of samples, \(\overline{y }\) is the mean of all reference values, \(\overline{\widehat{y} }\) is the mean of all predicted values, yi is the ith reference value, \({\widehat{y}}_{i}\) is the ith predicted value, Q1 the first quartile of all reference values, and Q3 the third quartile of all reference values.

The ratio of performance to interquartile range was displayed to describe the relationship between the spread of the data and error of prediction, e.g. an RPIQ of 2 means that the spread of 50% of the data around the median is twice the root mean square error of prediction, thus the higher RPIQ the better the prediction performance (Bellon-Maurel et al., 2010).

Results and discussion

The descriptive statistics of BNC parameters for all studied soil samples are shown in Table 1. The initial soil pH values in deionized water (pH0) varies from 4.5 to 8.0 with a mean value of 6.6 and a strong skewness of − 0.80 towards higher pH values. pH0 is strongly negatively correlated with BNC-based lime requirement (LRBNC) as given by a Pearson’s correlation coefficient of r = − 0.95 (Table 2). When pH0 is higher than the target pH value of 6.5, LRBNC becomes negative implicating that there is no need for lime treatment. With higher pH values it might even be considered to use acid fertilizers or other soil amendments to lower the pH value (Vogel et al., 2020). As the mean value of pH0 of 6.6 is close to the target pH, the mean value of LRBNC for all samples is close to zero. Moreover, the number of negative LRBNC values are higher than the positive ones because of a skewness of over 0.6. This means that there are less acidic soil samples than basic soil samples in the present dataset, i.e. 60 soil samples (35%) are below the target pH value and 110 (65%) are above.

Table 1 Statistical overview of soil parameters
Full size table
Table 2 Internal correlations of soil parameters using the Pearson’s correlation coefficient (r)
Full size table

The regression coefficient α represents the maximum pH value at the endpoint of the soil-base titration (pH5) having a mean value of 12.4. It varies less than all other parameters showing a coefficient of variation (CV) of only 2.6%. In contrast, β, which is the difference from the starting point (pH0) to the end point (pH5) of the titration curve, has a mean value of 5.4. Thus, pH0 should be equal to α–β. However, small differences between experimentally determined titration points and the fitted exponential function can occur (Vogel et al., 2020). There is strongly negative correlation of r = − 0.95 between pH0 and β (Table 2). The third regression coefficient γ varies between 0 and 1 and refers to the curvature of the titration curve, whereas 1 results in a linear and near 0 in an almost right-angled curve. In this sample set, γ varies more than all other BNC parameters showing a coefficient of variation of almost 25% (Table 1). This indicates a high variability in the pH buffer capacity of the investigated soils (Vogel et al., 2020).

The BNC parameters pH0, δpHtotal, β and LRBNC are strongly correlated with absolute Pearson’s r coefficients between 0.82 and 0.95 (Table 2). Although the correlation of γ to pH0, δpHtotal, α, β and LRBNC is less pronounced, a higher correlation to single titration points (pH0.5, pH1.25 and pH2.5) exists. Parameter α shows the lowest correlations. It only exhibits closer relationship with β (r = 0.79).

The descriptive statistics of the 170 reflectance spectra used in this study are presented in Fig. 2. The samples show a classical spectral signature of mineral soils in the visNIR region from 360 to 2450 nm with reflectances up to 50%. Three absorption bands can be recognized in the region of 1400, 1900 and 2200 nm. According to Bishop et al. (1994), the absorptions around 1400 and 1900 nm belong to combination vibrations of water bound as hydrated cations in the interlayer lattices or water adsorbed on particle surfaces. Moreover, the absorption at around 1400 nm can also occur due to overtone O–H stretch vibrations, as they exist in e.g. octahedral lattices of clay minerals. Absorptions near 2 200 nm result from O–H stretch combination vibrations or overtone bending vibrations of aluminium hydroxides (Al–OH) as they occur in kaolinite, illite and montmorillonite (Stenberg et al., 2010).

Fig. 2
figure 2

Descriptive statistics of 170 soil reflectance spectra from 360 to 2450 nm. Min and Max: minimum and maximum of all reflectance spectra; Q1 and Q3: first and third quartile of all reflectance spectra; LR-Min: spectrum of the sample with the lowest LR; LR-Max: spectrum of the sample with the highest LR

Full size image

In Fig. 2, a relationship between the original reflectance spectra and LR is not visible. In contrast minimum and maximum LR show only a small difference in reflectance. However, from the concavity of the spectral curves (Fig. 3) highlighted by the 2nd derivative, an inverse correlation at around 1914 nm can be recognized showing a more pronounced concavity of low LRs compared to high LRs. This is mathematically confirmed, because the 2nd derivative was the best performing pre-processing method for predicting LR with 1914 nm being the most important wavelength (Tables 3 and 4).

Fig. 3
figure 3

Second derivative of 170 soil reflectance spectra, coloured by the intensity of LR

Full size image
Table 3 Best predictions of soil parameters in the validation set
Full size table
Table 4 Importance ranking of different wavelength regions for the prediction of the BNC parameters
Full size table

The model performances of all examined pre-processing and regression methods are shown in Appendix 1. Table 3 displays the best-performing prediction models for each soil parameter including the pre-processing methods and the evaluation parameters (R2, RMSEP, RPIQ and the number of components used). The frequently published ratio of performance to deviation (RPD) with its classifications for quality control of prediction models (Chang et al., 2001) was not calculated because of its equivalence with R2 (Minasny & McBratney, 2013) and its inappropriateness for skewed parameters (Bellon-Maurel et al., 2010). Best predictions for four out of six parameters (LRBNC, β, pH0, δpHtotal) were achieved by the regression model FS-PLSR. As pre-processing method, the second Savitzky-Golay derivation (SG2) produced best results for three out of six models (LRBNC, pH0, δpHtotal). PLSR achieved best prediction performances for α and γ. The best pre-processing method for the BNC curve parameter β was Kubelka–Munk (KM) transformation combined with normalized differences (ND) and for γ the standard normal variate (SNV) transformation. For α, no pre-processing method optimized prediction performance. However, the removal of two outliers out of 43 samples in the validation set improved R2 for α from 0.49 to 0.67. R2 in the validation set for all parameters were in the range of 0.67 (α) to 0.82 (δpHtotal) demonstrating the good model performances. Except for α and γ, the predictive power explained by RPIQ was good for of all models with RPIQs higher than 2.4. The reason, why α has the lowest R2 and RPIQ could be attributed to the generally poor correlations to the other BNC parameters. Although the R2 of γ is good with 0.72, the RPIQ is lower due to a centred distribution of the reference values which is responsible for a smaller interquartile range (Q3-Q1). The prediction models required a small number of PLSR-components from 5 (pH0, δpHtotal, α) to 8 (β) with one exception for γ, which needed 16 components. In general, a high number of components can include too much noise of the data set and thus might lead to overfitting when applying the model to independent data (Gowen et al., 2011). However, in this case, a reduction by e.g. 5 components to 11 components would decrease R2 from 0.73 to 0.62 and is thus not recommended.

From the current knowledge, the soil pH value, i.e. the soil’s proton concentration, does not seem to be directly spectrally active but it is suspected to correlate with other parameters that show spectral response such as soil organic matter (SOM) content or clay minerals (Chang et al., 2001). However, there could be a problem of predicting soil properties being only correlated to spectrally active components because the type and intensity of these correlations may differ from site to site (Stenberg et al., 2010). Nevertheless, comparing the prediction results for the initial pH value (pH0) of the examined nine study fields with previous studies, a validated R2 of 0.76 and an RMSEP of 0.28 can be considered quite good since the published R2 values for pH at an equal scale (country or state scale) vary from 0.55 to 0.77 with an RMSEP of a third to half a pH unit (Stenberg et al., 2010). Only at a smaller spatial scales like farm or field scale, literature results are sometimes better because of less variations in soil characteristics receiving R2s from 0.54 to 0.92 with RMSEPs from 0.17 to 0.31 (Stenberg et al., 2010).

Figure 4 shows the plotted validation results for predicted and measured values of LR and the five BNC parameters. For the LRBNC, the two best direct modelling methods with two different prediction results are displayed (model A, Fig. 4A and model B, Fig. 4B). Although R2 is almost the same (0.67 vs. 0.68), the RMSEP of the predictions calculated with FS-PLSR and SNV pre-processing (model A) is with 144 kg CaCO3 (ha × dm)−1 much better than the other one calculated with FS-PLSR and SG 2nd derivation (model B) with 221 kg CaCO3 (ha × dm)−1. One reason is that the slope factor of model A is with 0.59 poorer than of model B having 0.91. This results in a smaller range of the predicted values of model A compared to model B. A comparison of the ratios of slope and RMSEP shows that they are the same for the two methods (0.65). This indicates that both models can be considered of equal quality if the prediction results of model A are corrected by the slope factor. Nevertheless, model B was regarded as best-performing prediction model because it avoids slope factor correction. Comparing the performance of direct prediction of LRBNC with the indirect prediction using the BNC parameters α, β and γ (Fig. 4C), the indirect method achieved slightly higher correlation coefficients of 0.75 versus 0.68, respectively. This is the opposite behaviour as one would expect, because the accumulated error of the three predicted variables α, β and γ was assumed to be higher than the error of the directly predicted variable LRBNC. Nevertheless, with the present data, it is not possible to give a general preference of either the direct or the indirect prediction method due to the relatively small number of samples in the validation set. This means that the two methods may just mathematically vary around the same expected value.

Fig. 4
figure 4

Validation results for predicted and measured values for A LRBNC (directly predicted, Model A); B LRBNC (directly predicted, Model B); C LRBNC (indirectly predicted by BNC parameters α, β, and γ); D initial pH value of the soil sample (pH0); E total pH increase (δpHtotal); F,G,H regression coefficients of the exponential model α, β and γ including 1:1-line (grey, regression line (black) as well as R.2, RMSEP, bias and slope of the regression model)

Full size image

Figures 4D–H show the regression lines of the five BNC parameters. The reader should be aware that for alpha, two outliers, which were removed for modelling, are still displayed in the plot (Fig. 4F). The overall bias (intercept of the regression line and representative of the systematic error) for the best predictions of all BNC parameters is low and the slope factor is good with values close to 1 (0.79 to 0.95).

The influence of each individual wavelength by variable importance projection (VIP) scores on the PLSR prediction models, as well as the selected wavelengths by forward stagewise subset algorithm for each response variable are displayed in Table 4. VIP scores for wavelengths higher than 1 have a higher influence on the prediction results than the average. The prediction model for α mainly used wavelengths (i) < 500 nm, (ii) 1 000–1 200 nm and (iii) at the water band around 1900 nm, which were also important for δpHtotal, γ, and LRBNC. The model for γ mainly used wavelengths (i) 500–700 nm (ii) at the water band around 1400 nm, and between 2250 and 2450 nm, which were also particularly important for β. The region between 610 and 710 nm was important for all BNC parameters except α. Further less important wavelengths with VIP-scores can be seen in Appendix 2. The attempt to assign the important wavelengths to optically active components was not accomplished in this study as the potential compounds of soil organic matter or clay minerals have many different, wide and strongly overlapping vibrations, such as first, second, or third overtone vibrations as well as combination vibrations (e.g. of O–H, C–H, CH2, CH3, C = C, C-O, C = O, S–H, N–H etc.) including Fermi resonance shifts (Fang et al. (2018). Finding corresponding vibrations for the important wavelengths could be the scope of further research.

In the literature, LR predictions by spectroscopic methods are relatively rare. Reported R2 values are 0.73 (Janik et al., 1998), 0.77 (Merry & Janik, 1999), and 0.50 (Viscarra Rossel et al., 2006a, 2006b) using NIR as well as 0.75 (Viscarra Rossel et al., 2001 republished in Viscarra Rossel et al., 2006a, 2006b), 0.86 (Janik et al., 1998), 0.45 to 0.92 (Leenen et al., 2019), and 0.41 to 0.76 (Metzger et al., 2020) using MIR spectroscopy. However, in contrast to the present study, except for Metzger et al. (2020), they did not validate their models with a new and independent validation set, where R2 values are generally lower compared to the calibration set. In this regard, an R2 of 0.68 or 0.76 obtained in this study can be considered as good model performance.

An important aspect, which was neglected in this study, is the error of the laboratory method for determining BNC data. Janik et al. (1998) determined LR twice for 224 samples using a soil–lime incubation method based on a 14-day incubation (Richards, 1992) similar to this procedure. They found rather low replication R2 values of 0.73 for the slope and 0.85 for the intercept of the titration curves. Currently, it is unclear whether prediction errors of visNIR spectroscopy for LR are considerably different from laboratory errors. Hence, this topic needs further research.

Conclusions

VisNIR spectroscopy (visNIRS) was successfully used to predict the soil’s lime requirement and five parameters of the base neutralizing capacity by testing seven pre-processing techniques and six different multivariate regression methods on 170 soil spectra of nine agricultural fields, thus confirming the main hypothesis H1. The prediction performance in terms of a validated R2 ranged from 0.67 (α) to 0.82 (δpHtotal). In view of the presumption that the soil pH value does not seem to be directly spectrally active, this results indicate a good model performance. The good correlation seems to be attributed to other parameters that show spectral response such as soil organic matter (SOM) content or clay minerals. For four out of six response variables (LRBNC, β, pH0, δpHtotal), best predictions were obtained by a forward stagewise subset selection combined with PLSR as regression model whereas for three out of six models (LRBNC, pH0, δpHtotal) the second Savitzky-Golay derivation was the best pre-processing method. Important wavelength regions could be identified between 300 and 500 for α, between 610 and 710 nm for all BNC parameters except α, around 1400 for γ, around 1900 nm for δpHtotal, α, γ, and LRBNC, and between 2250 and 2450 nm for α and β. Hence, hypothesis H2 was confirmed. Even though the indirect prediction of LRBNC (by using the predicted BNC parameters α, β and γ) performed slightly better (R2 = 0.75) than the direct one (R2 = 0.68), it was not possible to generally prefer one of the two methods. Thus, hypothesis H3 could not be confirmed nor refuted. Field-dependent prediction models may improve the accuracy in comparison to field-independent models as presented in this study. However, practical agriculture demands rather general field-independent prediction models to reduce total soil mapping costs. The effect on the prediction performance of iteratively excluding one or more agricultural fields from the modelling data was not examined in this study due to a low number of samples. However, it is the subject of current research activities.

It can be concluded that visNIRS is a fast and cheap method for predicting lime requirements. Hence, much more field samples can be investigated compared to the standard lab method. This can make visNIR spectroscopy a very efficient method for soil acidity management particularly in combination with precision agriculture applications enabling a within field spatial analysis. Regarding its applicability for a site-specific soil acidity management, in a next step, visNIR spectroscopy will also be tested in the field on moist soils using an on-the-go sensing system. This transfer from lab to field is also encouraged by the promising results of Metzger et al. (2021) using a handheld MIR spectrometer to predict LR on field-moist samples.