Guidelines for precise lime management based on high-resolution soil pH, texture and SOM maps generated from proximal soil sensing data

Soil acidification is caused by natural paedogenetic processes and anthropogenic impacts but can be counteracted by regular lime application. Although sensors and applicators for variable-rate liming (VRL) exist, there are no established strategies for using these tools or helping to implement VRL in practice. Therefore, this study aimed to provide guidelines for site-specific liming based on proximal soil sensing. First, high-resolution soil maps of the liming-relevant indicators (pH, soil texture and soil organic matter content) were generated using on-the-go sensors. The soil acidity was predicted by two ion-selective antimony electrodes (RMSEpH: 0.37); the soil texture was predicted by a combination of apparent electrical resistivity measurements and natural soil-borne gamma emissions (RMSEclay: 0.046 kg kg−1); and the soil organic matter (SOM) status was predicted by a combination of red (660 nm) and near-infrared (NIR, 970 nm) optical reflection measurements (RMSESOM: 6.4 g kg−1). Second, to address the high within-field soil variability (pH varied by 2.9 units, clay content by 0.44 kg kg−1 and SOM by 5.5 g kg−1), a well-established empirical lime recommendation algorithm that represents the best management practices for liming in Germany was adapted, and the lime requirements (LRs) were determined. The generated workflow was applied to a 25.6 ha test field in north-eastern Germany, and the variable LR was compared to the conventional uniform LR. The comparison showed that under the uniform liming approach, 63% of the field would be over-fertilized by approximately 12 t of lime, 6% would receive approximately 6 t too little lime and 31% would still be adequately limed.


