Search

Digital soil mapping (DSM) approaches provide soil information by utilising the relationship between soil properties and environmental variables. Calibration of DSM models requires measurements that may often have substantial measurement errors which propagate to the DSM outputs and need to be accounted for. This study applied a geostatistical-based DSM approach that incorporates measurement error variances in the covariance structure of the spatial model, weights measurements in accordance with their measurement accuracies and assesses the effects of measurement errors on the accuracies of DSM outputs. The method was applied in the Western Cameroon, where soil samples from 480 locations were collected and analysed for pH, clay and soil organic carbon (SOC) using conventional and mid-infrared spectroscopy methods. Variogram parameters and regression coefficients were estimated using residual maximum likelihood under two scenarios: with and without taking measurement errors into account. Performance of the spatial models in the two scenarios was compared using validation metrics obtained with three types of cross-validation. Acknowledging measurement errors impacted the regression coefficients and influenced the variogram parameters by reducing the nugget and sill variance for the three soil properties. Validation metrics including mean error, root mean square error and model efficiency coefficient were quite similar in both scenarios, but the prediction uncertainties were more realistically quantified by the models that account for measurement errors, as indicated by accuracy plots. There were relatively small absolute differences in predicted values of soil properties of up to 0.1 for pH, 1.6% for clay and 2 g/kg for SOC between the two scenarios. We emphasised the need of incorporating measurement errors in DSM approaches to improve uncertainty quantification, particularly when applying spectroscopy for estimating soil properties. Further development of the approach is the extension to non-linear machine learning regression methods. Highlights: Errors in soil measurements are usually not accounted for and may affect DSM results. Measurement error variances were incorporated in the geostatistical models of three soil properties. Quantifying measurement errors in DSM allows to weigh measurements in accordance with their accuracy. Accounting for measurement errors in DSM better assesses prediction accuracy.

There is a growing demand for high quality soil data to model soil processes and map soil properties. However, wet chemistry measurements on soil properties are subjected to many error sources, such as the observer, the instrument and lack of standardised methodologies. Consequently, soil data are imperfect and uncertain because of these error sources. Uncertainties in measurements of fundamental soil properties can propagate through, e.g., pedotransfer functions, spectroscopic models and digital soil mapping algorithms. Therefore, it is important to provide detailed uncertainty information about soil measurements to potential data users. In practice, uncertainty estimates are rarely specified by providers of analytical soil data.In this research, we aimed to quantify uncertainties in synthetic and real-world pH (1:1 soil-water suspension) and Total Organic Carbon (TOC) measurements. We assumed that uncertainty can be represented by a normal distribution. A linear mixed-effects model was applied to estimate the parameters of the normal distribution, i.e., mean and standard deviation, of both synthetic and real-world datasets. The model included ‘sample ID’ as a fixed effect, and ‘batch’ and ‘laboratory’ as random effects. The use of synthetic datasets allowed us to investigate how well the model parameters could be estimated given a specific experimental measurement design, whereas the real-world case served to explore if the parameter estimates were still accurate for such unbalanced datasets.For a balanced dataset (n=20, n=100, n=200 and n=500), using synthetic pH data for three hypothetical laboratories (two batches per laboratory), the mean estimated standard deviations (σ) of the random effects were σbatch=0.10, σlaboratory=0.24 and σresidual=0.2. These estimates were in agreement with the σ for the respective random effects used to generate the synthetic dataset, meaning that the model could accurately estimate the model parameters. Subsequently, changes were made to the experimental measurement design by randomly removing 20%, 50% and 80% of the data, resulting in unbalanced datasets. In general, the interquartile range (IQR) of σ for each random effect increased with a larger percentage of removed data. However, the increase in IQR was larger for n=20 compared to, e.g., n=200. When comparing 0% and 80% randomly removed data, the IQR for the batch effect increased with 60.3%. Conversely, for n=200 an increase of only 23.5% was observed.Subsequently, the same model was fitted on real-world pH and TOC data, provided by the Wageningen Evaluating Programs for Analytical Laboratories (WEPAL). The unbalanced dataset structure was first reconstructed and filled with synthetically generated data, based on sample means and standard deviations derived from the measured data. The model was fitted on both datasets. For measured pH, the model yielded σbatch=0.27, σlaboratory=0.17 and σresidual=0.10. The IQRs of the estimated σ from synthetic WEPAL data were 0.04 (batch), 0.06 (laboratory) and 0.02 (residual). The model fitted on the measured TOC data estimated σbatch=5.3%, σlaboratory=2.8% and σresidual=2.1%. For the synthetic WEPAL data, IQRs of 1.3% (batch), 1.4% (laboratory) and 0.4% (residual) were determined for the estimated σ. These findings suggest that despite having a highly unbalanced dataset, realistic model parameter estimates can still be obtained.

Spatial soil applications frequently involve binomial variables. If relevant environmental covariates are available, using a Bayesian generalized linear model (BGLM) might be a solution for mapping such discrete soil properties. The geostatistical extension, a Bayesian generalized linear geostatistical model (BGLGM), adds spatial dependence and is thus potentially better equipped. The objective of this work was to evaluate whether it pays off to extend from a BGLM to a BGLGM for mapping binary soil properties, evaluated in terms of prediction accuracy and modelling complexity. As motivating example, we mapped the presence/absence of the Pleistocene sand layer within 120 cm from the land surface in the Dutch province of Flevoland, using the BGLGM implementation in the R-package geoRglm. We found that BGLGM yields considerably better statistical validation metrics compared to a BGLM, especially with – as in our case – a large (n = 1,000) observation sample and few relevant covariates available. Also, the inferred posterior BGLGM parameters enable the quantification of spatial relationships. However, calibrating and applying a BGLGM is quite demanding with respect to the minimal required sample size, tuning the algorithm, and computational costs. We replaced manual tuning by an automated tuning algorithm (which eases implementing applications) and found a sample composition that delivers meaningful results within 50 h calculation time. With the gained insights and shared scripts spatial soil practitioners and researchers can – for their specific cases – evaluate if using BGLGM is feasible and if the extra gain is worth the extra effort.

