## 2 Materials and Methods
### Multivariate Calibration  
Using the ParLeS software (Vers. 3.1) (The University of Sidney, Australia) [36], multivariate calibrations were carried out to predict the investigated soil properties from reflectance spectra. The spectral data were calibrated with reference laboratory soil data through PLSR. Before data processing, spectroscopic data were pre-treated to reduce undesirable variance in the data to improve the predictive capacity of the multivariate calibration models. Before multivariate calibration, all of the common pre-processing methods were compared, either alone or in combination. The following methods were considered: signal correction (MSC) [37], standard normal variance (SNV) correction [38; 39], wavelet detrending (WD) [40], first and second derivative transformation, and median and Savitzky-Golayfiltering [41]. Before multivariate modelling, spectra were also pre-treated by mean centring the data. The abovementioned procedures were applied to the original reflectance spectra (R) and their absorbances (A = log 1/R). The number of components to keep in the calibration models was then determined using leave-one-out cross-validation [42]. We calculated the root mean squared error (RMSE) of predictions to find the best cross-validated calibration model:  
\[\text{RMSE}=\sqrt{\frac{1}{\text{N}}\sum_{\text{i=1}}^{\text{N}}\left(\text {y}_{\text{pred}}-\text{y}_{\text{ob}}\right)^{2}} \tag{1}\]  
where N is the sample size, ypred is the predicted value, and yob is the observed value. Typically, the model with the lowest RMSE is chosen [4]. However, it would be recommended to use a more parsimonious PLSR model (i.e., one with fewer components) that accurately describes data variability without overfitting [4]. To this end, the Akaike Information Criterion (AIC) [43,44] can be used to determine the best factor selection:  
\[\text{AIC}\ =\ \text{log}\ (\text{RMSE})\text{N}\ +2\text{m} \tag{2}\]  
where N is the sample size and m is the number of model parameters (i.e., the number of factors to be estimated). The adjusted coefficient of determination (R2\({}_{\text{adj}}\)) [45] and the relative percent deviation (RPD), i.e., the ratio of the standard deviation of analysed data (that is, soil properties) to RMSE, were used to assess model accuracy [4]. According to previous studies in the literature [46,47], the quality of RPD predictions was classified as follows: RPD less than or equal to 1 designates a very poor predictive model, whose use is not recommended; RPD between 1.0 and 1.4 designates a poor model, where only high and low values can be distinguished; RPD between 1.4 and 1.8 designates a fair model, which can be used for evaluation and correlation; RPD between 1.8 and 2.0 designates a good predictive model, which can be used for predictions; RPD higher than 2.0 designates an excellent model, which can be used for quantitative prediction. Due to the availability of a limited number of samples, it was not possible to validate the model with an independent set of samples [6]. The PLSR \(b\) regression coefficients [48] and the Variable Importance in the Projection (VIP) [49, 50, 51] were computed to find important wavelengths for the prediction of the investigated variables. PLS regression coefficients offer a compact representation of the X-Y (i.e., wavelength--soil contaminants) relationships of the PLS model [45]. VIP scores summarise the impact of individual X variables (wavelengths) on the PLSR model. They are calculated as the weighted sum of squares of the PLS weights (w*), which take the amount of explained y variance in each extracted latent variable (dimension) into account [45]. VIP scores give a valuable measure for determining which variables mostly contribute to the y variance explanation. There will always be only one VIP scores-vector for a given model and dataset, summarising all selected components and y variables [52]. The threshold for the b-coefficients was based on their standard deviation [53].