Introduction
The productivity of agricultural soils is highly controlled by their acidity and buffering capacity. Soil acidity results from the release of H + from dissolved and solid acids to form H 3 O + ions in the soil solution and is measured as pH. Soil acidity is a key factor in soil fertility that concurrently influences several yield-relevant soil properties, such as: (i) Nutrient availability (particularly P) and pollutant mobility (especially Al, Mn, Cd) (Dahiya and Singh 1982;Goulding and Blake 1998;Gray et al. 2006), (ii) Nutrient utilization and use efficiency (particularly N) (Ahmad et al. 2016;Edmeades et al. 1986), (iii) Biological activity (Cheng et al. 2013;Ekenler and Tabatabai 2003;Larink and Joschko 2014;Stöven and Schnug 2005), (iv) Soil humus content and type (Briedis et al. 2012;Haynes and Naidu 1998;Paradelo et al. 2015), (v) Soil structure, porosity and aggregate stability (aeration, water availability, root growth) (Fiedler and Bergmann 1955;Hartge 1959;Schachtschabel and Hartge 1958), and (vi) Water infiltration, water storage and soil erosion (Ahn et al. 2013;Cuisinier et al. 2011;Horsnell 1984).
For these reasons, farmers strive to obtain and maintain an optimal soil pH to improve crop growth in their fields (Tunney et al. 2010). As soil acidification is a pedogenetic process in humid climates, more protons (H + ions) are added or liberated by precipitation and internal soil processes over time than the soil is able to neutralize (Fujii et al. 2012;Blume et al. 2016). The physico-chemical processes that are relevant to acidification include the dissociation of carbonic acids, the atmospheric deposition of acidic gases and/or acidic precipitation, microbial respiration and/or root exudates, oxidation reactions and the formation of organic acids and anthropogenic activities, e.g., fertilization, or the removal of alkalis by harvesting crops (Holland et al. 2018;Goulding 2016). Hence, in soils that do not contain geogenic carbonates, farmers need to apply lime to their fields to maintain soil fertility.
However, even in countries with intensive agricultural production, such as Germany, the soil pH of agricultural fields is often not within the optimum range. According to a recent national soil pH survey by Jacobs et al. (2018) in Germany, only 35% of the arable soils and 24% of the grassland soils were in the optimum range, whereas the pH of approximately 42% of the mineral soils under arable farming and 57% of the grassland soils was too low. Apparently, lime management on farms in Germany is not sufficient. One reason is that most farmers do not manage soil heterogeneity at the field scale. They try to avoid (i) the additional effort required for soil sampling, (ii) the uncertainties concerning the interpretation of soil information and fertilization decision making and, (iii) the problems related to the availability and use of appropriate fertilizer application technology.
Since crops vary in their tolerance to soil acidity, the optimum pH at which maximum yields are achieved ranges between 5.3 and 6.6 (Goulding 2016). Below this range, yields of crops with high lime demand may decrease by approximately 20-40% (Holland et al. 2018;Kerschberger 1996;Kerschberger and Marks 2007;Manna et al. 2007). Hence, the main goal of liming is to reduce the total acidity of a specified soil volume (e.g., the plough layer) by increasing the pH value to a target value that is optimal for crop growth (Sims 1996). In contrast, pH values that are too high may also have negative effects on nutrient 1 3 availability and reduce crop yields by 5-10% (von Wulffen et al. 2008). To determine the lime requirement (LR) of a soil to achieve its target pH value, several practical techniques have been developed. The most commonly used LR tests are as follows: i Soil-lime incubations involving increasing rates of liming material applied to a fixed quantity of soil, equilibration for a certain duration and deriving a lime-response curve from the pH changes, ii Soil-base titrations with the titration of a soil suspension with a basic solution (e.g., Ca(OH) 2 or NaOH) (McLean 1978;Alley and Zelazny 1987) and pH measurement after a certain equilibration time, followed by the conversion of the added basic solution into a lime requirement, iii Soil-buffer equilibration (the most widespread approach in the USA), adding a chemical buffer solution to a soil sample, allowing them to equilibrate and measuring the buffer pH decrease to assess the amount of soil acidity to be neutralized by liming (McLean 1978), and. iv Estimates based on algorithms developed in empirical studies that use soil pH and other soil properties such as soil texture, soil organic matter, soil type or CEC as indicators of the soil carbon buffer capacity; this method is mainly used in the UK and Germany.
In this study, an empirical algorithm (LR test type iv, above) was used as a standard and adapted to precision farming by including mappings from proximal soil sensors. The empirical algorithm was developed by the Association of German Agricultural Investigation and Research Institutions (VDLUFA) and has been established as the best management practice for liming in Germany (von Wulffen et al. 2008). The procedure is based on 30 years of fertilization trials studying the correlation between soil pH and agricultural yield, brought into a simplified management structure (Kerschberger 1996;Kerschberger et al. 2000;Kerschberger and Marks 2007). The approach involves two steps: (i) a soil sampling of one mixed soil sample that is composed of several sub-samples from either the whole field or from sub-plots of 3-5 ha of assumed soil homogeneity and (ii) a look-up table system that defines the target pH value for the management unit from the analysed soil texture, soil organic matter (SOM) content and the current pH value (Methods). However, the VDLUFA guidelines for liming are limited because they are based on relatively rough classifications of soil texture and SOM into five and four classes, respectively. However, the algorithms that are needed in the context of the present-day requirements of precision farming should be continuous and stepless. Furthermore, site-specific and variable-rate liming (VRL), which is a precondition for optimizing soil acidity management, requires soil data at a very fine spatial scale (von Cossel et al. 2019). High-resolution maps can therefore help to assess internal field variations in soil properties and reduce the decision uncertainty caused by this unknown spatial variation (Schellberg et al. 2008;Zhang et al. 2016). Various soil proximal sensors are available that can provide information on relevant input parameters for lime requirement calculations, including geo-electrical and gamma-ray sensors for soil texture, optical sensors for organic matter content and ion-selective pH electrodes for pH values (Adamchuk et al. 2018;Gebbers 2018).
Most sensors do not measure the soil property of interest directly but provide readings from a proxy that can be related to the soil property of interest by analysing reference soil samples and establishing statistical models. Sensors for measuring electrical resistivity (ERa) and its reciprocal, bulk electrical conductivity (ECa), are commonly 1 3 used for mapping soil properties that are affected by soil texture, water content and bulk density as well as by mineralogy, porosity, salinity, temperature and organic matter (Corwin and Lesch 2005). The natural variation in total γ-activity in soils is mainly related to the decay of K, U and Th isotopes. Since K is usually associated with clay minerals, γ-activity is a good indicator of clay content and soil texture. Compared to the spatial variations in ERa (ECa), the spatial variation in soil moisture has little effect on γ-activity. Thus, a multiple-sensor approach combining electrical and γ measurements can improve the determination of soil properties (Castrignano et al. 2012;Mahmood et al. 2013). Optical sensors that obtain visible and near-infrared (Vis-NIR) spectra can provide information on soil properties such as the clay, iron oxide, SOM content and carbon mineralogy (Rossel and Chen 2011), and electrochemical sensors that use ion-selective membranes can detect the activity of ions such as hydrogen, potassium or nitrate (Gebbers and Adamchuk 2010;Adamchuk and Viscarra Rossel 2011).
However, the successful adoption of these systems in practice is often hindered by the lack of knowledge on (i) how the sensors work and how reliable they are, (ii) how the sensor data should be calibrated, and (iii) how the sensor data should be further processed to produce site-specific liming recommendations that are in line with best management practices. These questions are related to the scientific foundations of measurement principles, soil buffering, technical possibilities and restrictions, and socio-economic aspects, including cost efficiency and official regulations.
Moreover, only a few studies have compared VRL with conventional approaches (Borgelt et al. 1994;Zaman et al. 2003;Bianchini and Mallarino 2002). For North American soils, Borgelt et al. (1994) found that mean liming rates would have resulted in over-fertilization of 9 to 12% and under-fertilization of 37 to 41%, whereas Bianchini and Mallarino (2002) found that much less lime (56-61%) needed to be applied with the VRL approach. In a similar study in the UK, Zaman et al. (2003) found that 35% of the tested field required more than the average liming rate, 56% required less and only 9% was adequately limed. However, none of these studies used high-resolution soil maps based on proximal soil sensing. This kind of sensor-based approach was explored by Kuang et al. 2014Kuang et al. , 2015. They used on-the-go visible and near infrared (vis-NIR) spectroscopy sensors and two statistical methods (artificial neuronal networks and partial least square regression) to generate high-resolution SOC, pH and clay content maps as inputs for VRL on two fields in Denmark. For one of these fields, Kuang et al. 2014 compared sensor based VRL with uniform liming and observed increase in spring barley yields under VRL. However, Kuang et al. 2014 used a high number of soil reference samples (132 samples on 18 ha) and the proximal soil sensing system was operating at a slow speed of 2 km/h. This might not be accepted for practical farm management. Lime recommendations were calculated by an algorithm from the "Danish Centre for Food and Agriculture", but no bibliographical references or other details were provided.
With practical application in mind, the overarching objective of this paper is to provide guidelines/a protocol for deriving high-resolution lime recommendation maps from the following mobile proximal soil sensor systems: pH electrode, electrical conductivity, gamma ray and optical dual wavelength systems. The specific objectives were (i) to test different proximal soil sensors and sensor combinations to predict the target parameters of soil pH, texture and soil organic matter (SOM) content, (ii) to apply an adapted and currently well-accepted lime recommendation algorithm to the demands of site-specific acidity management and (iii) to compare the results of the novel variable-rate liming approach with a uniform-rate liming strategy developed with the conventional protocol.

Workflow for producing the variable lime requirement maps
To produce the variable lime requirement maps, extensive guidelines were established, including the proximal soil sensing as well as the whole data processing method, from generating maps of soil pH, texture and SOM to the calculation of the precise lime demand to the aggregation of the data to potential working widths (Fig. 1). All data processing and statistical analyses were carried out in the free R software environment for statistical computing and graphics (version 3.6.1) (R Core Team 2018).

Applied on-the-go sensors for generating high-resolution maps of pH, texture and SOM
In this study, the non-commercially available Geophilus measurement system (Lück and Rühlmann 2013) and the commercially available Veris MSP (VERIS Technologies, Salinas, KS, USA) ( Fig. 2) to generate high-resolution soil ancillary data and subsequent predictions of the parameters soil pH, texture and SOM.