Soils provide a variety of goods and services and are a key natural resource, non-renewable on the time scale of a human life-span, to realise several UN-Sustainable-Development–Goals, including zero hunger (SDG2), climate action (SDG13) and life on land (SDG15). Consistent soil information is required to underpin a large range of global assessments, such as soil and land degradation, sustainable land management, and environmental conservation. In this study, we present an application addressing the modelling of soil functions at global scale, using erosivity risk and soil carbon sequestration potential as examples. We used SoilGridsv2.0, as set of soil property maps at 250m resolution, derived from a large set of standardized soil profile observations (WoSIS), an updated set of environmental covariates, and improved machine learning models. It provides global assessments of prediction uncertainty, quantified with a 90% prediction interval, and considers an evaluation procedure that provides more realistic metrics of map accuracy. We used simplified models to derive soil functions (e.g. erosivity or carbon sequestration potential) from basic soil properties, meaningful for different pedo-climatic regions. We provide an indication of areas of low/high risk of soil degradation to support sustainable soil management planning. The uncertainty limits of the input soil layers were used to provide a preliminary assessment of the uncertainty of the derived layers. The present modelling framework, using soil properties maps to derive soil functions, offers great flexibility and may be applied to a diverse set of models to generate soil information products tailored to specific applications. We highlight some of the challenges of assessing soil functions at global scale.

Since the establishment of Digital Soil Mapping (DSM) as a research field, the main focus has been on implementing new methods to improve the predictive performance of soil maps. However, considerably less effort has been invested in investigating the best way to communicate the quality of soil mapping products with users. This is essential for soil maps to be adopted by a broader community, future research guidance and most importantly, to ensure that they are used correctly. We introduce a high-resolution 3D soil modelling and mapping platform for the Netherlands (BIS-3D) using a quantile regression forest (QRF) for spatial interpolation approach that includes an assessment of the map quality using GlobalSoilMap (GSM) accuracy thresholds. Our objectives are twofold: a) providing accurate and high-resolution (25m) soil pH, soil organic carbon, and soil texture (clay, silt, and sand) maps over 3D space including prediction uncertainty; and b) providing an intuitive way to communicate accuracy of soil maps for users by means of accuracy thresholds. In this work, the first outputs of the modelling and mapping platform BIS-3D are being presented.QRF models were trained and validated, yielding average predictions for each target location and depth as well as the 90% prediction interval. Predicted soil maps were evaluated using an independent validation data set based on a stratified random sampling design covering the entire Netherlands (1151 locations). Furthermore, at every validation location, predictions were assessed as A, AA or AAA quality using the GSM specifications.First results for soil pH (KCl) using 15887 soil observations between depths 0-2 m and 180 covariates reveal a mean square error skill score (SSmse) = 0.88, RMSE = 0.49 and bias = 0.01 for out of bag predictions. Model evaluation using the independent validation set resulted in SSmse = 0.66, RMSE = 0.81 and bias = 0.12 across all depths. Prediction accuracy was highest for depths between 0-15 cm (SSmse = 0.66, RMSE = 0.76) and 60-100 cm (SSmse = 0.69, RMSE = 0.78) and lowest for 100-200 cm (SSmse = 0.61, RMSE = 0.86). The soil measurement (observation) was within the 90% prediction interval of model predictions in 83% of the cases, indicating that QRF is slightly over-optimistic in quantifying the prediction uncertainty. 61% of predictions that were independently validated over all depths were within the highest GSM accuracy threshold (AAA = +/- 0.5 pH), 23% were AA (+/- 1.0 pH), 9% were A (+/- 1.5 pH) and the remaining 7% were below A. A categorical physical geography map was the most important covariate, although other covariates associated with relief, geomorphology, land use and temperature were also effective. However, such variable importance measurements are merely indications and should be handled with care. The BIS-3D can easily be extended for predicting additional soil properties and it may provide a basis for decision makers to easily assess to what extent and in which areas soil maps can be used for their applications.

Digital soil mapping (DSM) may be defined as the use of a statistical model to quantify the relationship between a certain observed soil property at various geographic locations, and a collection of environmental covariates, and then using this relationship to predict the soil property at locations where the property was not measured. It is also important to quantify the uncertainty with regards to prediction of these soil maps. An important source of uncertainty in DSM is measurement error which is considered as the difference between a measured and true value of a soil property.The use of machine learning (ML) models such as random forests (RF) has become a popular trend in DSM. This is because ML models tend to be capable of accommodating highly non-linear relationships between the soil property and covariates. However, it is not clear how to incorporate measurement error into ML models. In this presentation we will discuss how to incorporate measurement error into some popular ML models, starting with incorporating weights into the objective function of ML models that implicitly assume a Gaussian error. We will discuss the effect that these modifications have on prediction accuracy, with reference to simulation studies.