Geophilus measurement system
The Geophilus system is merely built for scientific purposes and includes a multi-depth electrical resistivity sensor and a gamma ray sensor (Lück and Rühlmann 2013) (Fig. 2a). The Geophilus system consists of seven pairs of rolling electrodes. One pair directs an electrical current into the soil, and the other six pairs measure the voltage drop. Electrical resistivity (ERa) is explored at six depth levels, from the soil surface to a depth of investigation of 1.5 m. The γ-sensor measures the soil-borne γ-radiation activity as the total counts per second (cps) in approximately the upper 0.3 m soil layer. The system logs the sensor values each second along with the co-ordinates tracked with a differential Global Navigation Satellite System (dGNSS). When mapping with a typical speed of 10 km/h, the sampling interval is approximately 3 m. When the distance between the tracks is 18 m (Fig. 3a), approximately 200 data points are measured per hectare.
The Geophilus system enables the fusion of sensor data to produce additional information. Because the γ-radiation is less sensitive to soil moisture than the ERa readings, the ratio between the γ-activity and the ERa of the array with the smallest electrode spacing (investigation depth: 0-0.25 m) represents the influence of the soil water on the ERa readings. This ratio is expressed as the dimensionless soil water index (SWI): where SWI increases with increasing soil moisture.

Veris multisensor platform
Currently, there is only one commercially available automated on-the-go pH sensor system. The Soil pH Manager™ is part of the Veris MSP (Fig. 2b, c) and was developed based on the work of Adamchuk et al. (1999), and described and applied by Lund et al. (2005) and Schirrmann et al. (2011a, b).

Soil pH manager
The pH value was measured on-the-go by two ion selective antimony electrodes on naturally moist soil material. While driving across the field, a sampler was lowered into the soil to approximately 0.12 m depth, and the soil flowed through the sampler's orifice. When the soil sampler was raised out of the soil, the soil inside the sampler was pressed against the two ion-selective antimony electrodes. Measurements were recorded if they were sufficiently stable within a maximum time of 20 s. A logger recorded the raw potential data along with the dGNSS co-ordinates. Additionally, an online conversion (1) SWI = ERa ⋅ 100 Fig. 2 Applied soil sensing platforms: a The non-commercially available Geophilus system with 7 rolling electrode pairs (1) and a γ probe (2), b the commercially available Veris mobile sensor platform (Veris MSP) by VERIS Technologies with the ERa instrument (3), c the Soil pH Manager™ (water tank (4), soil sampler (5) and pH electrodes (6)) and d the OpticMapper (opening coulter (7) and optical module between depth-sensing side wheels (8)) of the voltage data into pH values was conducted based on a preceding calibration with pH 4 and 7 standard solutions. After each measurement, the sampler was pushed into the soil again, and the old soil sample was replaced by new material that entered the sampler trough. In the meantime, the pH electrodes were cleaned with tap water from two spray nozzles to prepare them for the next measurement cycle. Typically, pH values were recorded every 10-12 s. Geographic co-ordinates were recorded when the sampler shank was raised out of the soil. This sensor can be operated at an approximate speed of 7.5 km/h. With measurements taken every 10 s and a track distance of 12 m, approximately 30 measurements per hectare can be obtained (Fig. 3b). After calibration, the estimated total error of the soil pH maps is less than 0.3 pH (Adamchuk and Lund 2008). In addition, ERa is measured by the sensor platform at a rate of 1 Hz with a galvanic coupled resistivity instrument using six parallel rolling coulter electrodes. This electrode configuration provided readings from two depths with a median depth of exploration of 0.12 and 0.37 m, and the data are expressed as the apparent electrical conductivity (ECa) (Gebbers et al. 2009).

OpticMapper
The soil organic matter (SOM) content was estimated using data generated from the Optic-Mapper (Veris Technologies, Salinas, KS, USA) ( Fig. 2d). It is an on-the-go optical soil sensor that basically consists of a single photodiode and two light sources (LED) that enable reflectance measurements at 660 nm (red) and 940 nm (near-infrared NIR), each with a bandwidth of 20 nm. According to  and , absorption at these two wavelengths is particularly sensitive to organic matter content. At the front, the OpticMapper has an opening coulter that cuts crop residues. The optical module is mounted on the bottom of a furrow 'shoe' between two side wheels that control the sensing depth. The wear plate is pressed against the bottom of the furrow approximately 0.04 m below the soil surface with a consistent pressure to provide a self-cleaning function.
Light is emitted alternately from the two LEDs and passes through a sapphire window onto the soil. The reflected light is captured by a photodiode, and the light intensity is stored in dimensionless values. The digital reflectance data and GNSS co-ordinates are recorded at a rate of 1 Hz. At a speed of 10-12 km/h and 12 m track distance, an average of 260 reflection data points per hectare can be collected (Fig. 3c).

Test field
The selected test field is part of the farm Komturei Lietzen (KL) and is located approximately 40 km east of Berlin (Germany) in the eastern North German Plain (5831100N, 450100E; UTM ETRS89 33N). While the Geophilus system and the Soil pH Manager were applied in April 2018, the OpticMapper campaign was conducted in August 2018. Records were taken along the field's working tracks. This driving path caused fewer errors than driving against the actual working tracks and addressing their spatial irregularities.
The soils at the agricultural study site developed on morainic landscapes shaped by the Pleistocene glaciation processes as well as by fluvial processes in the river valley of the Oder River (Krbetschek et al. 2008). The patterns of spring wheat yield from 2018 therefore reflect the natural geological conditions of the current test field (Fig. 3d). In accordance with the German soil classification system KA5 (Eckelmann et al. 2005), the soil textures at the study site range from pure sand (class: Ss) to loamy sands (classes: St2, Su2, Sl2 and Sl3) and loams (classes: Lt2, Ls2, Ls3 and Ls4). Hence, the soil cover is highly heterogeneous at the test site and in the selected field (field 60 of KL, henceforth called KL60) and is therefore a good example for demonstrating the potential of proximal soil sensing for site-specific lime management. Climatically, the test site is located in the transition zone of the humid oceanic and dry continental climates. The annual mean temperature is 9 °C, and the total mean annual precipitation is 550 mm.