Spatially resolved estimates of change in soil organic carbon (SOC) stocks are necessary for supporting national and international policies aimed at achieving land degradation neutrality and climate change mitigation. In this work we report on the development, implementation and application of a data-driven, statistical method for mapping SOC stocks in space and time, using Argentina as a pilot. We used quantile regression forest machine learning to predict annual SOC stock at 0–30 cm depth at 250 m resolution for Argentina between 1982 and 2017. The model was calibrated using over 5,000 SOC stock values from the 36-year time period and 35 environmental covariates. We preprocessed normalized difference vegetation index (NDVI) dynamic covariates using a temporal low-pass filter to allow the SOC stock for a given year to depend on the NDVI of the current as well as preceding years. Predictions had modest temporal variation, with an average decrease for the entire country from 2.55 to 2.48 kg C m−2 over the 36-year period (equivalent to a decline of 211 Gg C, 3.0% of the total 0–30 cm SOC stock in Argentina). The Pampa region had a larger estimated SOC stock decrease from 4.62 to 4.34 kg C m−2 (5.9%) during the same period. For the 2001–2015 period, predicted temporal variation was seven-fold larger than that obtained using the Tier 1 approach of the Intergovernmental Panel on Climate Change and United Nations Convention to Combat Desertification. Prediction uncertainties turned out to be substantial, mainly due to the limited number and poor spatial and static, whereas SOC is dynamic and SOC dynamics are of particular interest to carbon sequestration and land degradation studies. Thus, there is a clear need to extend spatial SOC mapping to space–time SOC mapping. temporal distribution of the calibration data, and the limited explanatory power of the covariates. Cross-validation confirmed that SOC stock prediction accuracy was limited, with a mean error of 0.03 kg C m−2 and a root mean squared error of 2.04 kg C m−2. In spite of the large uncertainties, this work showed that machine learning methods can be used for space–time SOC mapping and may yield valuable information to land managers and policymakers, provided that SOC observation density in space and time is sufficiently large. Fil: Heuvelink, Gerard B.M. ISRIC - World soil information; Holanda. Wageningen University. Soil Geography and Landscape Group; Holanda Fil: Angelici, Marcos E. Instituto Nacional de Tecnología Agropecuaria (INTA). Instituto de Suelos; Argentina Fil: Poggio, Laura ISRIC - World soil information, Wageningen; Holanda Fil: Bai, Zhanguo ISRIC - World soil information, Wageningen, The Netherlands Fil: Batjes, Niels H. ISRIC - World soil information, Wageningen, The Netherlands Fil: an den Bosch, Rik ISRIC - World soil information, Wageningen, The Netherlands Fil: Bossio, Deborah The Nature Conservancy; Estados Unidos Fil: Estella, Sergio Vizzuality; España Fil: Lehmann, Jhoannes. Cornell University. Soil and Crop Sciences; Estados Unidos Fil: Olmedo, Guillermo F. Instituto Nacional de Tecnología Agropecuaria (INTA). Estación Experimental Agropecuaria Mendoza; Argentina Fil: Sandermann, Jonathan. Woods Hole Research Center; Estados Unidos

Identifying poverty determinants in a region is crucial for taking effective poverty reduction measures. This paper utilizes two variable importance analysis methods to identify the relative importance of different geographic factors to explain the spatial distribution of poverty: the Lindeman, Merenda, and Gold (LMG) method used in multiple linear regression (MLR) and variable importance used in random forest (RF) machine learning. A case study was conducted in Yunyang, a poverty-stricken county in China, to evaluate the performances of the two methods for identifying village-level poverty determinants. The results indicated that: (1) MLR and RF had similar explanation accuracy; (2) LMG and RF were consistent in the three main determinants of poverty; (3) LMG better identified the importance of variables that were highly related to poverty but correlated with other variables, while RF better identified the non-linear relationships between poverty and explanatory variables; (4) accessibility metrics are the most important variables influencing poverty in Yunyang and have a linear relationship with poverty.

Nematodes are indicators of soil quality and soil health. Knowledge of the relationships between nematode-based soil quality indices and environmental properties is beneficial for assessing environmental threats on soil biota. This study evaluated the spatial distribution of nematode-based soil quality indices in a 23-ha heavy metal-polluted nature reserve using geostatistical methods. We expected that a selection of abiotic soil properties (pH and moisture, clay, organic matter, cadmium (Cd), and zinc (Zn) contents) could explain a significant portion of the spatial variation of the indices and that regression kriging could more accurately model their spatial distribution than ordinary kriging. A stratified simple random sampling scheme was used to select 80 locations where soil samples were taken to extract nematodes and derive the indices. The area had a distinct gradient in soil properties with Cd and Zn content ranging from 0.07 to 68.9 and 5.3 to 1 329 mg kg−1, respectively. Linear regression models were fitted to describe the relationships between the indices and soil properties. By also modelling the spatial correlation structure of regression residuals using spherical semivariograms, regression kriging was used to produce maps of the indices. The regression models explained between 21% and 44% of the total original variance in the indices. Soil pH was a significant explanatory variable in almost all cases, while heavy metal conent had a remarkably low effect. In some cases, the regression residuals had spatial structure. Independent validation indicated that in all cases, regression kriging performed slightly better because of having lower values of the root mean square prediction error and a mean prediction error closer to zero than ordinary kriging. This study showed the importance of soil properties in explaining the spatial distribution of biological soil quality indices in ecological risk assessment.

International audience; Pedometrics, the application of mathematical and statistical methods to the study of the distribution and genesis of soils, has broadened its scope over the past two decades. The primary focus of pedometricians has traditionally been on spatial and spatio-temporal soil inventories with numerical soil classification, geostatistical modelling of spatial variation and mapping. The rapid development of remote and proximal soil sensing as well as data-driven statistical modelling techniques have had a major impact on pedometrics over the past decades. During this time, a general demand for quantitative digital soil information for environmental modelling and management has compelled pedometricians to address other soil-related questions from a quantitative point of view: soil genesis and utility and quality of soil. While scientific progress is largely an autonomous process that is difficult to steer, research efforts could benefit from an agenda with pressing pedometric research topics. This paper defines and discusses ten recent or longstanding pedometrics challenges, with the attempt to identify knowledge gaps and suggest new concepts and methods to overcome them. The ten challenges were selected through a collaborative effort and may serve as a guidance for future pedometrics research and to foster collaboration among soil scientists. The challenges discussed in this paper are also indicators of the current understanding and state of knowledge from which progress can be measured in the future.

For decades scientists have produced maps of biological, ecological and environmental variables. These studies commonly evaluate the map accuracy through cross-validation with the data used for calibrating the underlying mapping model. Recent studies, however, have argued that cross-validation statistics of most mapping studies are optimistically biased. They attribute these overoptimistic results to a supposed serious methodological flaw in standard cross-validation methods, namely that these methods ignore spatial autocorrelation in the data. They argue that spatial cross-validation should be used instead, and contend that standard cross-validation methods are inherently invalid in a geospatial context because of the autocorrelation present in most spatial data. Here we argue that these studies propagate a widespread misconception of statistical validation of maps. We explain that unbiased estimates of map accuracy indices can be obtained by probability sampling and design-based inference and illustrate this with a numerical experiment on large-scale above-ground biomass mapping. In our experiment, standard cross-validation (i.e., ignoring autocorrelation) led to smaller bias than spatial cross-validation. Standard cross-validation was deficient in case of a strongly clustered dataset that had large differences in sampling density, but less so than spatial cross-validation. We conclude that spatial cross-validation methods have no theoretical underpinning and should not be used for assessing map accuracy, while standard cross-validation is deficient in case of clustered data. Model-free, design-unbiased and valid accuracy assessment is achieved with probability sampling and design-based inference. It is valid without the need to explicitly incorporate or adjust for spatial autocorrelation and perfectly suited for the validation of large scale biological, ecological and environmental maps.

Many national and international initiatives rely on spatially explicit information on soil organic carbon (SOC) stock change at multiple scales to support policies aiming at land degradation neutrality and climate change mitigation. In this study, we used regression cokriging with random forest and spatial stochastic cosimulation to predict the SOC stock change between two years (i.e. 1992 and 2010) in Hungary at multiple aggregation levels (i.e. point support, 1 × 1 km, 10 × 10 km square blocks, Hungarian counties and entire Hungary). We also quantified the uncertainty associated with these predictions in order to identify and delimit areas with statistically significant SOC stock change. Our study highlighted that prediction of spatial totals and averages with quantified uncertainty requires a geostatistical approach and cannot be solved by machine learning alone, because it does not account for spatial correlation in prediction errors. The total topsoil SOC stock for Hungary was predicted to increase between 1992 and 2010 with 14.9 Tg, with lower and upper limits of a 90% prediction interval equal to 11.2 Tg and 18.2 Tg, respectively. Results also showed that both the predictions and uncertainties of the average SOC stock change were smaller for larger spatial supports, while spatial aggregation also made it easier to obtain statistically significant SOC stock changes. The latter is important for carbon accounting studies that need to prove in Measurement, Reporting and Verification protocols that observed SOC stock changes are not only practically but also statistically significant.

With increasing discrepancies between population growth and food production in China, the monitoring of crop yield is essential to support food security policies. However, current studies about spatio-temporal variation of yield mainly focus on the influence of climatic factors on grain crops, and do not explore the contributions of agricultural, environmental and economic factors on crop yield in China. In this study, a large yield dataset, covering 31 provinces and a 38-year period from 1978 to 2015, and related explanatory variables were collected for analyzing the spatio-temporal variation of different yield aggregations using stepwise multiple linear regression. At the national scale, the average aggregate yield increased from 3.04 Mg ha−1 in 1978 to 10.04 Mg ha−1 in 2015. Overall, the average aggregate yield increased in all provinces but the average annual growth rates varied: it was smaller than 2.5 % in Heilongjiang, Guizhou, Beijing, Qinghai and Jilin, more than 4.0 % in Hainan, Guangxi, Ningxia, Hebei and Shaanxi, and between 2.5 % and 4.0 % in other provinces. The spatial patterns of the average yield from 1978 to 2015 were different for different crop aggregations. Most of the regression models explained more than 60 % of the yield variance, except for rice, potato and cotton models. Agricultural management practices, soil and economic covariates were important explanatory variables in all models. Topography and climatic covariates were also important for some of the crop models. The regression model of the aggregate yield for all crops explained 95 % of the yield variance, which was mainly explained by planting area index of vegetables (20 %), followed by farmer income (14 %), planting area index of other crops (orchards 11 %, melons 8 %, sugar 6 %, cereals 6 %), and density of agricultural diesel engines (5 %). Although the regression residual of the aggregate yield model was zero on average, the trends were different in different provinces: most provinces demonstrated a small negative or positive residual; the yield was substantially lower (< -0.20 Mg ha−1 y−1) than predicted by the regression model in three provinces in central China (Hebei, Shaanxi and Anhui) and substantially higher (> 0.20 Mg ha−1 y−1) in four provinces (Shanxi, Shandong, Sichuan and Guangdong). These systematic over- and underpredictions may be caused by other factors, such as plagues, pests, natural hazards, market structures (such as competition for labor or impediments to market access) and farmer's management skills. With the increasing population and limited agricultural land resources, enhancing economic growth might be an adequate solution to meet the growing demand for food. It can also promote agricultural efficiency in China, certainly when combined with better management practices, crop composition, breeding and planting technologies.

Uncertainty is often ignored in urban water systems modelling. Commercial software used in engineering practice often ignores uncertainties of input variables and their propagation because of a lack of user-friendly implementations. This can have serious consequences, such as the wrong dimensioning of urban drainage systems (UDS) and the inaccurate estimation of pollution released to the environment. This paper introduces an uncertainty analysis framework in urban drainage modelling and applies it to a case study in the Haute-Sûre catchment in Luxembourg. The framework makes use of the EmiStatR model which simulates the volume and substance flows in UDS using simplified representations of the drainage system and processes. A Monte Carlo uncertainty propagation analysis showed that uncertainties in chemical oxygen demand (COD) and ammonium (NH4) loads and concentrations can be large and have a high temporal variability. Further, a stochastic sensitivity analysis that assesses the uncertainty contributions of input variables to the model output response showed that precipitation has the largest contribution to output uncertainty related with water quantity variables, such as volume in the chamber, overflow volume and flow. Regarding the water quality variables, the input variable related to COD in the wastewater has an important contribution to the uncertainty for COD load (66 %) and COD concentration (62 %). Similarly, the input variable related to NH4 in the wastewater plays an important role in the contribution of total uncertainty for NH4 load (34 %) and NH4 concentration (35 %). The Monte Carlo simulation procedure used to propagate input uncertainty showed that among the water quantity output variables, the overflow flow is the most uncertain output variable with a coefficient of variation (cv) of 1.59. Among water quality variables, the annual average spill COD oncentration and the average spill NH4 concentration were the most uncertain model outputs (coefficients of variation of 0.99 and 0.82, respectively). Also, low standard errors for the coefficient of variation were obtained for all seven outputs. These were never greater than 0.05, which indicates that the selected MC replication size (1,500 simulations) was sufficient. We also evaluated how uncertainty propagation can explain more comprehensively the impact of water quality indicators for the receiving river. While the mean model water quality outputs for COD and NH4 concentrations were slightly above the threshold, the 0.95 quantile was 2.7 times above the mean value for COD concentration, and 2.4 times above the mean value for NH4. This implies that there is a considerable probability that these concentrations in the spilled CSO are substantially larger than the threshold. However, COD and NH4 concentration levels of the river water will likely stay below the water quality threshold, due to rapid dilution after CSO spill enters the river.