Data cleaning and pre-processing
Before data interpolation, the raw sensor data were observed visually in advance (e.g., for points with strong deviations from the surrounding observations), and obvious measurement errors were removed if necessary. These errors may occur due to insufficient sensor connectivity to the soil or recording issues related to the handling of the sensor platforms.
For example, the OpticMapper still records measurements while the sensor shoe is being lifted out of the soil. Thus, records from residual soil -for example, those taken while driving from one tracking line to the next at the end of the field -were removed from the overall data set. To avoid errors in building the covariance matrices used in kriging, observations that shared identical spatial locations were identified, and duplicates were removed in advance.

Variography
The theoretical semivariogram models were fitted as global variograms to the empirical semivariograms, which provided the spatial weighting function for the subsequent kriging interpolation. The empirical semivariogram calculations were performed by selecting robust variogram estimates to prevent effects from extreme outliers (Cressie 1993). The theoretical variograms were additionally fitted with localized cut-offs to meet the criteria of obtaining good fits at distances smaller than the whole plot. Furthermore, the model fitting was performed by weighted least squares approximation (fit method 7 in gstat), dividing the number of pairs in one bin by the square of the bin's metric distance (Pebesma 2004). After an initial fit of the semivariogram model, a leave-one-out cross-validation procedure was applied (Webster and Oliver 2007) using the initial semivariogram model to predict the values by ordinary kriging at each measurement location after excluding the sample value at that particular point.

Interpolation
Two geostatistical methods were applied: The Geophilus' point-based sensor data were interpolated using the geostatistical method of ordinary kriging, and block kriging was applied for the soil pH and OpticMapper data (R package 'gstat'; Pebesma (2004)). The smoothing procedure block kriging eliminates spatial outliers that show a strong deviation from the surrounding observations. Block kriging produces averaged values within a predefined neighbourhood (block) around the prediction location (Olea 2012). To maintain an appropriate ratio between the prediction of real spatial micro-patterns and the elimination of erroneous sensor measurements, different block sizes were tested, and a suitable block size of 20 × 20 m with low root mean square error (RMSE) values was chosen. This block size allows for inclusion of measurements along one track and measurements from neighbouring tracks.
Due to the sensor high measurement intervals and the consequent high spatial resolution, two criteria (whichever applies first) were established to reduce the number of neighbouring points in the ordinary kriging procedure to considerably reduce the computation time: (i) the maximum distance from the prediction location was set to 100 m and (ii) the maximum number of nearest observations was set to 100. To facilitate automation of the applied processes, this localized kriging approach allows the computation time to be reduced and avoids the complexity of filtering model variograms for local prediction models (Hengl 2009). The final raster data sets created had a spatial resolution of 2 × 2 m for each parameter and were clipped to the boundary of the field.
More advanced geostatistical methods could have been used (e.g., kriging with local variograms, external drift kriging or modelling of spatial anisotropy). However, the geostatistical methods were restricted to simple methods here to keep the focus on the main topic of this research, which is the complete workflow of sensor-based site-specific liming.
Moreover, considering their feasibility under practical conditions, geostatistical modelling efforts should be reasonable regarding the extent to which interpolation errors can be minimized. For experts, it would be easy to integrate more advanced geostatistical methods. However, in addition to being time intensive, more advanced methods introduce some problems (e.g., overfitting). Furthermore, ordinary kriging with a large amount of data is relatively robust according to Webster and Oliver (2007) and Goovaerts (1997). That is, kriging results will not differ substantially regardless of the variogram modelling and kriging approaches.

Reference soil sampling
Reference soil sampling locations were selected based on the proximal soil sensor data (sensor-guided sampling). To relate the sensor data to the target parameters, 30 reference soil samples were taken from the field for soil texture analysis, 15 for pH and 25 for SOM at locations that met the following three criteria ): (i) The targeted samples cover the entire range of the sensor data (feature space): (ii) From the sensor data, high and low values were selected using the 30% and 70% quantiles. This was in order to have the calibration model be based on a wide range of values. (iii) The location is spatially homogeneous: (iv) To avoid the sampling of outliers or erroneous sensor measurements, high and low values should be clustered within a radius of 30 m around the reference sampling point. (v) The samples are well distributed throughout the area of investigation: Conditioned Latin hypercube sampling by Minasny and McBratney (2006) (using R package 'clhs'; Roudier (2011)) was applied to spread the sampling points evenly over the field by maximizing the distance between them. This was done by means of stratified sampling using x and y co-ordinates and the consecutive point ID of the sensor measurements.
Soil samples were taken along the soil sensing trajectory. At each reference sampling point, five subsamples were taken with an auger from 0 to 0.3 m depth within a radius of 0.5 m. The bulked samples were oven-dried at 75 °C and sieved to less than 2 mm in the laboratory. The pH value was measured in 10 g of soil and 25 ml of 0.01 M CaCl 2 solution according to DIN ISO 10390 with a glass electrode after 60 min. The particle size distribution of the < 2 mm fraction was determined according to the German standard for soil science (DIN ISO 11277) by wet sieving and sedimentation after the removal of organic matter with H 2 O 2 and dispersion by 0.2 N Na 4 P 2 O 7 . The soil organic carbon (SOC) was analysed by elementary analysis using the dry combustion method (DIN ISO 10694) after removing the inorganic carbon with hydrochloric acid. To calculate the amount of SOM, the SOC was multiplied by 1.72 (Peverill et al. 1999).