Heavy metal contamination in soil is a major environmental issue intensified by rapid industrial and population growth. Understanding the spatial distribution of soil contamination by heavy metals in the ecosystem is a necessary precondition to monitor soil health and to assess the ecological risks. The main sources of heavy metals in soil are natural and anthropogenic sources. Natural sources are typically released of heavy metals from rock by weathering and atmospheric precipitation. Anthropogenic sources are related to industrialization, rapid urbanization, agricultural practices, and military activities. We analyzed a total of 358 topsoil samples (0–30 cm) collected in Golestan province in the northeast of Iran based on a regular square grid networks with 1,700 squares each sized 2.5 km²(random sampling within the grid). From these samples, we determined the spatial distribution of Cd, Cu, Ni, Zn, and Pb using random forest (RF). A multi-spectral image (Landsat 8), and environmental derivatives calculated from terrain attributes, climatic parameters, parent material, land use maps, distances to mine sectors, main roads, industrial sites, and rivers were used as covariates to predict the spatial distribution of concentrations of heavy metals. The multi-collinearity of the predictors was examined by the variance inflation factor (VIF), and a feature selection process (genetic algorithm) was applied to avoid noise and optimize the selected input variables for the final model. The predictive accuracy of RF model was assessed by the mean prediction error (ME), root mean squared error (RMSE), and coefficient of determination (R2) using 5-fold cross-validation technique. The results showed that the concentration levels (mg kg-1) of Cd, Cu, Pb, Ni, and Zn varied from 0.02 to 2.75, 9.70 to 93.70, 6.80 to 114.20, 9.50 to 93.20, and 25.10 to 417.4, respectively. The best prediction performance was for Ni (RMSE=9.9 mg kg-1 and R2=56.6%), and the lowest prediction performance for Cd (RMSE=0.4 mg kg-1 and R2=28.0%). Environmental covariates that control soil moisture and water flow along with climatic factors were the most important variables to define the spatial distribution of soil heavy metals. We conclude that the RF model using easily accessible environmental covariates is a promising, cost-effective and fast approach to monitor the spatial distribution of heavy metal contamination in soils.Keywords: Heavy metals; digital soil mapping; machine learning; random forest; spatial variation; soil pollution.

Spatially resolved estimates of change in soil organic carbon (SOC) stocks are necessary for supporting national and international policies aimed at achieving land degradation neutrality and climate mitigation through better land management. In this work we report on the development, implementation and application of a data-driven, statistical space-time method for mapping SOC stocks, using Argentina as a pilot area. We used the Quantile Regression Forest machine-learning algorithm to predict SOC stock at 0-30 cm depth at 250 m resolution for Argentina between 1982 and 2017, on an annual basis. The model was calibrated using over 5,000 SOC stock values from the 36-year time period and 35 environmental covariates. Most covariates were static and could only explain the spatial SOC distribution. SOC change over time was modelled using time series maps of the AVHRR NDVI vegetation index. These NDVI time series maps were pre-processed using a temporal low-pass filter to allow the SOC stock for a given year to depend on the NDVI of the current as well as preceding years. Spatial patterns of SOC stock predictions were persistent over time and comparable to baseline SOC stock maps of Argentina. Predictions had modest temporal variation with an average decrease for the entire country from 2.55 kg C m‑2 to 2.48 kg C m‑2 over the 36-year period (equivalent to a decline of 210.7 Gg C, 3.0% of the total 0‑30 cm SOC stock in Argentina). The Pampa region had a larger estimated SOC stock decrease from 4.62 kg C m‑2 to 4.34 kg C m‑2 (5.9%) during the same period. For the 2001-2015 period, predicted temporal variation was 7-fold larger than that obtained using the Tier 1 approach of the Intergovernmental Panel on Climate Change and the United Nations Convention to Combat Desertification. Prediction uncertainties turned out to be substantial, mainly due to the limited number and poor spatial and temporal distribution of the calibration data, and the limited explanatory power of the covariates. Cross-validation confirmed that SOC stock prediction accuracy was limited, with a Mean Error of 0.03 kg C m-2 and a Root Mean Squared Error of 2.04 kg C m-2. The model explained 45% of the SOC stock variation. In spite of the large uncertainties, this work showed that machine learning methods can be used for space-time SOC mapping and may yield valuable information to land managers and policy makers, provided that SOC observation density in space and time is sufficiently large.

Soil information is fundamental for many global applications, such as food security, land degradation, water resources, hydrology, climate change and ecological conservation. To address these diverse needs, it is important to provide free, consistent, easily accessible and standardized soil information. SoilGrids meets these requirements being a global product supporting global modelling and providing complementary information for the development of regional and national products in data-poor areas. This presentation will focus on the methodological aspects for modelling and mapping of global soil information. We describe the selection of models for global mapping using quantile random forest and recursive feature elimination to obtain a parsimonious model. We also use a refined cross-validation procedure to account for bias caused by spatial differences in sampling density at different depths. SoilGrids also quantifies location-specific uncertainty at global level by computing 90% prediction interval limits.

Regression kriging is one of the most popular spatial interpolation techniques. Its main strength is that it exploits both spatial autocorrelation as well as information contained in environmental covariates. While regression kriging is still dominant, in the past few years machine learning, especially random forest, is increasingly being used for mapping. Machine learning is more flexible than multiple linear regression and can thus make better use of environmental covariates. But machine learning typically ignores spatial autocorrelation. Several attempts have been made to include spatial autocorrelation in random forest, by adding distances to observation locations and other geometries to the set of covariates. But none of these studies have tried the obvious solution to include the nearest observations themselves and the distances to the nearest observations as covariates. In this study we tried this solution by introducing and testing Random Forest for Spatial Interpolation (RFSI). RFSI trains a random forest model on environmental covariates as well as nearest observations and their distances from the prediction point. We applied and evaluated RFSI for mapping daily precipitation in Catalonia for the 2016-2018 period. We trained four models (random forests, RFsp, pooled regression kriging and RFSI) using 63,927 daily precipitation observations from 87 GHCN-stations located in Catalonia. Maximum and minimum daily temperatures and IMERG precipitation estimates (derived from the GPM mission) were used as environmental covariates for all four models. Results based on 5-fold cross validation showed that RFSI (R-square 69.4%, RMSE 3.8 mm) significantly outperformed all random forest (R-square 50.6%, RMSE 3.8 mm), RFsp (R-square 55.5%, RMSE 4.6 mm) and pooled regression kriging (R-square 65.3%, RMSE 4.0 mm). Finetuning RFSI could potentially improve prediction accuracy even more. In addition to improved prediction accuracy, RFSI has the advantage that it uses much fewer spatial covariates than RFsp.

The soil organic carbon (SOC) pool is the largest terrestrial carbon (C) pool and is two to three times larger than the C stored in vegetation and the atmosphere. SOC is a crucial component within the C cycle, and an accurate baseline of SOC is required, especially for biogeochemical and earth system modelling. This baseline will allow better monitoring of SOC dynamics due to land use change and climate change. However, current estimates of SOC stock and its spatial distribution have large uncertainties. In this study, we test whether we can improve the accuracy of the three existing SOC maps of France obtained at national (IGCS), continental (LUCAS), and global (SoilGrids) scales using statistical model averaging approaches. Soil data from the French Soil Monitoring Network (RMQS) were used to calibrate and evaluate five model averaging approaches, i.e., Granger-Ramanathan, Bias-corrected Variance Weighted (BC-VW), Bayesian Modelling Averaging, Cubist and Residual-based Cubist. Cross-validation showed that with a calibration size larger than 100 observations, the five model averaging approaches performed better than individual SOC maps. The BC-VW approach performed best and is recommended for model averaging. Our results show that 200 calibration observations were an acceptable calibration strategy for model averaging in France, showing that a fairly small number of spatially stratified observations (sampling density of 1 sample per 2500 km2) provides sufficient calibration data. We also tested the use of model averaging in data-poor situations by reproducing national SOC maps using various sized subsets of the IGCS dataset for model calibration. The results show that model averaging always performs better than the national SOC map. However, the Modelling Efficiency dropped substantially when the national SOC map was excluded in model averaging. This indicates the necessity of including a national SOC map for model averaging, even if produced with a small dataset (i.e., 200 samples). This study provides a reference for data-poor countries to improve national SOC maps using existing continental and global SOC maps.

Countrywide estimates of soil organic carbon stock (SOCS) are useful to set up national strategies for sustainable land use management as well as to enhance the accuracy of global SOCS inventories. We appraised the spatial distribution of SOCS at five depth layers (0–15 cm, 15–30 cm, 30–100 cm, 0–30 cm and 0–100 cm) in Cameroon at 100 m spatial resolution, using a national harmonized legacy soil database (Camsodat 0.1) with 1432 georeferenced soil profiles. We assessed the prediction performances of random forest (RF) and generalized boosted regression (GBR), combined with two hybridization approaches of spatial interpolation of the residuals using ordinary kriging (OK) and inverse distance weighting (IDW). The estimates were compared to two global estimates derived from the Harmonized World Soil Database (HWSD) and SoilGrids250m. The SOCS distribution across the country showed a moderate spatial heterogeneity at all depth layers with coefficients of variation between 35% and 47%, and values ranging from 6 to 108 Mg C ha−1 at 0–15 cm, from 4 to 107 Mg C ha−1 at 15–30 cm, from 10 to 276 Mg C ha−1 at 30–100 cm, from 11 to 210 Mg C ha−1 at 0–30 cm and from 21 to 468 Mg C ha−1 at the 0–100 cm layer. Of the selected environmental covariates, terrain and climate attributes were the most relevant to predict the SOCS spatial distribution at country level. The RF model outperformed the GBR model, with about 10% improvement on prediction performance (R2) for most soil depths. The hybridization further slightly improved performance. However, OK was only slightly better than IDW in the overall assessment. Compared to national estimates, SoilGrids overestimated the SOCS by 15% at 0–30 cm depth, while HWSD underestimated SOCS by 26% at the same depth. Overall, about 5.7 Pg C are stored in the top 1 m of soils in Cameroon, with about 50% of that in the top 30 cm. The national distribution of SOCS is consistent with the pattern of agro-ecological zones. Our assessment provides baseline information for sustainable land management and climate change mitigation, as well as for improving the understanding of the spatial distribution of SOCS in Cameroon.

Nitrogen use efficiency (NUE) is crucial to establish efficient fertilizer application guidelines that balance crop yield, economic return and environmental sustainability. Although there are quite a few researches about the spatial and temporal variation of NUE, little work has been done on modelling NUE through deriving empirical relationships with explanatory environmental variables and exploring their relative importance quantitatively. The space-time patterns of NUE indicators (i.e., the Partial Factor Productivity of nitrogen, PFPN, and the Partial Nutrient Balance of nitrogen, PNBN) at provincial scale in China were derived and related to environmental covariates using stepwise multiple linear regression. PFPN was higher in east and south China than in central and west China and was smaller than 30 kg kg−1 yr−1 in most provinces, while PNBN was moderate in most provinces (0.41–0.50 kg kg−1 yr−1) and low (< 0.40 kg kg−1 yr−1) in south China. The national PFPN declined slightly from 32 kg kg−1 in 1978 to 27 kg kg−1 in 1995 and went up gradually to reach 38 kg kg−1 in 2015. The national PNBN decreased from 0.53 to 0.36 kg kg−1 from 1978 to 2003, thereafter stabilizing at around 0.40 kg kg−1 yr−1 between 2004 and 2015. The multiple linear regression models explained 74 % of the variance of PFPN and PNBN. The main explanatory variables of PFPN were planting area index of sugar crop (32 % of the R-square), followed by Arenosols (12 %), planting area index of oil crop (8 %), planting area index of vegetables (5 %), silt content (5 %) and total potassium (5 %). For PNBN, the variation was mainly attributed to mean annual daytime surface temperature (28 % of the R-square), planting area index of crops (beans 20 %, orchards 10 % and vegetables 9 %) and wet day frequency (5 %). The results of this study indicate that crop types, temperature and soil properties are important variables that determine NUE. These should be considered by policy makers when agricultural land development decisions are made in order to balance NUE and productivity (i.e., agronomy and environment).