Spatial prediction of soil texture, pH and SOM
To construct relationship models between the sensor data and the lab-analysed soil properties, the interpolated on-the-go sensor data were extracted at the reference sampling locations. Calibrating the interpolated sensor data (particularly the pH data) resulted in better models (lower RMSE, higher R 2 ) than the models developed by calibrating the sensor point data first and interpolating afterwards because pre-processing and interpolation removes some noise from the sensor data.
Since the calibration of the pH sensor data is solely related to the lab-analysed pH values, a univariate linear regression (ULR) model was generated. The predictions of the three soil texture fractions and SOM, on the other hand, were based on the Geophilus (ERa, γ, DEM, SWI) and OpticMapper (Red, IR) data, respectively. Hence, multi-variate linear regression (MLR) models were applied as: where z is the dependent variable at the ith site; X 1 ,X 2 , …, X n are the ancillary data measured at the same site; b 0 , b 1 , b 2 , …, b n are the n + 1 regression coefficients; and is the random error. Before MLR modelling was applied, the interpolated sensor data were checked for their predictive power. If Pearson's correlation coefficient (R) of two variables was found to be larger than 0.65, the variable that correlated best with the target soil property was chosen. Based on the reduced data set of independent variables, a backward stepwise selection (R package 'caret'; Kuhn et al. (2019)) was conducted to find the best set of predictive variables for the MLR model. To assess the accuracy of the MLR models, a k-fold cross-validation was applied one hundred times with k = 3 for SOM and k = 4 for the soil texture prediction. The accuracy of each model was determined using the root mean square error (RMSE) and the coefficient of determination (R 2 ).
Here, clay, silt and sand were considered as fractions summing to 100% or 1 kg kg −1 and having non-negative values (De Gruijter et al. 1997). Hence, when the soil fractions are estimated individually from MLR models, compositional data rules apply to the predicted values (Huang et al. 2014;Muzzamal et al. 2018). To meet these requirements, an additive log-ratio (ALR) transformation was performed (R package 'compositions'; van den Boogaart and Tolosana-Delgado (2008)) following the approaches of Chayes (1960) and Aitchison (1982). In ALR, no fraction is interpreted as isolated from the others. The two advantages of this approach are (i) the removal of closure effects and (ii) the production of suitable data for classical statistical analysis, such as MLR, because the transformed values may be closer to a normal distribution than the untransformed data through perturbation (Odeh et al. 2003).

Determination of variable lime requirements (CaO amounts)
In this study, an empirical lime requirement algorithm was utilized and was adapted to the needs of high-resolution soil data. The conventional VDLUFA approach consists of a look-up table system that allows farmers in Germany to very easily determine the lime requirement (LR) as the amount of CaO that needs to be applied to adjust the soil pH value towards the optimum level and maintain that level until the next fertilization cycle (von Wulffen et al. 2008). This approach defines five pH/lime supply classes for five mineral soil texture classes (Table 1) and for a peat soil class as well as four SOM classes (≤ 4 g kg −1 , 4.1…8 g kg −1 , 8.1…15 g kg −1 , 15.1…30 g kg −1 ) for arable land. The current pH values in classes A and B are further subdivided into small 1/10 pH unit steps.
This rough soil texture and SOM classification system contrasts with the sensitivity and density of the information mapped with mobile on-field sensor systems. Thus, the conventional VDLUFA approach was improved by deriving a continuous or 'stepless' algorithm, i.e., using real values for the three soil properties instead of classified integer values. The adaptation is briefly summarized here. First, a central value was defined for each VDLUFA  (Table 2) were related to both the five clay contents and the four SOM contents as reported above. Third, non-linear regressions were used to calculate the functional relationships that allow the estimation of the respective lime supply level (A-E) for any combination of clay and SOM content. Finally, the lime fertilization recommendation can be calculated depending on the difference between the current and the target pH (lime supply level C, Table 2) as well as the actual clay and SOM content.

Data aggregation and evaluation of the variable lime requirements
Because accurate GNSS receivers and auto-guidance systems are available at reasonable prices, controlled trafficking has gained much popularity and can almost be seen as an integral part of precision agriculture in practice. Consequently, prescription maps for liming should consider the fixed tramlines and working widths used in controlled trafficking. The results were therefore not only shown for the potentially highest resolution of 2 × 2 m but were also aggregated for possible lime spreader working widths of 18 × 18 m and 36 × 36 m for management purposes.
To evaluate the novel VRL approach, the lime amounts from the generalized variable LR maps were compared with possible LRs from a uniform liming strategy, and each management unit (e.g., 18 × 18 m) was determined to be either under-, adequately or overfertilized by the uniform liming approach. Therefore, the pH range for each management unit was computed by subtracting or adding the pH RMSE from each modelled pH value. Afterwards, these pH ranges were used to calculate CaO threshold values for over-and under-fertilization using the stepless algorithm described above. For simplicity, the estimates are based on the error (RMSE) of the derived pH map only, as pH has been determined to be the most important soil property for LR estimation in the investigated soils (Vogel et al. 2020).
The uniform lime demand was determined from the VDLUFA look-up table system using the average pH, SOM and soil texture values for the corresponding SOM and soil texture groups as derived from the sensor-based soil maps of the test field. To account for the coarse classification system of the conventional VDLUFA approach, the uniform estimated LR based on the known soil texture group was additionally compared to estimated LRs based on other potentially selectable soil texture groups.

Geostatistics
The empirical semivariograms and the fitted models for all on-the-go sensor data formed the basis for the interpolation of the sensor point measurements by ordinary kriging (Fig. 4). The selected semivariogram models and the derived variogram parameters sill, nugget and range are summarized in Table 3. The nugget indicates that the sensor data show no or very low spatial micro-variance and random error in their measurements. The spatial correlation structure of the sensor data on the test field can be best characterized by circular (γ, elevation, pH, ECa), exponential (Red, IR) and Gaussian (ERa) models. Cutoffs were set at a distance when a first local maxima is reached or became slightly visible. Due to the exponential character of the fitted semivariogram model for the OpticMapper sensor data, the sill, i.e., the limit of spatial correlation, is reached at rather low ranges of  Table 3). Dashed lines (red) represent the spatial separation distance up to which point pairs are included in the semivariance estimates (Color figure online) 1 3 55 (Infrared) and 60 m (Red), indicating very high spatial variability in that optical soil characteristic. The remaining sensor data showed slightly lower spatial variability with higher ranges of 170 (ERa) to 316 m (pH).