For many decades, kriging and deterministic interpolation techniques, such as inverse distance weighting and nearest neighbour interpolation, have been the most popular spatial interpolation techniques. Kriging with external drift and regression kriging have become basic techniques that benefit both from spatial autocorrelation and covariate information. More recently, machine learning techniques, such as random forest and gradient boosting, have become increasingly popular and are now often used for spatial interpolation. Some attempts have been made to explicitly take the spatial component into account in machine learning, but so far, none of these approaches have taken the natural route of incorporating the nearest observations and their distances to the prediction location as covariates. In this research, we explored the value of including observations at the nearest locations and their distances from the prediction location by introducing Random Forest Spatial Interpolation (RFSI). We compared RFSI with deterministic interpolation methods, ordinary kriging, regression kriging, Random Forest and Random Forest for spatial prediction (RFsp) in three case studies. The first case study made use of synthetic data, i.e., simulations from normally distributed stationary random fields with a known semivariogram, for which ordinary kriging is known to be optimal. The second and third case studies evaluated the performance of the various interpolation methods using daily precipitation data for the 2016–2018 period in Catalonia, Spain, and mean daily temperature for the year 2008 in Croatia. Results of the synthetic case study showed that RFSI outperformed most simple deterministic interpolation techniques and had similar performance as inverse distance weighting and RFsp. As expected, kriging was the most accurate technique in the synthetic case study. In the precipitation and temperature case studies, RFSI mostly outperformed regression kriging, inverse distance weighting, random forest, and RFsp. Moreover, RFSI was substantially faster than RFsp, particularly when the training dataset was large and high-resolution prediction maps were made.

Digital soil mapping techniques usually take an entirely data-driven approach and model soil properties individually and layer by layer, without consideration of interactions. In previous studies we implemented a structural equation modelling (SEM) approach to include pedological knowledge and between-properties and between-layer interactions in the mapping process. However, it typically does not consider spatial correlation. Our goal was to extend SEM by accounting for residual spatial correlation using a geostatistical approach. We assumed second-order stationary and estimated the semivariogram parameters, together with the usual SEM parameters, using maximum likelihood estimation. Spatial prediction was done using regression kriging. The methodology is applied to mapping cation exchange capacity, clay content and soil organic carbon for three soil horizons in a 150100-km 2 study area in the Great Plains of the United States. The calibration process included all parameters used in lavaan , a SEM software, plus two extra parameters to model residual spatial correlation. The residuals showed substantial spatial correlation, which indicates that including spatial correlation yields more accurate predictions. We also compared the standard SEM and the spatial SEM approaches in terms of SEM model coefficients. Differences were substantial but none of the coefficients changed sign. Presence of residual spatial correlation suggests that some of the causal factors that explain soil variation were not captured by the set of covariates. Thus, it is worthwhile to search for additional covariates leaving only unstructured residual noise, but provided that this is not achieved, it is beneficial to include residual spatial correlation in mapping using SEM.

The functional potential of soil ecosystems can be predicted from the activity and abundance of the microbial community in relation to key soil properties. When describing microbial community dynamics, soil physicochemical properties have traditionally been used. The extent of correlations between properties, however, differs between studies, especially across larger spatial scales. In this research we analysed soil microbial biomass and substrate-induced respiration of 156 samples from Irish grasslands. In addition to the standard physicochemical, soil type and land management variables, soil diagnostic properties were included to identify if these important soil-landscape genesis classes affected microbial biomass and respiration dynamics in Irish soil. Apart from physicochemical properties, soil drainage class was identified as having an important effect on microbial properties. In particular, biomass-specific basal (qCO2) and substrate-induced respiration (SIR:CFE) were explained best by the soil drainage. Poorly drained soil had smaller values of these respiration measures than well-drained soil. We concluded that this resulted from different groups within the microbial community that could use readily available carbon sources, which suggests a change in microbial community dynamics associated with soil texture and periods of water stress. Overall, our results indicate that soil quality assessments should include both physicochemical properties and diagnostic classes, to provide a better understanding of the behaviour of soil microbial communities. Highlights: Assessing the effect of soil diagnostic features and properties on microbial biomass and respiration A soil biological survey from 156 grassland sites in Ireland Soil drainage class has an important effect on microbial properties Soil quality assessments should include both physicochemical properties and diagnostic classes

An approach for using lasso (Least Absolute Shrinkage and Selection Operator) regression in creating sparse 3D models of soil properties for spatial prediction at multiple depths is presented. Modeling soil properties in 3D benefits from interactions of spatial predictors with soil depth and its polynomial expansion, which yields a large number of model variables (and corresponding model parameters). Lasso is able to perform variable selection, hence reducing the number of model parameters and making the model more easily interpretable. This also prevents overfitting, which makes the model more accurate. The presented approach was tested using four variable selection approaches – none, stepwise, lasso and hierarchical lasso, on four kinds of models – standard linear model, linear model with polynomial expansion of depth, linear model with interactions of covariates with depth and linear model with interactions of covariates with depth and its polynomial expansion. This framework was used to predict Soil Organic Carbon (SOC) in three contrasting study areas: Bor (Serbia), Edgeroi (Australia) and the Netherlands. Results show that lasso yields substantial improvements in accuracy over standard and stepwise regression — up to 50 % of total variance. It yields models which contain up to five times less nonzero parameters than the full models and that are usually more sparse than models obtained by stepwise regression, up to three times. Extension of the standard linear model by including interactions typically improves the accuracy of models produced by lasso, but is detrimental to standard and stepwise regression. Regarding computation time, it was demonstrated that lasso is several orders of magnitude more efficient than stepwise regression for models with tens or hundreds of variables (including interactions). Proper model evaluation is emphasized. Considering the fact that lasso requires meta-parameter tuning, standard cross-validation does not suffice for adequate model evaluation, hence a nested cross-validation was employed. The presented approach is implemented as publicly available sparsereg3D R package.