Regionalized sensor data
The interpolated mapping results are shown in Fig. 5. They have a spatial resolution of 2 m and show distinct spatial patterns. The colour scales (and displayed value ranges) of the ERa, ECa and γ data as well as those for the SWI indicate the moisture and/or textural condition at the specific location. For example, for ERa, while values of < 100 Ω m indicate areas with high soil moisture and/or higher clay content, values of > 150 Ω m represent the driest and/or most sandy areas (Fig. 5). Since both ERa and soil γ are strongly related to soil texture, the low-resistivity areas correspond well to the high γ-activity areas, and vice versa. Differences between the patterns in the two maps can be explained by the different soil moisture sensitivities of the two sensors, as shown in the SWI map, with lower values indicating dryer areas and higher values indicating wetter areas. ERa and ECa represent the same content, as they are reciprocal values, and the scales and colours are arranged accordingly to provide similar interpretations. Lower values of ECa (< 3 mS m −1 ) indicate dryer and more sandy areas, and higher ECa values (> 6 mS m −1 ) indicate higher soil moisture and clay contents. The OpticMapper sensor data are characterised by a large cluster of high red and IR values (dimensionless) in the southern and south-western parts of the field. Lower values can be found in the immediate surroundings to the north and to the east as well as in the northernmost part of the field. The spatial patterns of the IR map show more contrasts than those of the red map, whereas the IR/red ratio map shows patterns that are almost identical to the IR patterns. The sensor pH values in field KL60 show four different zones. The northern part of the field is characterized by the highest pH values. To the southeast, there are intermediate pH values and, farther to the south, the pH increases slightly. In the southernmost part of the field, however, the pH values reach their minimum.

Calibration of sensor data
The sensor-based spatial prediction models for pH, soil texture and SOM were calibrated and validated using the lab analysed soil samples selected based on the sensor maps. The descriptive statistics for the reference soil samples are summarized in Table 4. The prediction performance of the univariate linear regression model for pH was very good, with an R 2 of 0.91 and an RMSE of 0.37 (Table 5). The calibrated pH values are lower than the field-measured sensor pH values (Fig. 6). This occurred for the following three reasons: i. The field pH was measured with antimony electrodes instead of with the glass electrodes that are standard in the laboratory. ii. The field pH was measured in tap water, which has a neutral to slightly alkaline pH value, whereas the lab analysis was performed with 0.01 M CaCl 2 solution. Due to the exchange processes of Al 3+ by Ca 2+ at the surface of soil colloids, the pH measured in salt solution is generally lower by 0.6 (± 0.2) pH units. Furthermore, in salt solution, there is no suspension effect to balance the diffusion potential between the pH electrode and the soil solution (Blume et al. 2016). iii. The exposure time between the soil and the solution in which the pH value is measured is a maximum of 20 s (Lund et al. 2005) in the field compared to 60 min during MLR models were used to regionalize the SOM content and soil texture with the sensor data. After testing the proxy variables for independence, the SOM content was predicted using the covariates IR, SWI and ECa. These and the lab-analysed SOC results multiplied by 1.72 were used to calibrate the sensor data. The prediction performance is shown in Fig. 6, showing that the RMSE for SOM was 6.4 g kg −1 with a range of approximately 55 g kg −1 .  After analysing ERa, γ, SWI and elevation for interdependence, only γ and ERa were used as independent variables for predicting the soil texture fractions of sand and clay in the combined MLR and ALR approach. The soil texture prediction results are shown in Figs. 6c-e. The good performance of the models is reflected by, e.g., the prediction of the clay and sand fraction; 87 and 88% of the variability could be explained, and the corresponding RMSE values were 0.046 kg kg −1 and 0.072 kg kg −1 , respectively. Due to the log-ratio transformation of the two predictors, the sand and clay fractions, the prediction of the silt fraction was poor, with an RMSE of 0.039 kg kg −1 .

Generated soil maps
The soil maps of the pH, clay and SOM in KL60 are the basis for the calculation of the site-specific lime requirement of the field following the VDLUFA guidelines for liming in Germany. The descriptive statistics for the predicted soil properties can be found in Table 5.
Regarding the soil texture, derived regression models were applied to the interpolated Geophilus raster data sets, and the soil texture fractions were predicted for the entire field (Fig. 7). Sand is the dominant soil texture fraction, followed by the clay fraction and a comparatively low mean silt content. However, the sand and clay fraction values have relatively large ranges of 0.7 kg kg −1 and 0.44 kg kg −1 , respectively, whereas the silt fraction values have a comparatively low range of 0.26 kg kg −1 . Larger portions of sand were found in the more elevated areas to the south, the south-western part of the field and near the eastern and north-eastern borders. A more linear sandy zone stretches out from the southeast to the northwest in the centre of the field. This might have been formed by streams as part of the (post-) Palaeozoic glacial landscapes, which are well known for their high soil and sediment variability (Krbetschek et al. 2008). Glacial, periglacial and interglacial processes created a mosaic of landforms and unconsolidated sediments that vary over short distances. Clayey areas dominate the lower elevated flanks along the sandy areas from the southeast to the northwest, indicating lower water drainage. The silt fraction in this field shows less pronounced variation than the sand and clay fractions. The classified soil texture map (Fig. 8a) shows that the classes, according to the German soil classification system KA5, cover a relatively wide range from pure sand (Ss) to slightly loamy sand (St2), highly loamy sand (Sl4), medium clayey sand (St3), highly sandy loam (Ls4) and clayey sandy loam (Lts). The distribution of the classes corresponds well to the findings of the clay, silt and sand distribution and provides some clarification. An area of approximately 4 ha is covered by pure sand only in the southern part of the field. Slightly loamy sand (St2) covers a total area of approximately 7 ha. The classes Sl4 and Ls4 are only visible in tiny patches of less than 5 ha combined. They occur as transition areas to clayey sandy loam (Lts, which covers three larger areas of approximately 5.5 ha in the centre, the west and the north of the field) and very slightly loamy sand (St2, which covers approximately 7.7 ha in total). The KA5 soil classification was chosen to avoid the insufficient spatial resolution of the VDLUFA soil texture classification system. Fig. 8 Geophilus mapped soil texture classes (derived from the German KA5 soil texture classification (Eckelmann et al. 2005)) Ss pure sand, St2 slightly loamy sand, St3 medium clayey sand, Sl4 highly sandy loam, Ls4 highly sandy loam, Lts clayey sandy loam (a), the MSP-mapped pH values (b) and SOM content (c) and the lime supply level at 2 × 2 m resolution In other studies, Boenecke et al. (2018) and Meyer et al. (2019) used data from the Geophilus system to successfully generate predictive soil texture maps of the clay, silt and sand fractions of the topsoil for practical purposes. Meyer et al. (2019) achieved the best prediction results by deriving the soil texture of the topsoil using the gamma mapping results and by calculating the dimensionless relationship between the gamma and electrical resistivity mapping results.
The calibrated pH values in field KL60 ranged between 4.2 and 7.1 (Fig. 8b) with a median of 6.4. Since the pH calibration is based on a univariate linear regression model, the spatial patterns of the calibrated data were identical to the sensor data and indicated four different zones of soil acidity. The error of the pH measurements was 0.07 pH units larger than that determined by Adamchuk and Lund (2008). It is striking that the lowest pH values were found in the sandy regions of the field, which also showed the lowest amounts of SOM. In contrast, pH values were highest in the loamy parts of the field that had higher SOM contents.
The field KL60 is characterized by low SOM content, showing a mean of 30 g kg −1 , a range of 55.5 g kg −1 and a standard deviation of 10.5 g kg −1 . The spatial patterns of the SOM map show many similarities to those of the soil texture map (Fig. 8c). The slightly elevated sandy hilltops in the southern and central parts of the field are characterized by lower amounts of SOM. In contrast, higher amounts can be found in the lower-lying areas, coinciding with a loamy soil texture and higher pH values.
In a study by Kuang et al. (2015), artificial neuronal network (ANN) and partial least square regression (PLSR) were used for the calibration of a visible and near infrared (vis-NIR) sensor data to generate high-resolution maps of pH, SOC and clay. Using a high number of soil reference samples for their two test fields (n = 132 and n = 80), calibration with ANN outperformed the PLSR method. For example, the root mean square error (RMSE) for the ANN calibrated sensor data was 12.5 g kg −1 for SOC, 0.12 for pH and 0.0096 kg kg −1 for the clay content (PLSR: 14.8 g kg −1 for SOC, 0.13 for pH, 0.0105 kg kg −1 for clay content).
Regarding the lime supply status of KL60, approximately 20% of the field requires recovery or build-up liming (levels A and B), and 16.5% requires no liming at all (levels D and E). Nearly two-thirds of the field is within the optimal lime supply range (level C) and only requires maintenance liming (Fig. 8d). By decreasing the resolution to management conform units, the areal percentages did not change considerably (Table 6).