Simple and ordinary kriging assume a constant mean and variance of the soil variable of interest. This assumption is often implausible because the mean and/or variance are linked to terrain attributes, parent material or other soil forming factors. In kriging with external drift (KED) non-stationarity in the mean is accounted for by modelling it as a linear combination of covariates. In this study, we applied an extension of KED that also accounts for non-stationary variance. Similar to the mean, the variance is modelled as a linear combination of covariates. The set of covariates for the mean may differ from the set for the variance. The best combinations of covariates for the mean and variance are selected using Akaike's information criterion. Model parameters of the selected model are then estimated by differential evolution using the Restricted Maximum Likelihood (REML) in the objective function. The methodology was tested in a small area of the Hunter Valley, NSW Australia, where samples from a fine grid with gamma K measurements were treated as measurements of the variable of interest. Terrain attributes were used as covariates. Both a non-stationary variance and a stationary variance model were calibrated. The mean squared prediction errors of the two models were somewhat comparable. However, the uncertainty about the predictions was much better quantified by the non-stationary variance model, as indicated by the mean and median of the standardized squared prediction error and by accuracy plots. We conclude that the non-stationary variance model is more flexible and better suited for uncertainty quantification of a mapped soil property. However, parameter estimation of the non-stationary variance model requires more attention due to possible singularity of the covariance matrix.

Superresolution mapping (SRM) is a widely used technique to address the mixed pixel problem in pixel-based classification. Advanced object-based classification will face a similar mixed phenomenon--a mixed object that contains different land-cover classes. Currently, most SRM approaches focus on estimating the spatial location of classes within mixed pixels in pixel-based classification. Little if any consideration has been given to predicting where classes spatially distribute within mixed objects. This paper, therefore, proposes a new object-based SRM strategy (OSRM) to deal with mixed objects in object-based classification. First, it uses the deconvolution technique to estimate the semivariograms at target subpixel scale from the class proportions of irregular objects. Then, an area-to-point kriging method is applied to predict the soft class values of subpixels within each object according to the estimated semivariograms and the class proportions of objects. Finally, a linear optimization model at object level is built to determine the optimal class labels of subpixels within each object. Two synthetic images and a real remote sensing image were used to evaluate the performance of OSRM. The experimental results demonstrated that OSRM generated more land-cover details within mixed objects than did the traditional object-based hard classification and performed better than an existing pixel-based SRM method. Hence, OSRM provides a valuable solution to mixed objects in object-based classification.

Summary In a previous study we introduced structural equation modelling (SEM) for digital soil mapping in the Argentine Pampas. An attractive property of SEM is that it incorporates pedological knowledge explicitly through a mathematical implementation of a conceptual model. Many soil processes operate within the soil profile; therefore, SEM might be suitable for simultaneous prediction of soil properties for multiple soil layers. In this way, relations between soil properties in different horizons can be included that might result in more consistent predictions. The objectives of this study were therefore to apply SEM to multi-layer and multivariate soil mapping, and to test SEM functionality for suggestions to improve the modelling. We applied SEM to model and predict the lateral and vertical distribution of the cation exchange capacity (CEC), organic carbon (OC) and clay content of three major soil horizons, A, B and C, for a 23 000-km2 region in the Argentine Pampas. We developed a conceptual model based on pedological hypotheses. Next, we derived a mathematical model and calibrated it with environmental covariates and soil data from 320 soil profiles. Cross-validation of predicted soil properties showed that SEM explained only marginally more of the variance than a linear regression model. However, assessment of the covariation showed that SEM reproduces the covariance between variables much more accurately than linear regression. We concluded that SEM can be used to predict several soil properties in multiple layers by considering the interrelations between soil properties and layers. Highlights We tested structural equation modelling (SEM) for multi-layer and multivariate soil mapping. SEM models soil property covariation better than multiple linear regression. The SEM re-specification step improves prediction accuracy. SEM supports learning about soil processes from data.

The accuracy of spatial predictions of rainfall by merging rain-gauge and radar data is partly determined by the sampling design of the rain-gauge network. Optimising the locations of the rain-gauges may increase the accuracy of the predictions. Existing spatial sampling design optimisation methods are based on minimisation of the spatially averaged prediction error variance under the assumption of intrinsic stationarity. Over the past years, substantial progress has been made to deal with non-stationary spatial processes in kriging. Various well-documented geostatistical models relax the assumption of stationarity in the mean, while recent studies show the importance of considering non-stationarity in the variance for environmental processes occurring in complex landscapes. We optimised the sampling locations of rain-gauges using an extension of the Kriging with External Drift (KED) model for prediction of rainfall fields. The model incorporates both non-stationarity in the mean and in the variance, which are modelled as functions of external covariates such as radar imagery, distance to radar station and radar beam blockage. Spatial predictions are made repeatedly over time, each time recalibrating the model. The space-time averaged KED variance was minimised by Spatial Simulated Annealing (SSA). The methodology was tested using a case study predicting daily rainfall in the north of England for a one-year period. Results show that (i) the proposed non-stationary variance model outperforms the stationary variance model, and (ii) a small but significant decrease of the rainfall prediction error variance is obtained with the optimised rain-gauge network. In particular, it pays off to place rain-gauges at locations where the radar imagery is inaccurate, while keeping the distribution over the study area sufficiently uniform.