Determined lime requirements and data aggregation
The high-resolution LR map (expressed in t CaO ha −1 ) was used to generate a prescription map for liming whose spatial resolution was adapted to the working width of a lime spreader (Fig. 9). For that, the CaO data were resampled to an 18 × 18 m raster grid and Table 6 Areal percentages of the VDLUFA lime supply levels (A, B, C, D, E, description in  aligned to the appropriate management direction of the field (Fig. 9b). For spreaders with larger working widths, e.g., 36 × 36, the creation of maps with larger raster widths leads to information losses and increasingly over-or under-limed portions of the field (Fig. 9c). The map shows that the total LR of field KL60 is rather low. Nevertheless, relatively high spatial variability exists that can only be well explained by the combination of all three soil maps: pH, soil texture and SOM (Table 7). The 2 × 2 m resolution map shows that CaO values had a range of more than 7 t ha −1 , with some areas requiring no lime at all and some areas showing very high demand. Low lime requirements were identified in the northern and central parts of the field where pH values are high, SOM content is low and the soil texture is dominated by sand. In contrast, the higher lime requirements in the northern central areas are particularly caused by the loamy soil textures that increase the target pH value according to the VDLUFA guidelines for liming. As mentioned above, this is because soils with a higher clay content require more lime to stabilize their soil structure. Furthermore, clayey soils can be prone to aluminium toxicity, which can be counteracted by liming (Schilling 2000;Blume et al. 2016). The highest lime requirements on field KL60 can be seen in the south and southeast. These areas coincide with the lowest pH values and sandy soil textures. The aggregation of the 2 × 2 m resolution map to management units of 18 × 18 m and 36 × 36 m revealed that the maximum LRs decreased by approximately 6% and 13%, respectively. The cumulative LR amounts increased by 11% and 25%, respectively.

Evaluation and comparison of the variable-rate lime requirements
The variable liming results were compared with a uniform liming approach in which the lime demands were determined following the conventional VDLUFA guidelines. Assuming that the sensor based LR map most closely reflects the actual LR conditions in the field, uniform liming would result in certain areas of the field being under-, adequately or over-fertilized (Table 8, Fig. 10). According to the mean values for the clay, silt and sand fractions in the test field (Table 5), VDLUFA soil texture group 2 (Table 1) would be most suitable for assessing the LR in the uniform fertilization strategy. The errors and uncertainties of the conventional approach can hardly be quantified. While the pH and SOM content can be more easily assessed by laboratory analysis, the texture values are often measured by quick on-field methods in conventional farming practices. However, this method is highly prone to error (Stocker and Walthert 2013) and may lead to the potential selection of adjacent soil groups (e.g., in this study, soil group 1 or 3). Moreover, these errors and variances are usually not considered in practice. Within this study, the uncertainties of the conventional approach were therefore expressed by comparing the lime demands of the potentially selectable soil groups as per the VDL-UFA soil classification system (Table 1). It also needs to be emphasized that a soil mapping process is also not free from error. A digital soil map is the outcome of several consecutive steps that are associated with uncertainties. These steps include soil sampling, laboratory analysis and final digital soil mapping (comprising sensor data interpolation and parameter calibration). Uncertainties can be caused by several factors, such as the measurement methods and tools as well as the natural variability of the soil. For example, even high-resolution surveys suffer from the fact that a field cannot be measured at each individual point. Such uncertainties are discussed widely in the literature (Brendan et al. 2017;Piiki and Söderström 2019). Due to the complexity of error propagation, only a few studies have tried to compare the impacts of sources of errors to optimize the soil mapping process (Mueller et al. 2004;Gebbers and De Bruin 2010). Gebbers and De Bruin (2010) have shown how the relevance of factor uncertainties (e.g., sampling design and interval, positioning error, regionalisation method) can be quantified by global sensitivity analyses of a stochastic simulation model of the soil mapping process. They found that uncertainties due to the sampling density, the spatial variation of soil properties and the prediction method had the greatest influence on the results. Compared to the errors of these factors, the error of the soil chemical analysis had little impact when it was increased from 0 to 20%. Hence, the accuracy of the entire soil mapping process can most efficiently be improved by increasing the sampling density (e.g., by using mobile sensor platforms), while improving the precision of chemical and physical laboratory analyses has a smaller effect.
Regarding the estimates of the uniform LR map based on VDLUFA soil group 2 (SG2), 1.1 t CaO ha −1 should be applied according to the VDLUFA look-up scheme. The total CaO demand based on this soil group is only approximately 1 t different from the total lime demand of the variable LR map (Table 8). While this difference is relatively low, a uniform LR determination using soil group 1 or 3 would lead to nearly half the CaO demand or an almost three times higher CaO demand, respectively, than the total LR determined by the variable approach. While approximately one-third of the field would be adequately fertilized using SG2 for lime demand determination, the adequately fertilized area would increase to approximately half of the field with SG1 and decrease to less than 10% with SG3 ( Fig. 10). In contrast, nearly two-thirds of the field would be over-fertilized by approximately 12 t CaO with SG2, and slightly more than 10% would be under-fertilized. With SG1, in comparison, only one-third of the field would be over-fertilized, and 11% would be under-fertilized. Interestingly, by choosing SG3, the over-fertilized areas of the field would increase to more than 90% of the field, and a total of approximately 55 t too much CaO would be applied to the field. Making up less than 1% of the area, the under-fertilized areas might be neglected. Although SG4 and SG5 would most likely not be chosen in this example, the over-fertilized areas would increase by merely 5% and 7%, respectively, in comparison to those under SG3. However, the determined CaO amounts would therefore double or triple.
Only a few studies have compared variable-rate LR approaches with uniform LR approaches based on mean values. In North America, Borgelt et al. (1994), for example, compared variable lime rates with a uniform mean liming approach as well as with LRs estimated by a soil-buffer and a rule-based method. The latter was based on the parameters crop type, soil pH and soil texture. They produced variable-rate liming maps using geostatistical analysis and soil samples taken from a modified soil sampling design. Although the 8.8 ha test field in their study was less heterogeneous than that in this study, showing only two soil texture classes compared to KL60 (5 main groups, 6 subclasses), the mean liming 1 3 rates would have resulted in over-fertilization of 9 to 12% and under-fertilization of 37 to 41% of the field. Overall, the uniform (mean-based) liming approach would have resulted in an 8 to 28% lower total lime application. In another study by Zaman et al. (2003) from the UK, 35% of the tested field required more lime under the variable approach than the average amount, 56% less and 9% of the field was adequately limed. The field, however, had low ranges of sand (mean = 0.5 kg kg −1 , range = 0.1 kg kg −1 ), silt (0.24 kg kg −1 , 0.11 kg kg −1 ) and clay (0.26 kg kg −1 , 0.11 kg kg −1 ) contents, and the LRs were estimated based on pH only. In two fields in North America, Bianchini and Mallarino (2002) found that 56-61% less lime needed to be applied under a variable lime rate approach based on a very dense sampling grid than under uniform application. Determination of variable LR in all these studies was, however, based on regular sampling grids for soil texture and pH and none of these studies used spatial interpolation methods or high-resolution soil maps based on proximal soil sensing. In a study by Kuang et al. (2014), conducted on an 18 ha field in Denmark, variable liming rates were derived from high resolution mapping with a Vis-NIR spectrometer system and a recommendation algorithm from the Danish Centre for Food and Agriculture (DCA). They found that the VRL consumed the largest amount of lime (37 t) for the entire field, while the uniform treatment required just 25.3 t for entire field. However, for variable lime management the yield of spring barley of 7.57 t ha −1 was slightly but significantly better than the yield of 7.51 t ha −1 under uniform lime management.

Outlook
The success of lime applications based on this approach is currently being studied in field trials at different study sites and will be verified within the project period by repeated sensor campaigns with the multi-sensor platform and yield measurements. Moreover, the fusion of the sensor data within the project will be tested to enhance the predictability of the required soil parameters for liming. The financial aspects of this approach have not been addressed so far. It is evident that the costs to calibrate the sensors through soil sampling and soil analysis should be as low as possible. The number of reference samples taken in this study was relatively high. Part of the ongoing research is to reduce the number of reference samples to a maximum of 5 to 10 per field for pH, soil texture and SOC together while keeping prediction accuracy at a sufficient level. On the other hand, the economic and environmental value of precise liming must be highlighted. This value will only be perceived in the long term, and the adoption of precision liming will likely be supported by the relevant authorities. Farmers, advisors and service providers need training and accessible software tools to obtain the full benefits of existing soil mapping systems.

Conclusions
The present study presents a developed procedure that allows the easy and semi-automated generation of topsoil pH, texture and soil organic matter maps based on proximal soil sensing. This can be used for soil acidity management practices that respect the natural soil variability at a high level of detail and improve the currently available best practices as described above. Moreover, this study provides guidelines for implementation in practice 1 3 and, for scientists and advisors, information for the comparison and further development of approaches to variable-rate liming.
It was shown that high-resolution soil maps of pH, soil texture and soil organic matter could be generated through sensor-based digital soil mapping using two multi-sensor platforms and semi-automated (geo-) statistical methods. These soil maps exhibited smallscale spatial patterns and spatial interrelations between the target variables. More elevated parts of the field are characterized by a sandy soil texture, low amounts of SOM and low pH values. Zones at low elevations, which most likely developed from fluvial processes, are characterized by loamy soil textures, higher SOM contents and higher pH values.
Based on the high-resolution soil maps, a lime requirement map at 2 × 2 m spatial resolution was calculated following an adapted approach to the conventional VDLUFA guidelines for liming in Germany. However, to generate a lime prescription map that can be processed by a currently available lime spreader as used in the present study, the spatial resolution needed to be changed to 18 × 18 m. Given the high resolution of input data that proximal soil sensors can provide, the lack of precision in the currently available lime applicators is probably a bottleneck for the improvement of soil acidity management.
The combined average soil map error was used as a threshold value for identifying over, adequately and under-limed areas, and the LRs of the precision and the uniform liming approaches were compared. The results showed that 59% of the field would be over-fertilized by approximately 12 t of lime, 24% would receive approximately 10 t too little lime and merely 17% would be adequately limed with the uniform approach.