Articles | Volume 16, issue 6
Data description paper
25 Jun 2024
Data description paper |  | 25 Jun 2024

BIS-4D: mapping soil properties and their uncertainties at 25 m resolution in the Netherlands

Anatol Helfenstein, Vera L. Mulder, Mirjam J. D. Hack-ten Broeke, Maarten van Doorn, Kees Teuling, Dennis J. J. Walvoort, and Gerard B. M. Heuvelink

In response to the growing societal awareness of the critical role of healthy soils, there has been an increasing demand for accurate and high-resolution soil information to inform national policies and support sustainable land management decisions. Despite advancements in digital soil mapping and initiatives like GlobalSoilMap, quantifying soil variability and its uncertainty across space, depth and time remains a challenge. Therefore, maps of key soil properties are often still missing on a national scale, which is also the case in the Netherlands. To meet this challenge and fill this data gap, we introduce BIS-4D, a high-resolution soil modeling and mapping platform for the Netherlands. BIS-4D delivers maps of soil texture (clay, silt and sand content), bulk density, pH, total nitrogen, oxalate-extractable phosphorus, cation exchange capacity and their uncertainties at 25 m resolution between 0 and 2 m depth in 3D space. Additionally, it provides maps of soil organic matter and its uncertainty in 3D space and time between 1953 and 2023 at the same resolution and depth range. The statistical model uses machine learning informed by soil observations amounting to between 3815 and 855 950, depending on the soil property, and 366 environmental covariates. We assess the accuracy of mean and median predictions using design-based statistical inference of a probability sample and location-grouped 10-fold cross validation (CV) and prediction uncertainty using the prediction interval coverage probability.

We found that the accuracy of clay, sand and pH maps was the highest, with the model efficiency coefficient (MEC) ranging between 0.6 and 0.92 depending on depth. Silt, bulk density, soil organic matter, total nitrogen and cation exchange capacity (MEC of 0.27 to 0.78), and especially oxalate-extractable phosphorus (MEC of −0.11 to 0.38) were more difficult to predict. One of the main limitations of BIS-4D is that prediction maps cannot be used to quantify the uncertainty in spatial aggregates. We provide an example of good practice to help users decide whether BIS-4D is suitable for their intended purpose. An overview of all maps and their uncertainties can be found in the Supplement. Openly available code and input data enhance reproducibility and help with future updates. BIS-4D prediction maps can be readily downloaded at (Helfenstein et al.2024a). BIS-4D fills the previous data gap of the national-scale GlobalSoilMap product in the Netherlands and will hopefully facilitate the inclusion of soil spatial variability as a routine and integral part of decision support systems.

1 Introduction

Life on Earth, including that of humans, relies fundamentally on the availability and quality of air, water and soil. These essential resources exhibit spatial variations in accordance with Tobler's first law of geography, asserting that “everything is related to everything else, but near things are more related than distant things” (Tobler1970). However, the spatial heterogeneity of soil properties stands out prominently over short distances compared to air and water. This disparity arises from the multifaceted nature of soil, comprising solid, liquid and gaseous phases, making it less mobile and unable to form homogeneous mixtures akin to air or water. Moreover, soil formation is a gradual process unfolding over hundreds to millions of years, shaped by intricate interactions between the climate, organisms (including humans), topography and parent material (Dokuchaev1899; Jenny1941). Some of these soil-forming factors themselves exhibit high heterogeneity over short distances. Fully grasping soil spatial variability requires dense sampling, but this is hindered by the difficulty, time and expense of collecting soil samples. These challenges underscore the complexity of quantifying soil variation, highlighting the formidable task of mapping soils in 3D space and time (3D+T).

With the rising awareness of soil health among diverse stakeholders and within value chains (Lehmann et al.2020), soil scientists have been increasingly dedicated to delivering high-resolution, accurate soil maps. Internationally prominent examples of policies for which spatio-temporal soil information is essential include several of the Sustainable Development Goals, such as zero hunger and life on land (United Nations2015), and, in Europe, the Green Deal, Common Agricultural Policy and Zero Pollution Action Plan (Panagos et al.2022). The importance of soil information for these policies has led to the EU soil strategy for 2030; the Soil Deal (European Commission2021); and, most recently, the Proposal for a Directive on Soil Monitoring and Resilience (European Commission2023). For such policies to have an impact, it is essential that soil scientists deliver information required to facilitate land use decisions and management practices at multiple scales.

In the Netherlands (land area = 33 481 km2), the demand for soil information is also large. Located in the midst of Europe's largest delta, soils in the Netherlands are naturally very fertile (Edelmann1950; Römkens and Oenema2004). As one of Europe's most densely populated countries, multi-functional land use decisions made at national or regional level need to be implemented at the field level, involving a broad range of diverse stakeholders. This spectrum of stakeholders collaborates on initiatives like the Smart Land Use project, which aims to sequester an additional 0.5 Mt CO2-eq per year to Dutch mineral agricultural soils (Slier et al.2023). Spatial information on soil properties can be used to evaluate soil health on Dutch agricultural fields using tools such as the Open Soil Index (OSI; Ros et al.2022; Ros2023) and Soil Indicators for Agriculture (BLN 2.0; Ros et al.2023) and for assessing soil functions at different scales (Schulte et al.2015). Information on soil texture and soil organic matter (SOM) is necessary for greenhouse gas reporting of the Land Use, Land Use Change and Forestry (LULUCF) sector for the United Nations Framework Convention on Climate Change and the Dutch LULUCF submission under the Kyoto Protocol (KP-LULUCF; Arets et al.2020). Data on basic soil properties serve as inputs for modeling agricultural suitability (Mulder et al.2022), precision agriculture (Been et al.2023) and soil–water–atmosphere–plant interactions (SWAP; van Dam et al.1997; Kroes et al.2017). Furthermore, soil property maps contribute to initiatives such as the WaterVision Agriculture and Nature (Hack-ten Broeke et al.2019), Hydrological Instrumentations of the Netherlands (NHI2023), and Delta Program 2024 (Delta Programme2023).

Soil maps can also be used to identify and prioritize threats to soil health, as reviewed for the Netherlands by Römkens and Oenema (2004) and Hack-ten Broeke et al. (2009). Specific threats to soil health in the Netherlands include soil compaction (van den Akker and Hoogland2011; van den Akker et al.2012), subsidence of peat due to oxidation and compaction (Brouwer et al.2018; van Asselen et al.2018), subsidence of young clay soils due to ripening on reclaimed land (Brouwer et al.2018), and soil erosion (Hessel et al.2011). Recently, Helfenstein et al. (2024c) mapped SOM in 3D+T, which identified decreases in SOM at high resolution in 3D space. Spatial soil information is also crucial for not only agricultural businesses, for optimizing both fertilizer and manure applications for crop growth, but also environmental accounting. The demand for such information is especially high in the Netherlands (Stokstad2019; Erisman2021; Aarts and Leeuwis2023), as it has the highest livestock density in the EU (Eurostat2022, p. 32) and ranks as the world's second-largest agricultural exporter (Jukema et al.2023). An estimated 1 300 000 ha are phosphate saturated soils, where phosphate loss due to leaching exceeds ecologically tolerable limits (Römkens and Oenema2004). Hence, providing spatially explicit soil information is crucial to adhere to targets 4.2 and 4.3 of the Soil Deal for Europe, which aim to reduce fertilizer use by at least 20 % and reduce nutrient losses by at least 50 % by 2030 (European Commission2021). In summary, the pressure of using soils sustainably in the Netherlands is immense.

Between the 1950s and 2000, conventional soil maps were completed in many countries. Today, the well-established discipline of digital soil mapping (DSM) has been widely adopted to meet the demands for accurate and high-resolution soil information for a wide range of purposes. Since DSM was first conceptualized (McBratney et al.2003; Scull et al.2003), maps of soil properties and soil types have been produced from local to global scales. These advances were propelled by initiatives like GlobalSoilMap (GSM) under the support of the International Union of Soil Sciences (Arrouays et al.2014a, b, 2015) and the availability of openly accessible tutorials elucidating standard DSM workflows (Malone et al.2017; Hengl and MacMillan2019; Brus et al.2017; Brus2019, 2022).

Historically, the Netherlands was at the forefront of soil mapping. Scientific soil investigations in the Netherlands were started by Winand C. H. Staring in the mid-1800s followed by Jan van Baren and David J. Hissink in the early 1900s (Bouma and Hartemink2003). The first publication on the spatial distribution of soil properties in the Netherlands dates back to the 19th century (Felix1995). Systematic soil mapping became institutionalized with the establishment of the Dutch Soil Survey institute or Stichting voor Bodemkartering (StiBoKa) in 1945 (Hartemink and Sonneveld2013). From 1950 to 1995, StiBoKa conducted conventional soil surveys (Buringh et al.1962; de Bakker and Schelling1989; ten Cate et al.1995) and produced regional maps (1:10 000 and 1:25 000 scale) and a national map (1:50 000 scale) of soil types (de Vries et al.2003). After the development of DSM as a research field, various studies used (geo-)statistical methods to map qualitative and quantitative soil properties using the data collected by StiBoKa (Brus and Heuvelink2007; Brus et al.2009; Kempen et al.2014; van den Berg et al.2017; Helfenstein et al.2022, 2024c). Several regions of the national soil map have since been updated (Kempen et al.2009, 2011, 2012a; de Vries et al.2014, 2017, 2018; Brouwer et al.2018; Brouwer and Walvoort2019, 2020; Brouwer et al.2021, 2023) and a variety of thematic maps were derived, such as a map of re-worked soils (Brouwer and van der Werff2012), a peat thickness map (Brouwer et al.2018), a map of soil landscapes (van Delft and Maas2022, 2023) and the soil physical units map of the Netherlands (BOFEK; Heinen et al.2022).

DSM has established itself and is routinely implemented across the world, but various challenges remain (Chen et al.2022; Wadoux et al.2021b). Maps of basic chemical, physical and especially biological soil properties are often missing (Chen et al.2022; Wadoux et al.2021b). Approximately 78 % of articles reviewed by Chen et al. (2022) mapped SOM, carbon content and carbon stocks. If a DSM product is available, predictions are often only made for one depth layer. Half of the studies reviewed by Chen et al. (2022) focused on soil properties at less than 30 cm depth only. However, users also require soil information at greater depths and could benefit from models being able to predict at any desired depth in 3D and, for dynamic soil properties, in 3D+T (Chen et al.2022; Wadoux et al.2021b).

There are numerous challenges related to the accuracy of soil maps (Wadoux et al.2021b). One major challenge is that the uncertainty in soil maps is often not quantified. A recent review showed that only 35 % of studies mapping continuous soil properties estimated prediction uncertainty (Piikki et al.2021). Without providing the uncertainty in a map, users cannot determine its fitness for use. Moreover, assessing map accuracy is not straightforward and involves many demanding pre-requisites – for example, the sampling design of the locations used for statistical validation. According to Piikki et al. (2021), only 13 % of studies used probability sampling for map validation, which is the best approach for assessing map accuracy (Brus et al.2011; Wadoux et al.2021a; de Bruin et al.2022). When using a soil map in a model or analysis, the uncertainty may be so large that it compromises the quality of the outputs of the model or analysis, posing risks of erroneous conclusions and decisions for end users (Knotters and Vroon2015; Knotters et al.2015a, b; Heuvelink2018). The efficacy of uncertainty propagation analysis relies on quantifying input uncertainty realistically, emphasizing the need to quantify uncertainty in soil maps. The above challenges also apply to the Netherlands, where there is not yet a product that meets all these requirements.

To meet these challenges and demands, we introduce a high-resolution soil modeling and mapping platform for the Netherlands called BIS-4D (Fig. 1). It delivers maps of key soil properties according to GSM specifications and assesses their accuracy using prediction uncertainty and statistical validation. The platform provides maps of soil texture (clay, silt and sand content), bulk density (BD), pH, total nitrogen (Ntot), oxalate-extractable phosphorus (Pox) and cation exchange capacity (CEC) at 25 m resolution between 0 and 2 m depth in 3D space (Table 1). Furthermore, we provide maps of SOM in 3D+T between 1953 and 2023 at the same resolution and depth range since SOM has changed substantially over time. Note that for soil pH and SOM, specific updates were made compared to previous versions (Helfenstein et al.2022, 2024c, Sect. 2.7). These nine soil properties were chosen based on those prioritized by GSM (Arrouays et al.2014a, b, 2015), end-user needs in the Netherlands and data availability. In collaboration with soil surveyors, database maintainers, and experts on Dutch soils from Wageningen University & Research, we assessed the strengths and limitations of the BIS-4D maps and recommended potential map applications. Finally, model inputs, outputs (BIS-4D maps) and code, using free and open-source software, were made available, easily accessible and well documented so that BIS-4D can be updated for future applications.

Figure 1Graphical abstract of the BIS-4D soil modeling and mapping platform, where Y is a target soil property and X are covariates that vary in 2D space (s), depth (d) and, for SOM, time (t). CLORPT stands for the soil-forming factors, i.e., climate, organisms, relief, parent material and time (Dokuchaev1899; Jenny1941). RFE: recursive feature elimination, QRF: quantile regression forest, PI90: 90th prediction interval width, PICP: prediction interval coverage probability (Sect. 2.32.6).


(NEN 57532020)ten Cate et al. (1995)de Bakker and Schelling (1966)de Bakker and Schelling (1989)(NEN 57532020)(NEN 57532020)ten Cate et al. (1995)de Bakker and Schelling (1966)de Bakker and Schelling (1989)(Maring et al.2009, Appendix E, p. 79)(Maring et al.2009, Appendix E, p. 81)(Maring et al.2009, Appendix E, p. 81)

Table 1Abbreviations, units and description of methods used for laboratory measurements and field estimates of target soil properties. Mineral soil is defined as the dried soil fraction (105 °C) put through a 2 mm sieve after the removal of SOM and CaCO3.

Download Print Version | Download XLSX

2 Materials and Methods

We predicted soil properties Y^ in 3D space and SOM in 3D+T using well-established DSM methods (Fig. 1). BIS-4D uses machine learning to model the relationship between a soil property measured or estimated in the field at point locations as the model response Y (Tables 13) and environmental covariates as the explanatory variables X (Table 5).

2.1 Soil point data

BIS-4D uses laboratory measurements and field estimates of soil properties from point locations collected in the Dutch soil database, or Bodemkundig informatie systeem (BIS,, last access: 23 January 2024). Definitions and laboratory measurement and field estimation methods for the soil properties mapped using BIS-4D are described in Table 1. We only included observations between 0 and 2 m depth excluding the O horizon or the layer with dead plant material, leaves, branches and other decomposing organic material on top of mineral soils. As the majority of the soil point data were collected before modern Global Positioning Systems (Table 2), soil surveyors marked the point locations on a 1:25 000 topographic map.

Note that clay, silt and sand content are particle size fractions (PSFs) which together constitute soil texture. Thus, soil texture is a compositional variable: each PSF must be non-negative and together they must add up to 100 % (Pawlowsky-Glahn and Buccianti2011; Pawlowsky-Glahn et al.2015). In order to achieve this, soil texture can be spatially interpolated as a compositional variable using geostatistical models (Odeh et al.2003; Lark and Bishop2007; Wang and Shi2017), e.g., compositional kriging (de Gruijter et al.1997; Walvoort and de Gruijter2001), machine learning (Akpa et al.2014; Amirian-Chakan et al.2019; Poggio and Gimona2017; Poggio et al.2021; Malone et al.2021; Varón-Ramírez et al.2022) and other techniques (Buchanan et al.2012; Román Dobarco et al.2017). Most commonly, these studies used the additive log-ratio transformation with the Gauss–Hermite quadrature (Aitchison1986). When not modeled as a compositional variable, other approaches include estimating two of the three PSFs and calculating the third by subtracting the sum of the two estimates from 100 % (Adhikari et al.2013) or modeling all three PSFs separately (Viscarra Rossel et al.2015; Chagas et al.2016; Mulder et al.2016; Taghizadeh-mehrjardi et al.2016; Pahlavan-Rad and Akbarimoghaddam2018) and post-processing the predictions to ensure that they are all non-negative and sum to 100 %. For BIS-4D, we decided to model PSFs separately followed by post-processing (Sect. 2.5) for three reasons. Firstly, we wanted to use the additional locations where only one or two PSFs were observed (Table 2). Secondly, modeling soil texture as a compositional variable does not necessarily improve model performance (Amirian-Chakan et al.2019). Thirdly, modeling separately followed by post-processing is easy to implement.

2.1.1 Soil point data for model calibration

We used laboratory measurements and field estimates from the Boring Bodemkundig pakket (BPK) and Profielbeschrijving (PFB) datasets in BIS for model selection, tuning and calibration (Tables 2 and 3, Fig. 3). Observations in BPK and PFB were made by soil horizon and the range of years of soil sampling are given in Table 2. Laboratory measurements and field estimates were available for all depths between 0 and 2 m (Table 3). All laboratory measurements were made at PFB locations. These locations are arranged in a purposive sampling design selected in the past to create the national 1:50 000-scale soil type map (de Vries et al.2003). For the majority of the target soil properties, these locations covered soil variability in the Netherlands well (Fig. 2). The majority of field estimates are part of the BPK dataset and are spatially clustered in specific areas for regional soil mapping purposes or specific projects (Fig. 2 in Helfenstein et al.2024c). Most soil properties follow a skewed distribution, especially SOM, Ntot, Pox and CEC (Fig. 3). However, pH, sand and, to a lesser extent, silt followed bimodal distributions. The distributions of the target soil properties likely affected model predictions (Sect. 4.1).

Table 2Descriptive statistics of soil point data used for model calibration (field estimates and laboratory measurements) across all depths. Locations: number of locations, Obs.: number of observations, min: minimum, max: maximum, year: years during which observations were made, lab: laboratory measurements, field: field estimates. Minimum, median, mean and maximum values are in units of measurement of each soil property (Table 1). Soil point data used for model calibration are publicly available (see Sect. 3).

Download Print Version | Download XLSX

Figure 2Observation density for locations with laboratory measurements used for model calibration of all BIS-4D target soil properties. All of these locations are part of the PFB dataset.

Table 3Number of laboratory measurements (lab) and field estimates (field) used for model calibration per standard GSM depth layer for each soil property.

Download Print Version | Download XLSX

Figure 3Histograms of soil property observations used for model calibration, colored by observation type.


The laboratory measurements were deemed more important than field estimates because they are more accurate and locations with laboratory measurements were less spatially clustered. Nevertheless, field estimates from BPK and PFB also provide valuable information, expanding spatial coverage and, for SOM, temporal coverage from 1953–2022 as well (Table 2). In addition, since around 2000, most observations that were added to the BIS are field estimates, a trend which is likely to continue into the future due to limited budgets for laboratory measurements. Other national mapping studies have also used field estimates in the past (van den Berg et al.2017). We accounted for differences in data quality between laboratory measurements and field estimates using rigorous model tuning based on optimizing model performance (Sect. 2.3). Field estimates were removed if there was a laboratory measurement available from the same location and soil horizon (and year in case of SOM). Methods for estimating clay content, BD and SOM in the field are described in ten Cate et al. (1995), de Bakker and Schelling (1966) and de Bakker and Schelling (1989).

2.1.2 Soil point data for statistical validation

For clay, silt, sand and CEC, no separate dataset with laboratory measurements was available for statistical validation, meaning all observations were used for model calibration. Therefore, statistical validation of these four soil properties was conducted using cross validation (CV) of PFB laboratory measurements (Sect. 2.6).

For BD, pH, SOM, Ntot and Pox, laboratory measurements from either the Landelijke Steekproef Kaarteenheden (LSK) or the Carbon Content NL (CCNL) dataset were available for model validation (Table 4). LSK is a separate and independent dataset gathered between 1993 and 2000, where locations were determined using probability sampling. The LSK locations are shown in Fig. 1 in Helfenstein et al. (2022). The stratified simple random sample contains 94 strata defined based on soil type and groundwater class (Finke et al.2001; Visschers et al.2007), with the original purpose being to validate the national soil type map (de Vries et al.2003). Observations were made for each soil horizon. Statistical validation of BD, pH, SOM, Ntot and Pox maps was conducted using LSK because map accuracy should preferably be estimated with design-based statistical inference using a probability sample (Brus et al.2011). LSK data were also used to validate earlier versions of soil pH (Helfenstein et al.2022) and SOM maps (Helfenstein et al.2024c).

Table 4Descriptive statistics of separate soil point datasets used for statistical validation across all depths. Note that for statistical validation, only laboratory measurements were used. Separate datasets were not available for clay, silt, sand and CEC. Locations: number of locations, obs.: number of observations, min: minimum, max: maximum, year: periods during which observations were made.

Download Print Version | Download XLSX

For SOM and Ntot, the CCNL dataset was used for statistical validation (Table 4). The CCNL dataset consists of laboratory measurements from re-visited LSK locations in 2018, excluding locations that were no longer accessible. In contrast to LSK, during which soil samples were taken by soil horizon, CCNL locations were re-sampled at two fixed depth layers (0–30 and 30–100 cm). LSK and CCNL datasets were also used and their methodological sampling differences were explained in van Tol-Leenders et al. (2019), van den Elsen et al. (2020), and Knotters et al. (2022). Since LSK was sampled by soil horizon at more locations and also below 1 m depth, it is preferable to use it rather than CCNL.

For 3D+T maps of SOM, four different datasets were used for statistical validation with the specific purpose to assess SOM maps for specific years (Helfenstein et al.2024c): location-grouped 10-fold cross validation of PFB data (1953–2011; lab measurements shown in Table 2 and Fig. 3), design-based inference using LSK (1993–2000), design-based inference using CCNL (2018), and a separate set of PFB locations that were re-sampled in 2022 and used to assess changes in SOM over time (Table 4). Design-based inference and cross-validation procedures are explained in Sect. 2.6.

2.2 Covariates

In line with the DSM methodology (McBratney et al.2003; Scull et al.2003), we used 366 covariates as explanatory variables that were representative of the soil-forming factors: climate, organisms, relief (topography), parent material (geology) and time (Dokuchaev1899; Jenny1941). Accounting for Tobler's first law of Geography (Tobler1970) and spatial autocorrelation, easting (x coordinate) and northing (y coordinate) were also included as covariates. Numerous studies have used spatial position and geographical distances as covariates (Li et al.2011; Behrens et al.2018b; Hengl et al.2018; Møller et al.2020; Sekulić et al.2020). Sampling depth information, more specifically the upper and lower boundary and midpoint of each sampled horizon, were included as covariates so that predictions could be made at any chosen depth and depth interval. See Ma et al. (2021) for an overview of models using depth as a covariate in comparison to non-3D DSM methods. The majority of static covariates used in BIS-4D were previously used to map soil pH (Helfenstein et al.2022). Others, mainly derivations of monthly mosaics from Sentinel-2 RGB and NIR bands, were added to map SOM (Helfenstein et al.2024c). In order to map SOM in 3D+T, we extended upon established methods by also deriving the covariate variable in time (2D+T) and variable over depth and time (3D+T), as described in detail in Helfenstein et al. (2024c). All covariates were re-sampled at 25 m resolution according to the resolution of the national land use maps (WENR2020; Hazeu et al.2020).

We created a regression matrix containing the BIS-4D target soil property observations and static covariate values by performing a spatial overlay. For SOM, this was extended to a space–time overlay for 2D+T covariates and a space–depth–time overlay for 3D+T covariates (Helfenstein et al.2024c).

2.3 Model selection, tuning and calibration

For model selection as defined by Hastie et al. (2009), we removed covariates in a two-step procedure using decorrelation followed by recursive feature elimination (RFE), as in Poggio et al. (2021). From any pair of covariates for which the Pearson correlation coefficient was >0.85 or <-0.85, the covariate that was more correlated with all remaining covariates was removed. RFE (Guyon et al.2002) was implemented using the caret package (Kuhn2019) and the number of covariates was chosen with the lowest root mean squared error (RMSE; Eq. 3). From 366 covariates, this resulted in a set of 50, 30, 20, 15 or 10 covariates depending on the target soil property (Table 5), further used in model tuning, calibration and prediction.

(de Vries et al.2003)de Gruijter et al. (2004)Hoogland et al. (2014)Knotters et al. (2018)KNMI (2020)KNMI (2020)(Alterra2004)(WENR2020; Hazeu et al.2020)Loiseau et al. (2019)Roerink and Mücher (2023)(BGDM; RIVM2020)BIJ12 (2019)Bakker et al. (1989)de Vries and Al (1992)Clement (2001)Maas et al. (2019)AHN (2023)Knotters et al. (2018)Kombrink et al. (2012)van der Meulen et al. (2013)Koomen and Maas (2004)Maas et al. (2019)EZK (2013)Vos (2015)Vos et al. (2020)(Helfenstein et al.2024c)(Alterra2004)(WENR2020; Hazeu et al.2020)(Helfenstein et al.2024c)(de Vries et al.2003)

Table 5Covariates used during model calibration and prediction for different responses (soil properties), i.e., after covariate removal based on decorrelation and recursive feature elimination (RFE; Sect. 2.3). “All” implies that a covariate was used in tuning, calibration and prediction of all soil properties. Further information can be found in the metadata files and description of the provided covariates (Sect. 3).

Download Print Version | Download XLSX

For model tuning, we grew random forest (RF) models (Breiman2001) and optimized hyperparameters for mean predictions. We tuned the model using a location-grouped 10-fold cross validation, meaning that all measurements from the same soil profile location were forced to be in the same fold. Location-grouped cross validation was chosen because observations from the same profile in both model training and validation can lead to overly optimistic model accuracy metrics. Field estimates were excluded from the holdout fold. We assessed all combinations of the same hyperparameters as in Sect. 2.4 of Helfenstein et al. (2022) and chose the combination with the lowest RMSE (Eq. 3; Table 6).

For soil properties where both laboratory and field estimates were available (clay, silt, sand, BD and SOM), we also tuned whether designating a larger case weight for laboratory measurements improved model performance in order to account for the lower accuracy of field estimates compared to laboratory measurements. Values of 2, 5, 10 and 15 times the weight of field estimates were tested for laboratory measurements (Table 6). In addition, we also tested excluding field estimates entirely. The final set of hyper parameters was chosen based on the lowest RMSE (Eq. 3) across the cross validation. When the increase in RMSE was below 0.1 %, the model with fewer trees was chosen to reduce computation time during prediction. For silt and sand, model performance was the highest when using only laboratory measurements, so field estimates were excluded in model calibration (Table 6).

For model calibration and prediction, we used RF to predict the mean and quantile regression forest (QRF) due to its ability to predict the entire conditional distribution (Meinshausen2006). The final models used for prediction were fitted using all soil observations in the calibration set (Table 2), the selected covariates (Table 5) and the final set of hyperparameters (Table 6).

Table 6Final covariate count (post-decorrelation and RFE) and optimized hyperparameters for each modeled soil property. In instances without case weights, optimal performance was achieved excluding field estimates (silt and sand) or when the property was not estimated in the field (pH, Ntot, Pox, and CEC).

Download Print Version | Download XLSX

2.4 Variable importance

During model calibration, we assessed variable importance using the permutation method for pH, Ntot, Pox, CEC, silt and sand, and the impurity method for clay, BD and SOM. Permutation gives a better estimate of the variable importance than impurity because impurity has a bias towards covariates with more distinct values, making it negatively biased towards categorical covariates as they have a finite number of binary splits due to their limited number of classes (Sandri and Zuccolotto2008, 2010). However, the permutation measure is dependent on the out-of-bag error (Breiman2002). When case weights are high, out-of-bag estimation is not possible because the observations with high weights are selected in the bootstrap sample of all trees, regardless of the sample fraction. Hence, we could not compute the out-of-bag error and use the permutation variable importance measure for these observations in clay, BD and SOM models because they were never out of bag.

2.5 Prediction maps

The calibrated RF and QRF and final set of covariates were used to estimate the mean, median (0.50 quantile; q0.50), 0.05 quantile (q0.05) and 0.95 quantile (q0.95) at every 25 m pixel and each standard depth layer specified by GSM (0–5, 5–15, 15–30, 30–60, 60–100 and 100–200 cm) over the Netherlands. In addition, spatially explicit 90 % prediction interval widths (PI90) were obtained at every 25 m pixel as a measure of prediction uncertainty as follows:

(1) PI90 = q 0.95 - q 0.05 .

We post-processed the mean and median PSF prediction maps to ensure that the three PSF maps summed to 100 %. The predictions of clay, silt and sand were divided by the sum of the three at that location and multiplied by 100 for every 25 m pixel. We chose not to use kriging of the QRF prediction residuals (regression kriging) because there was no spatial autocorrelation in the residuals and to simplify the procedure.

2.6 Accuracy assessment

We evaluated map quality using internal (model-based) and external (model-free) accuracy assessment. At the location and depth, and year in the case of SOM, of a soil property measurement, all quantiles from 0 to 1 at steps of 0.02 were predicted to obtain the PI90 (Eq. 1) as well as the prediction interval coverage probability (PICP) of prediction intervals between 0.02 and 1. The PICP is the proportion of independent observations that fall into the corresponding prediction interval (Papadopoulos et al.2001). We refer to the PICP of the PI90 as the PICP90. The PICP is an indication of how accurately QRF quantifies uncertainty. Prediction uncertainty using PI90 is an example of a model-internal accuracy assessment since it is QRF-dependent, whereas PICP is an external accuracy metric.

Besides PICP, we used two different statistical validation methods for an external accuracy assessment: (1) design-based inference (Brus et al.2011; Brus2022) using either LSK or CCNL laboratory measurements and (2) non-design-based inference using PFB laboratory measurements (Sect. 2.1.2; Table 4). We used the same approach as described in detail in Helfenstein et al. (2022) to adapt design-based inference for statistical validation of prediction maps at different depth layers. However, design-based inference was not used to assess clay, silt, sand and CEC predictions as it was not measured in LSK or CCNL. For non-design-based inference, we used location-grouped 10-fold cross validation of the PFB laboratory measurements, similarly to during model tuning.

To obtain commonly used accuracy metrics, both mean and median predictions were used to calculate residuals. From these residuals, we estimated the mean error (ME or bias), the RMSE and the model efficiency coefficient (MEC):


where n is the number of validation observations; yi and y^i are the ith observation and prediction, respectively, at a certain location, depth and year (for SOM); and y is the mean of all validation observations. Equations (2)–(4) apply for non-design-based inference. The adapted equations for design-based inference are Eqs. (5), (8) and (11) in Helfenstein et al. (2022). We computed these accuracy metrics for all observations and separated into observations pertaining to each depth layer as the latter was necessary for design-based inference (Helfenstein et al.2022).

In addition to rigorous quantitative accuracy assessment, we also evaluated the spatial patterns of BIS-4D prediction maps qualitatively by comparing them to existing soil maps in the Netherlands (de Vries et al.2003; Brus et al.2009; Schoumans and Chardon2015; van den Berg et al.2017; Heinen et al.2022; Knotters et al.2022) and based on expert judgment. We acknowledge that qualitative evaluation was not definitive and is indicative only. Note that we did not compare visual patterns of the national soil map (de Vries et al.2003) and the soil physical units map (BOFEK; Heinen et al.2022) to BIS-4D predictions in peat areas, as covariates of peat classes were used in model calibration (Table 5 and Fig. 5 in Helfenstein et al.2024c).

2.7 BIS-4D updates: pH and SOM

Previous map versions of soil pH in 3D and SOM in 3D+T have recently been published using BIS-4D (Helfenstein et al.2022, 2024c). For soil pH, the version presented here contains several important updates. Firstly, covariates of peat classes (de Vries et al.2003), groundwater classes in agricultural areas (Knotters et al.2018), and Sentinel-2 RGB and NIR bands and spectral indices (Roerink and Mücher2023) were added, all of which were selected and thus used for model calibration and prediction of the updated version (Table 5). We also included decorrelation and RFE to increase the signal-to-noise ratio and make models more parsimonious (Sect. 2.3). For 3D+T maps of SOM, we included the latest national land use map (year 2022) to derive the dynamic 2D+T land use covariates and predict SOM for the year 2023.

2.8 Software and computational framework

The computational framework of BIS-4D is entirely based on open-source software and was operationalized on an Ubuntu 22.04 operating system with 48 cores and 128 GB working memory (RAM). Model input data (soil point data and covariates), scripts and model outputs (BIS-4D soil property prediction maps and their associated uncertainty maps) are openly accessible (Sect. 3).

BIS-4D is mostly based on R (version 4.3.1; R Core Team2023), although GDAL (version 3.7.2; GDAL/OGR contributors2023) and SAGA-GIS (version 7.8.4; Conrad et al.2015) were used during covariate preparation and processing because this massively decreased computation time compared to using similar functions in R. Further details about re-sampling, masking and processing of covariates and reclassification of categorical covariates can be found in Sect. 2.7 of Helfenstein et al. (2022). The indices necessary for the location-grouped 10-fold CV were made using the CAST R package (Meyer2023). The caret package (Kuhn2008, 2019, 2022) was used for the tuning and selection of hyperparameters. We used the ranger package (Wright and Ziegler2017) with the option quantreg to grow a QRF during calibration and without it to grow a RF during RFE and tuning. For predictions, the option quantiles was used to predict quantiles, while the option response was used to predict the mean. A combination of the ranger and terra packages was used for predicting at all locations and depths. We used QGIS (version 3.32.3; QGIS Development Team2023) and the rasterVis (Lamigueiro and Hijmans2023) and mapview (Appelhans et al.2023) R packages for exploratory and qualitative analysis and visualization of covariates and prediction maps. The computational workflow for all BIS-4D maps took approximately 5700 CPU hours.

3 Code and data availability

The BIS-4D soil property prediction maps at 25 m resolution can be downloaded at (Helfenstein et al.2024a). Prediction maps of the mean, median, 0.05 and 0.95 quantiles and the PI90 are available for each standard depth layer specified by GSM (0–5, 5–15, 15–30, 30–60, 60–100 and 100–200 cm). For SOM, maps at the same resolution and for the same depth layers are available for the years 1953, 1960, 1970, 1980, 1990, 2000, 2010, 2020 and 2023.

Regarding BIS-4D model inputs, the soil point data on laboratory measurements and field estimates used during model calibration (PFB and BPK data) are publicly available at (Helfenstein et al.2024d). The georeferenced soil point data of PFB and BPK can also be viewed at (last access: 23 January 2024). LSK and CCNL data used for design-based inference are not open due to privacy agreements. The pre-processed covariates that are openly available can be downloaded at 25 m resolution at (Helfenstein et al.2024b). This includes the majority of the covariates used for BIS-4D, with the main exception being the covariates related to the national forestry inventory, since these data are closed. The coordinate reference system of the spatial point and raster data is EPSG:28992 (Amersfoort/RD New).

A frozen version of the BIS-4D code is available at (Helfenstein2024b). The GitLab code repository is complete with the exception of BIS database credentials and the LSK and CCNL data. All data and code are available under the CC BY 4.0 license, except for the covariates (Helfenstein et al.2024b), which are available under the CC BY-NC-SA 4.0 license.

4 Results and discussion

BIS-4D prediction maps for every GSM depth layer at 25 m resolution can be downloaded at (Helfenstein et al.2024a). These include predictions of the mean, 0.05, 0.50 (median) and 0.95 quantiles and the PI90 of clay, silt, sand, BD, pH, Ntot, Pox and CEC. For SOM, these prediction maps are available for the years 1953, 1960, 1970, 1980, 1990, 2000, 2010, 2020 and 2023 (Sect. 3). An overview of all prediction maps together with the associated accuracy metrics (ME, RMSE, MEC and PICP) and variable important information can be found in the Supplement, which is organized by target soil property.

4.1 Accuracy assessment

4.1.1 Quantitative accuracy assessment

The accuracy of the produced maps varied considerably depending on the soil property (Tables 7 and  8; see the Supplement). Based on 10-fold cross validation (Table 7), the accuracy of mean predictions over all depths for clay, sand, BD, pH and Ntot maps was the highest (MEC > 0.70), followed by SOM and silt (MEC > 0.60). Mean predictions for Pox and CEC were the least accurate (MEC of 0.54 and 0.49, respectively). Design-based inference separated by the depth layer confirms the high accuracy of pH prediction maps (Table 8). MEC values computed for mean and median predictions using design-based inference were lower for BD (0.34–0.78) and Ntot (0.27–0.52) than when using 10-fold cross validation. Mean and median Pox maps were very inaccurate (MEC of −0.11–0.38) based on design-based inference. The large differences in accuracy between 10-fold cross validation using PFB laboratory measurements and design-based inference using LSK laboratory measurements for BD and Pox may be due to the clustered and limited spatial distribution of calibration data for those soil properties (Fig. 2d and h). Therefore, for BD and Pox, metrics using 10-fold cross validation are likely overly optimistic.

Table 7Accuracy metrics of BIS-4D soil property maps using mean and median predictions, computed using 10-fold cross validation (Sect. 2.6). Units of ME and RMSE are in units of the measured soil property (Table 1).

Download Print Version | Download XLSX

Table 8MEC for mean and median predictions of BIS-4D soil property maps, separated by the depth layer and computed using either 10-fold cross validation (CV) of PFB laboratory measurements or design-based inference (DBI) using LSK or CCNL data (Table 4). DBI for Ntot at 100–200 cm depth was not possible because soil samples were not collected below 100 cm in CCNL (Sect. 2.1.2 and 2.6). However, for this depth layer, CV metrics are included in the Supplement (Table S7).

Download Print Version | Download XLSX

The RMSE and ME were low for most soil properties (Table 7 and the Supplement). The RMSE of sand was higher than for clay and silt, even though the MEC of sand indicates higher model performance for sand than for silt. This can be explained by the high proportion of regions in the Netherlands with very high sand content (>75 %), i.e., the Pleistocene sandy areas shown in pink in Fig. 4d and h. In comparison, laboratory measurements of clay and silt content were rarely >75 % (Fig. 3).

Figure 4Mean predicted clay [%] (a, e), silt [%] (b, f) and sand [%] content (c, g) at 60–100 cm depth and associated prediction uncertainty (PI90: 90th prediction interval width) and the soil physical units map of the Netherlands (BOFEK; Heinen et al.2022; d, h) in comparison. The soil physical unit codes can be found in Heinen et al. (2022); here grouped into the main categories (code beginning with 1 representing peat, with 2 peaty, with 3 sand, with 4 loam/clay and with 5 loess). The zoomed-in area around Wageningen was chosen since this area contains all main soil physical categories except loess.

The differences in accuracy between mean and median prediction maps varied slightly between soil properties. The low and high values of the mean predictions were systematically biased so that the low values of the observed soil property were overestimated and the high values underestimated, both to varying degrees for different target soil properties (Figs. S10, S21, S32, S43, S54, S65, S76, S87 and S98). Thus, while the mean predictions were slightly less biased than the median predictions when averaging over all values except for SOM (Table 7), they were more biased than the median predictions for the low and high values. For soil properties where calibration data were positively skewed (Fig. 3), i.e., all soil properties except sand, BD and pH, the bias of mean predictions was negative, whereas the bias of median predictions was positive (Table 7). In contrast to the findings based on 10-fold cross validation, design-based inference of Ntot revealed that median predictions were less biased (between −609 and 120 mg kg−1; see the Supplement) than mean predictions (between −511 and −1408 mg kg−1; see the Supplement). Higher accuracy of median predicted Ntot was also reflected in lower RMSE (Table S7) and higher MEC values (Table 8). In summary, although it depends on the use, we overall recommend to use median predictions since low and high values were less biased and ME, RMSE and MEC values for both mean and median predictions were similar.

Mean predictions are more sensitive to extreme values and outliers than median predictions. For instance, in mineral soils, the predicted conditional distribution of SOM, Ntot, Pox and CEC was positively skewed and median predictions were usually smaller than mean predictions (e.g., von Hippel2005, Fig. 1 therein). In peat soils, the opposite was the case. Here, the predicted conditional distribution of SOM, Ntot, Pox and CEC was negatively skewed and median predictions were larger than mean predictions. For these soil properties, mean predictions were thus systematically higher than median predictions in mineral soils, whereas mean predictions were systematically lower than median predictions in peat soils.

The maps of prediction uncertainty (PI90) for every GSM depth layer revealed that uncertainty was high when mean and median predictions fell within a range with limited calibration data (Supplement). This meant that for most soil properties, uncertainty was high in areas where predictions were high due to the positively skewed distribution of observation data (Fig. 3). For example, the positive correlation between increasing uncertainty with increasing predictions can be clearly observed for clay and silt in Fig. 4e, f, i and j. The same positive correlation between predictions and uncertainty was observed for Ntot over depth (Fig. 5e and f). We found a similar pattern of high uncertainty in peatlands due to high predictions in these areas for SOM, Pox and CEC. However, given its bimodal distribution, the uncertainty for sand was the highest in areas where predictions ranged between 25 % and 75 % (for example, in the river areas) and uncertainty was comparatively low in marine clay areas (<25 % sand) and Pleistocene areas (>75 % sand) (Fig. 4c, g and k).

Prediction uncertainty for most soil properties increased with increasing depth (e.g., Fig. 5f), except if mean and median predictions decreased substantially over depth, as was the case for Pox (Figs. S78–S85 in the Supplement). Higher uncertainty at lower depths is in line with worse accuracy metrics at lower depths (Table 8 in the Supplement) and this tendency was found in the majority of recently reviewed DSM studies (Chen et al.2022). Deeper soil layers are generally more difficult to predict because limited information about the subsoil can be derived from most covariates, and especially remote sensing products. However, for BD and pH, the accuracy from 15–30 cm depth may have been higher than from 0–15 cm depth because only 245 observations were available for statistical validation in LSK from 15–30 cm depth (Tables S4 and S6). Therefore, the metrics computed via design-based inference from 15–30 cm depth for BD and pH are likely less representative of map quality compared to metrics of the other depth layers, where many more observations were available.

Prediction uncertainty in most soil properties was also higher in urban areas, which can be attributed to limited soil samples and heavily disturbed soils in urban areas. With increasing population growth in an already densely populated country, this highlights the need to map urban soils (Römkens and Oenema2004; Vasenev et al.2014, 2021; Kortleve et al.2023).

Figure 5Median predicted BD [g cm−3(a), Ntot [mg kg−1(b), Pox [mmol kg−1(c), and CEC [mmol(c) kg−1(d) at 0–5 cm depth; and median predicted Ntot (e) and PI90 (90th prediction interval width) as a measure of the associated prediction uncertainty (f) along the depth transect shown in (b).

The PICP90 (Table 7) and the PICP (the Supplement) indicated that prediction uncertainty was estimated relatively accurately using QRF, but small differences were found among the predicted soil properties. For clay content, the PICP90 was between 0.82 and 0.86 (Table S1) and hence less than 0.90, indicating that the uncertainty in clay predictions was underestimated. The uncertainty in BD based on PFB laboratory measurements was slightly underestimated (0.86; Table 7) but was slightly overestimated based on LSK laboratory measurements (0.88–0.95; Table S4). For SOM, the PICP90 varied strongly with depth (0.75–0.96; Table S5), but the PICP was overall very accurate for all depths combined (Fig. S53). In our study, the soil properties for which field estimates were included during calibration were the only ones for which the PI90 was sometimes underestimated. Similarly, Chen et al. (2023) found that increasing the proportion of spectral estimates combined with conventional laboratory measurements decreased the PI90. Hence, if calibration data are a smoothed version of the truth, which may be the case with predictions of spectral models and field estimates, this tends to lead to underestimation of the “true” uncertainty. The aim of sharp, i.e., narrow conditional probability distributions by including various types of observational data is desirable only if ensuring that the uncertainty is still reliable, e.g., by computing the PICP (Schmidinger and Heuvelink2023). This is important to avoid presenting overoptimistic results to end users. Besides clay, BD and SOM, prediction uncertainty for the remaining target soil properties was accurate but marginally overestimated (0.89–0.97) based on the independent datasets used for statistical validation (LSK and CCNL; see the Supplement). Hence, the PICP indicates that silt, sand, pH, Ntot, Pox and CEC maps are somewhat more accurate than suggested by the prediction uncertainty (PI90).

4.1.2 Qualitative accuracy assessment

The BIS-4D maps of the nine predicted soil properties align with the national soil map of the Netherlands (de Vries et al.2003). This can be seen when comparing our maps (Figs. 4 and 5; see the Supplement) with the soil physical units map (BOFEK; Fig. 4d and h; Heinen et al.2022) derived from the national soil map of the Netherlands. Further information on the comparison of BIS-4D maps to previous soil pH maps (Brus et al.2009) can be found in Sect. 4.2 of Helfenstein et al. (2022) and to previous SOM maps (Brus et al.2009; van den Berg et al.2017; Knotters et al.2022) in Helfenstein et al. (2024c). Nonetheless, visual evaluation of the maps also revealed several limitations.

The maps of soil texture or particle size fractions (clay, silt and sand) of the mineral soil component should be used with caution in peatlands (approximately 15 % of the surface area) since natural peat only consists of organic matter without a mineral component. However, the low-lying fen peatlands, located mostly in the west and northwest of the Netherlands, typically also contain some clay, silt or sometimes even sand due to past flooding events (Edelmann1950; de Bakker and Schelling1966, 1989; Brouwer et al.2023). Drained organic soils, particularly when under agricultural use, can also contain mineral components introduced or mixed in from mineral soil horizons from below or above the organic soil horizon. Nonetheless, 30 % clay content in a soil composed mostly of peat in absolute terms contains less clay than a mineral soil with 30 % clay content.

Visual examination of the BIS-4D maps reveals artifacts from the covariates. Although water and buildings were cropped out, some mapping artifacts remained, such as small buildings, roads and railways. For instance, the road on top of the dike, parallel to and south of the Rhine River is clearly visible in Fig. 4e–k. This highlights the difficulty of spatial modeling approaches such as DSM that rely strongly on remote sensing products. Other artifacts were due to the combination of several Sentinel-2 images from different days in 1 month to obtain one monthly, cloud-free mosaic (Sect. 2.2). Image mosaicking created artificial lines due to alterations in the brightness, hue and colors from images of different days. In addition, a few 25 m pixels contain no data even though they were not water or buildings. This may be due to no data values in some covariates but should be explored further in an updated version. Finally, there were also orthogonal artifacts, most likely due to using northing and easting coordinates as covariates, which can be largely removed by also including oblique axes in many additional directions (Møller et al.2020).

4.2 Strengths

BIS-4D maps fill the missing data gap of spatial soil property information on a national scale in the Netherlands and bring substantial improvements to previously mapped soil properties. The main strengths of BIS-4D are (1) the ability to provide information of soil properties as opposed to soil types; (2) the high spatial resolution (25 m); (3) the accuracy and uncertainty assessment based on best practices; (4) the benefits of machine learning combined with large amounts of data; (5) the flexibility to predict in 3D and 3D+T; and (6) the open availability of model code and data, making BIS-4D fully reproducible and easy to update.

The BIS-4D maps have several advantages compared to previous soil maps of the Netherlands. While categorical maps of soil type (de Vries et al.2003) and derived thematic maps (Brouwer and van der Werff2012; Brouwer et al.2018; van Delft and Maas2022, 2023; Heinen et al.2022) are important and useful, many users require information on specific numerical soil properties (Sect. 1). We acknowledge that clay content (Brus et al.2009), SOM (Brus et al.2009; van den Berg et al.2017; Knotters et al.2022), pH (Brus et al.2009), and soil properties related to soil texture (Heinen et al.2022) and Pox (Schoumans and Chardon2015) have previously been mapped on a national scale in the Netherlands. However, these maps were at much coarser resolution, accuracy was either not assessed or not assessed using design-based statistical inference, quantification and evaluation of uncertainty were missing, mapping approaches did not include machine learning and used only a few covariates, and predictions for one or several depth layers were modeled separately and only Knotters et al. (2022) assessed changes over time. The only standard GSM soil properties that we did not map are SOC, plant exploitable (effective) depth, depth to rock and coarse fragments (Arrouays et al.2014a, b, 2015). We mapped SOM instead of SOC because, in the Netherlands, SOC was not included in routine soil analyses until recent years. However, SOC can be derived from SOM, as investigated in other studies in the Netherlands (van Tol-Leenders et al.2019; van den Elsen et al.2020; Teuling et al.2021; Knotters et al.2022). Plant exploitable (effective) depth is mostly limited by high groundwater levels in most regions of the country. Since groundwater levels have been extensively mapped in the Netherlands (de Gruijter et al.2004; Hoogland et al.2014; Knotters et al.2018), mapping plant exploitable (effective) depth was not deemed necessary. Depth to rock and coarse fragments are not relevant on a national scale in the Netherlands, as the substrate materials of Dutch soils are almost exclusively either Pleistocene sand, fine-grained Quaternary sediments or peat.

Another strength of BIS-4D is that maps are at a high spatial resolution of 25 m. As covariates such as remote sensing products and national maps of land use (Hazeu et al.2023) and digital elevation models (AHN2023) are nowadays available at 5–25 m resolution, useful information for modeling complex relationships between soil-forming factors such as land cover and topography and soil properties is provided on these scales. The increasing availability of high-resolution information in soil-related domains has also increased the demand for high-resolution soil maps. While high-resolution products such as BIS-4D bring many advantages, it is crucial to emphasize that resolution is not an indicator of accuracy and should not be used solely to determine a map's fitness for use (de Bruin et al.2001; Malone et al.2013; Knotters and Walvoort2020; Szatmári et al.2021).

One of the main advantages of BIS-4D is the rigorous map quality evaluation using design-based statistical inference and prediction uncertainty. Based on sampling theory (Cochran1977; de Gruijter et al.2006; Gregoire and Valentine2007), map accuracy should be assessed with design-based statistical inference using a probability sample whenever possible, as this provides a better estimate of the true map accuracy compared to non-design-based approaches (Brus et al.2011). Moreover, it also produces confidence intervals (Tables S4–S8) so that we know how close the estimate of the map accuracy is to the true map accuracy. We were able to use design-based inference for BD, pH, SOM, Ntot and Pox maps due to the availability of the LSK and CCNL datasets. We are not aware of any other GSM products that used design-based inference to evaluate map accuracy on a national scale. For soil properties for which design-based inference was not possible, i.e., for clay, silt, sand and CEC, we used location-grouped 10-fold cross validation, as recommended in the case of non-clustered data (Wadoux et al.2021a; de Bruin et al.2022). In addition, BIS-4D maps provide spatially explicit estimations of prediction uncertainty (PI90), including GSM accuracy thresholds for soil pH (Helfenstein et al.2022), and we evaluated the accuracy of the uncertainty using PICP.

Another strength of BIS-4D, for example, when compared to previous soil property maps in the Netherlands (e.g., Brus and Heuvelink2007; Brus et al.2009; van den Berg et al.2017), is that machine learning leads to more accurate predictions than other geostatistical and regression techniques. Ensemble decision tree models such as RF and QRF have repeatedly outperformed other spatial interpolation methods (e.g., Hengl et al.2015; Nussbaum et al.2017; Keskin et al.2019; Khaledian and Miller2020). Ensemble decision tree models are able to capture complex, non-linear relationships between the covariates and soil properties and are widely used in recent DSM studies (Vaysse and Lagacherie2017; Heuvelink et al.2020; Poggio et al.2021; Baltensweiler et al.2021; Nussbaum et al.2023).

BIS-4D maps for clay, silt, sand, BD, pH, Ntot, Pox and CEC are in 3D (between 0 and 2 m depth) and, for SOM, they are also dynamic (3D+T; see the Supplement and Helfenstein et al.2024c). This fills a largely missing gap of soil information in deeper layers (Chen et al.2022). In addition, BIS-4D can predict at any depth as opposed to recalibrating models when mapping individual depth layers separately (Ma et al.2021). This improves model flexibility and efficiency and a larger amount of data can be leveraged during model tuning and calibration. For example, routine agronomic soil sampling depths in the Netherlands are 0–10 cm for grasslands and 0–25 cm for croplands, thereby deviating from the GSM standard depths (Arrouays et al.2014a, b, 2015). Predictions and associated uncertainty for those depths can be provided using BIS-4D without recalibrating models. This is particularly useful for uncertainty, which, unlike mean and median predictions, cannot be aggregated using, e.g., weighted averaging over depth layers (Sect. 4.1.1). In addition, we developed innovative covariates explicit in 3D+T, presenting a novel opportunity to extend the predictive power of machine learning to 3D+T (Helfenstein et al.2024c). This provided a new opportunity for monitoring SOM-related soil health using a method that is explicit in 3D space.

Lastly, compared to the time-consuming effort of updating conventional soil maps, DSM products such as BIS-4D can be easily extended to other soil properties in BIS and can be updated and delivered on demand (Heuvelink et al.2010; Kempen et al.2009, 2012b, 2015). In comparison to an earlier version for soil pH (Helfenstein et al.2022), the number of covariates has been substantially decreased during model selection (Sect. 2.3), which benefits reproducibility and possibilities of updating maps. The model code, workflow, inputs and outputs are well documented and openly available, making procedures reproducible and easy to update (Sect. 3).

4.3 Limitations and improvements

Uncertainty in DSM products such as BIS-4D can be linked to three overarching sources: (1) the quantity and quality of soil point data, (2) the quantity and quality of covariates, and (3) the model structure (Heuvelink2014, 2018). In the following, we discuss the limitations of BIS-4D maps with regard to these three sources of uncertainty and suggest improvements.

4.3.1 Soil point data

Measurement errors and differences in measurement methods of the soil point data may have contributed to the uncertainty in BIS-4D maps. For example, Fe, Al and P extracted by oxalate extraction are considered to consist of amorphous iron and aluminum (hydr)oxides and P bound to those oxides. However, a fraction of oxalate-extractable P in peat soils likely consist of P bound to organically complexed Fe and Al since those are also partially extracted during the oxalate extraction (McKeague1967; McKeague et al.1971; van der Zee et al.1990; Schoumans2013; Schoumans and Chardon2015). Recent research has devised methods to quantify the uncertainty in soil laboratory measurements (van Leeuwen et al.2021) and to incorporate these errors into machine learning algorithms (van der Westhuizen et al.2022). Furthermore, several slightly different methods, standards and laboratory facilities were used to measure Ntot, Pox and CEC (Maring et al.2009, Appendix E). This introduced uncertainty that can be minimized by standardizing laboratory measurements and procedures. Positional uncertainty due to marking locations on a 1:25 000 topographic map most likely also contributed to overall uncertainty in the BIS-4D maps, as investigated in other studies (Carré et al.2007; Grimm and Behrens2010).

There were several limitations related to the spatial and spatio-temporal distribution of the soil point data used in BIS-4D. The calibration data of BD, Pox and, to a lesser extent, CEC were spatially clustered (Fig. 2), which most likely affected mapping accuracy of those soil properties (Sect. 4.1). In addition, no wet-chemical laboratory measurements were available as part of a probability sample (LSK and CCNL) for design-based statistical inference of clay, silt, sand and CEC prediction maps (Sect. 2.1.2). As most of the soil point data were collected between 1950 and 2000, soil measurement age and time should be addressed also for other soil properties besides SOM (Arrouays et al.2017). Ntot and CEC are strongly linked to SOM and thus temporal changes may be similar to mapped SOM changes (Helfenstein et al.2024c). BD, pH, Ntot, Pox and CEC likely changed due to land use and management. However, yearly variation in Pox is relatively small since P binds strongly to soil particles, and the plant available fractions of P with short turnover times are less than 15 % of the total reversibly bound P pool (Withers et al.2014, Fig. 3), which is what is measured with Pox (Lookman et al.1995; Neyroud and Lischer2003). Large quantities of topsoil data are collected for agronomic surveys every 4 years in the Netherlands (BZK2022; Eurofins Agro2024a, b), but only a small part of these are not privacy-protected, making it challenging to incorporate in DSM approaches. Although the point data suggest good spatial coverage of most of the basic soil properties in the Netherlands, there is a major lack of repeated laboratory measurements collected using identical sampling strategies over time, as discussed in Knotters et al. (2022) and Helfenstein et al. (2024c). A consistent national soil monitoring scheme would be beneficial for modelling dynamic soil properties in 3D+T, updating static BIS-4D maps and accuracy assessment with more recent data.

4.3.2 Covariates

Although BIS-4D was able to make use of a large range of high-quality, country-specific covariates (Table 5), the main variable missing in our modelling approach is more detailed land management data, which is a common challenge in DSM (Finke2012; Arrouays et al.2021). Land cover and land use covariates only indirectly provided information on land management. From 2005 onwards, annual data on the specific crop type for every agricultural parcel in the Netherlands were available (BRP Gewaspercelen; EZK2019), but these were never selected among the final covariates used for model calibration (and therefore are not shown in Table 5). This implies that they did not provide additional useful information for the spatial distribution of the target soil properties at the national scale, although different drivers may be relevant on regional and local scales (Sect. 4.4). However, regardless of the scale, national crop parcel data do not capture information on management decisions such as fertilizer inputs, liming and ploughing frequency on agricultural lands, and maintaining forests and nature areas. These management decisions are highly relevant for many of the mapped soil properties, with the exception of particle size fractions. For example, BD is strongly dependent on the size and driving frequency of agricultural machinery on the fields (Stettler et al.2014).

As another example, Pox exhibits considerable small-scale spatial variability, as discussed and made evident by the high nugget in the semivariogram in Fig. 6 of Lookman et al. (1995). As P in the form of phosphate is bound in the soil much stronger than N or other plant nutrients affected by the base cation saturation and CEC, there are large legacy effects due to historic management not captured in the covariates currently used in BIS-4D. In our study, the three most important covariates for modelling Pox were the covariates related to soil horizon sampling depth (Fig. S88). The relationship of Pox with soil depth is supported by empirical findings of the maximum P sorption capacity decreasing with soil depth, especially in sandy soils. P from fertilization largely stays in the upper soil layers. Moreover, given that Pox map quality was poor (Table 8 and the Supplement), the relative importance of depth suggests that the other covariates did not explain the spatial variation in Pox well, likely due to missing (historic) management data. This is supported by other studies of mapping soil P content at high resolution (Delmas et al.2015; Matos-Moreira et al.2017; Kull et al.2023). Although this does not solve the problem of missing management data, one easy step towards improving the accuracy of BD and Pox and other management-dependent soil properties is to only map them for agricultural areas, as was done in the Netherlands for amorphous iron and aluminum (hydr)oxides (van Doorn et al.2024). We expect that including dynamic covariates of land management and climate, as discussed in Helfenstein et al. (2024c), would likely also improve modelling dynamic soil properties in 3D+T.

4.3.3 Model structure

Despite the many advantages of using QRF for DSM (Sect. 4.2), predictions may be further improved using methods such as convolutional or recursive neural networks (deep learning; Behrens et al.2005, 2018a; Padarian et al.2019b; Wadoux2019; Wadoux et al.2019) or transfer learning (Liu et al.2018; Padarian et al.2019a; Seidel et al.2019; Helfenstein et al.2021; Baumann et al.2021), defined as the process of sharing intra-domain information and rules learned by general models to a local domain (Pan and Yang2010). We recommend future research to investigate the use of deep learning and transfer learning in the Netherlands for SOM due to the large amount of SOM data and more opportunities in accounting for differences in observational quality (field estimates and laboratory measurements) using more complex models. However, to the best of our knowledge, deep learning has only outperformed ensemble decision tree models when using a small number of covariates covering only some of the soil-forming factors, from which hyper-covariates are then derived (Wadoux2019). Hence, deep learning may not improve predictions in the Netherlands, where large amounts of high-quality covariates are readily available for all soil-forming factors. In addition, quantifying model-based uncertainty using deep learning remains a challenge. Although model-free approaches of estimating uncertainty using deep learning, e.g., involving bootstrapping (Padarian et al.2019b; Wadoux2019), have been used, we are not aware of studies that have compared the accuracy of these uncertainty estimations to QRF-based uncertainty (PI90).

One of the main limitations of the BIS-4D modelling approach is that QRF predictions cannot be used to compute the uncertainty in spatial aggregates, for example, when aggregating prediction maps of different depth layers or computing average values of a soil property for a specific land use or province. This requires quantifying cross and spatial correlation in prediction errors, which can be accounted for by taking a multivariate or geostatistical approach (Szatmári et al.2021; van der Westhuizen et al.2022; Wadoux and Heuvelink2023).

4.4 Assessment scale

We recommend using BIS-4D maps on a national scale as long as the map quality based on the provided accuracy metrics (Tables 7 and 8; see the Supplement) and prediction uncertainty (Figs. 4i–k and 5f; see the Supplement) is sufficient for the intended use. The model was developed for the national scale for multiple land uses. Foremost, BIS-4D maps contribute to the GSM project by delivering high-resolution 3D (and 3D+T for SOM) maps of key soil properties with quantified uncertainty according to GSM specifications for the Netherlands. The BIS-4D maps may be especially useful for initiatives that require spatially explicit soil information across all land uses and soil types of the Netherlands. This may include national contributions to United Nations and pan-European directives and policies (Panagos et al.2022), such as the Green Deal, Common Agricultural Policy, Zero Pollution Action Plan, EU Soil Strategy for 2030, Soil Deal (European Commission2021) and Proposal for a Directive on Soil Monitoring and Resilience (European Commission2023). For example, clay, silt, sand and SOM maps can be used to improve estimates of soil-derived greenhouse gas emissions from the LULUCF sector for the Netherlands (Arets et al.2020).

Many potential users of BIS-4D may require information specifically for one land use or soil type. Perhaps most commonly, users may need information for agricultural soils. For example, maps of clay, silt, sand and SOM content can provide information used to estimate the carbon sequestration potential for the Smart Land Use project (Slier et al.2023), which is focused specifically on mineral soils under agricultural use on a national scale. Policymakers are mostly interested in more complex soil information such as soil health, soil functions or soil-based ecosystem services. Complex soil information can be quantified more accurately using a combination of basic soil properties, such as provided using BIS-4D, employing methods like pedotransfer functions and existing tools. For agricultural soils in the Netherlands, such tools include OSI (Ros et al.2022; Ros2023) and BLN 2.0 (Ros et al.2023). Although pH, Ntot, Pox and CEC can be used as approximations on a national scale, pH, plant-available N and P, and CEC are part of routine agronomic soil analyses. Therefore, maps of these soil properties may be more useful for forests and nature areas, where the base cation occupation of the CEC and pH should generally remain high enough to prevent Al toxicity. For soil pH, BIS-4D predictions should be compared to predictions in Wamelink et al. (2019), who mapped soil pH in Dutch nature areas. As soil variability is linked to soil type (e.g., mineral vs. organic) and land use, we expect that BIS-4D model predictions would improve when modeling only one land use or soil type separately. However, this was not the scale of assessment aimed for with BIS-4D.

BIS-4D maps may also be used on a regional scale as long as the accuracy allows for it and no better product is available. By regional scale we mean the level of provinces; regional water authorities, which are typically composed of one or more polders or watersheds; or large municipalities. These recommendations hold true for clay and sand content and pH in particular, which were predicted with higher accuracy than the other soil properties. However, regional management decisions come with social and economical risks. The costs of poor management decisions due to the use of inaccurate or soil information that is not detailed enough are often several magnitudes larger than investments for conducting a more detailed regional soil survey (Knotters and Vroon2015; Keller et al.2018). In agreement with Chen et al. (2022), more research is necessary in relating DSM performance indicators such as uncertainty to cost–benefit and risk assessment analysis for improving decision support. We do not recommend the use of BIS-4D maps on a farm or field scale, as the uncertainty in predictions is most likely too high for the precision required by the farmer. Drivers of soil variation vary locally and were presumably not captured on this scale by the soil point data, covariates and model structure. As shown in Helfenstein et al. (2022), even for a soil property like soil pH, which was relatively easy to predict, less than 10 % of the map pixels were designated with one of the highest two GSM accuracy thresholds (AA and AAA). On such a local scale, we expect that the time and costs invested in a new soil survey outweigh the risks of using inaccurate soil data (Lemercier et al.2022).

4.5 Good practices for proper use

Based on the accuracy assessment of BIS-4D maps (Sect. 4.1.1), clay, sand and pH maps were the most accurate. This is mostly in agreement with Chen et al. (2022), which found that pH was the best predicted standard GSM soil property, and PSFs (i.e., clay, silt and sand) were predicted to be third-best, based on a review of 244 articles. BIS-4D map quality of silt, BD, SOM, Ntot and CEC was lower. For Pox, we only recommend the produced maps for a qualitative overview of Pox spatial distribution in the Netherlands. We have summarized the following simple chronological steps for users to help decide whether BIS-4D maps may be suitable for their intended purpose:

  1. Choose one or multiple soil properties of interest.

  2. Choose the depth layer(s) (0–5, 5–15, 15–30, 30–60, 60–100 and 100–200 cm) and, in the case of SOM, the year (1953, 1960, 1970, 1980, 1990, 2000, 2010, 2020, 2023), for which soil information is needed.

  3. Consult the accuracy metrics of mean and median predictions of the soil property and depth layer of interest (Tables 7 and 8; see the Supplement), keeping in mind that the accuracy in the intended area of use may differ from the overall accuracy of the map. If the overall map quality based on these accuracy metrics is within an acceptable range for the intended purpose, continue to step 4. Use accuracy metrics based on design-based statistical inference for the soil properties for which it is available, i.e., BD, pH, SOM, Ntot and Pox (Table 8; see the Supplement).

  4. Choose whether to use mean or median prediction maps by comparing accuracy metrics of mean and median predictions (Tables 7 and 8; see the Supplement). Consult Sect. 4.1.1 for differences in mean and median predictions found using BIS-4D.

  5. Download the mean and/or median prediction maps for the chosen soil property and depth layer as well as maps of the associated uncertainty (0.05 quantile, 0.95 quantile and/or PI90) and open them using GIS software. If soil information is required for a specific area, continue to step 6.

  6. Prediction uncertainty is only useful for end users if it is reliable (Schmidinger and Heuvelink2023). Therefore, check whether the prediction uncertainty is reliable by consulting the PICP. If the PICP90 is close to 0.90 (Table 7) and the PICP plot close to the 1:1 line (Supplement), then the provided prediction uncertainty map is reliable.

  7. Ideally, prediction uncertainty should also be sharp; i.e., the PI90 should be as narrow as possible (Schmidinger and Heuvelink2023). Decide whether the PI90 is within an acceptable range for the intended purpose. If possible, fitness for use can be determined by analyzing how uncertainties in BIS-4D maps propagate through the intended usage, for example, for an environmental model that uses BIS-4D maps as input. Commonly used uncertainty propagation methods include the Taylor series or Monte Carlo methods (Heuvelink2018).

Video supplement

A short research pitch of BIS-4D is available at (Helfenstein2024a).


The supplement related to this article is available online at:

Author contributions

All authors contributed to the scientific research and writing of this paper. AH developed BIS-4D, prepared the data, and wrote all R scripts and sections of the paper. VLM, GBMH and MJDHB provided valuable inputs and suggestions during all stages of the research, including conceptualization of the modeling framework, interpretation of results and improving the writing. MvD, KT and DJJW provided expert knowledge on various mapped soil properties within the Dutch context.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


We would like to thank Joop P. Okx for envisioning and realizing the need for spatially explicit soil property information in the Netherlands and his involvement in the early stages of this project. We thank Fokke de Brouwer for his expertise on soil surveying and mapping in the Netherlands. Furthermore, we thank Dorothee van Tol-Leenders and WENR colleagues for their ongoing efforts to ultimately integrate BIS-4D maps at (last access: 22 January 2024). We express our utmost gratitude to all data providers. This includes decades of rigorous work by soil surveyors and scientists from Stichting Bodemkartering (StiBoKa); DLO Winand Staring Centre for Integrated Land, Soil and Water Research, later known as Alterra; and Wageningen Environmental Research. And, lastly, this includes the providers of data used to derive hundreds of different covariates for this study, made available from different sources and by various institutions.

Financial support

This research has been supported by the Dutch Ministry of Agriculture, Nature and Food Quality (grant nos. WOT-04-013-010 and KB-36-001-014).

Review statement

This paper was edited by Conrad Jackisch and reviewed by David Rossiter and two anonymous referees.


Aarts, N. and Leeuwis, C.: The Politics of Changing the Dutch Agri-Food System, Journal of Political Sociology, 1, 1,, 2023. a

Adhikari, K., Kheir, R. B., Greve, M. B., Bøcher, P. K., Malone, B. P., Minasny, B., McBratney, A. B., and Greve, M. H.: High-Resolution 3-D Mapping of Soil Texture in Denmark, Soil Sci. Soc. Am. J., 77, 860–876,, 2013. a

AHN: Actueel Hoogtebestand Nederland (AHN), AHN, (last access: 23 January 2024), 2023. a, b

Aitchison, J.: The statistical analysis of compositional data, Chapman and Hall, London, 1986. a

Akpa, S. I. C., Odeh, I. O. A., Bishop, T. F. A., and Hartemink, A. E.: Digital Mapping of Soil Particle-Size Fractions for Nigeria, Soil Sci. Soc. Am. J., 78, 1953–1966,, 2014. a

Alterra: Historisch Grondgebruik Nederland (HGN), (last access: 15 January 2024), 2004. a, b

Amirian-Chakan, A., Minasny, B., Taghizadeh-Mehrjardi, R., Akbarifazli, R., Darvishpasand, Z., and Khordehbin, S.: Some practical aspects of predicting texture data in digital soil mapping, Soil Till. Res., 194, 104289,, 2019. a, b

Appelhans, T., Detsch, F., Reudenbach, C., and Woellauer, S.: mapview: Interactive viewing of spatial data in R, Tech. rep., (last access: 15 January 2024), 2023. a

Arets, E. J. M. M., Kolk, J. W. H. v. d., Hengeveld, G. M., Lesschen, J. P., Kramer, H., Kuikman, P. J., and Schelhaas, N. J.: Greenhouse gas reporting for the LULUCF sector in the Netherlands: methodological background, update 2020, WOt-technical report 168, Statutory Research Tasks Unit for Nature & the Environment (WOT Natuur & Milieu), Wageningen, the Netherlands, WOT Natuur & Milieu,, 2020. a, b

Arrouays, D., Grundy, M. G., Hartemink, A. E., Hempel, J. W., Heuvelink, G. B. M., Hong, S. Y., Lagacherie, P., Lelyk, G., McBratney, A. B., McKenzie, N. J., Mendonca-Santos, M. d. L., Minasny, B., Montanarella, L., Odeh, I. O. A., Sanchez, P. A., Thompson, J. A., and Zhang, G.-L.: Chapter Three – GlobalSoilMap: Toward a Fine-Resolution Global Grid of Soil Properties, in: Advances in Agronomy, edited by: Sparks, D. L., vol. 125, 93–134, Academic Press,, 2014a. a, b, c, d

Arrouays, D., McKenzie, N., Hempel, J., Forges, A. R. d., and McBratney, A. (Eds.): GlobalSoilMap: Basis of the global spatial soil information system, CRC Press Taylor & Francis Group, Boca Raton, ISBN 978-1-315-77558-6, 2014b. a, b, c, d

Arrouays, D., McBratney, A., Minasny, B., Hempel, J., Heuvelink, G. B. M., MacMillan, R. A., Hartemink, A., Lagacherie, P., and McKenzie, N.: The GlobalSoilMap project specifications, in: Proceedings of the 1st GlobalSoilMap Conference, 9–12,, 2015. a, b, c, d

Arrouays, D., Leenaars, J. G. B., Richer-de Forges, A. C., Adhikari, K., Ballabio, C., Greve, M., Grundy, M., Guerrero, E., Hempel, J., Hengl, T., Heuvelink, G. B. M., Batjes, N., Carvalho, E., Hartemink, A., Hewitt, A., Hong, S.-Y., Krasilnikov, P., Lagacherie, P., Lelyk, G., Libohova, Z., Lilly, A., McBratney, A., McKenzie, N., Vasquez, G. M., Mulder, V. L., Minasny, B., Montanarella, L., Odeh, I., Padarian, J., Poggio, L., Roudier, P., Saby, N., Savin, I., Searle, R., Solbovoy, V., Thompson, J., Smith, S., Sulaeman, Y., Vintila, R., Rossel, R. V., Wilson, P., Zhang, G.-L., Swerts, M., Oorts, K., Karklins, A., Feng, L., Ibelles Navarro, A. R., Levin, A., Laktionova, T., Dell'Acqua, M., Suvannang, N., Ruam, W., Prasad, J., Patil, N., Husnjak, S., Pásztor, L., Okx, J., Hallett, S., Keay, C., Farewell, T., Lilja, H., Juilleret, J., Marx, S., Takata, Y., Kazuyuki, Y., Mansuy, N., Panagos, P., Van Liedekerke, M., Skalsky, R., Sobocka, J., Kobza, J., Eftekhari, K., Alavipanah, S. K., Moussadek, R., Badraoui, M., Da Silva, M., Paterson, G., Gonçalves, M. d. C., Theocharopoulos, S., Yemefack, M., Tedou, S., Vrscaj, B., Grob, U., Kozák, J., Boruvka, L., Dobos, E., Taboada, M., Moretti, L., and Rodriguez, D.: Soil legacy data rescue via GlobalSoilMap and other international and national initiatives, GeoResJ, 14, 1–19,, 2017. a

Arrouays, D., Mulder, V. L., and Richer-de Forges, A. C.: Soil mapping, digital soil mapping and soil monitoring over large areas and the dimensions of soil security – A review, Soil Security, 5, 100018,, 2021. a

Bakker, J., Dessel, B. v., and Zadelhoff, F. V.: Natuurwaardenkaart 1988: natuurgebieden, bossen en natte gronden in Nederland, 266862, s-Gravenhage SDU, ISBN 90-12-06089-3, (last access: 16 January 2024), 1989. a

Baltensweiler, A., Walthert, L., Hanewinkel, M., Zimmermann, S., and Nussbaum, M.: Machine learning based soil maps for a wide range of soil properties for the forested area of Switzerland, Geoderma Regional, 27, e00437,, 2021. a

Baumann, P., Helfenstein, A., Gubler, A., Keller, A., Meuli, R. G., Wächter, D., Lee, J., Viscarra Rossel, R., and Six, J.: Developing the Swiss mid-infrared soil spectral library for local estimation and monitoring, SOIL, 7, 525–546,, 2021. a

Been, T. H., Kempenaar, C., van Evert, F. K., Hoving, I. E., Kessel, G. J. T., Dantuma, W., Booij, J. A., Molendijk, L. P. G., Sijbrandij, F. D., and van Boheemen, K.: Akkerweb and farmmaps: Development of Open Service Platforms for Precision Agriculture, in: Precision Agriculture: Modelling, edited by: Cammarano, D., van Evert, F. K., and Kempenaar, C., Progress in Precision Agriculture, 269–293, Springer International Publishing, Cham, ISBN 978-3-031-15258-0,, 2023. a

Behrens, T., Förster, H., Scholten, T., Steinrücken, U., Spies, E.-D., and Goldschmitt, M.: Digital soil mapping using artificial neural networks, J. Plant Nutr. Soil Sci., 168, 21–33,, 2005. a

Behrens, T., Schmidt, K., MacMillan, R. A., and Viscarra Rossel, R. A.: Multi-scale digital soil mapping with deep learning, Sci. Rep., 8, 15244,, 2018a. a

Behrens, T., Schmidt, K., Rossel, R. A. V., Gries, P., Scholten, T., and MacMillan, R. A.: Spatial modelling with Euclidean distance fields and machine learning, Eur. J. Soil Sci., 69, 757–770,, 2018b. a

BIJ12: Informatiemodel Natuur (IMNa), (last access: 17 January 2024), 2019. a

Bouma, J. and Hartemink, A. E.: Soil science and society in the Dutch context, Netherlands Journal of Agricultural Science, 50, 133–140, 2003. a

Breiman, L.: Random Forests, Mach. Learn., 45, 5–32,, 2001. a

Breiman, L.: Manual on Setting Up, Using, and Understanding Random Forests v3.1, Technical report (last access: 14 January 2024), 2002. a

Brouwer, F. and van der Werff, M. M.: Vergraven gronden: inventarisatie van “diepe” grondbewerkingen, ophogingen en afgravingen, Alterra-rapport 2336, Alterra, Wageningen, (last access: 14 January 2024), 2012. a, b

Brouwer, F. and Walvoort, D.: Basisregistratie Ondergrond (BRO) - actualisatie bodemkaart: Herkartering van de bodem in Eemland, WOt-technical report 155, WOT Natuur & Milieu, Wageningen,, 2019. a

Brouwer, F. and Walvoort, D.: Basisregistratie Ondergrond (BRO) Actualisatie bodemkaart: Herkartering van de veengebieden aan de flanken van de Utrechtse Heuvelrug, WOt-technical report 177, WOT Natuur & Milieu, Wageningen,, 2020. a

Brouwer, F., de Vries, F. d., and Walvoort, D. J. J.: Basisregistratie Ondergrond (BRO) actualisatie bodemkaart: Herkartering van de bodem in Flevoland, WOt technical report 143, WOT Natuur & Milieu, WUR, Wageningen, (last access: 13 January 2024), 2018. a, b, c, d, e

Brouwer, F., Maas, G., Teuling, K., Harkema, T., and Verzandvoort, S.: Bodemkaart en Geomorfologische Kaart van Nederland: actualisatie 2020-2021 en toepassing: deelgebieden Gelderse Vallei-Zuid en -West en Veluwe-Zuid, WOt-rapport 134, WOT Natuur & Milieu, Wageningen,, 2021. a

Brouwer, F., Assinck, F., Harkema, T., Teuling, K., and Walvoort, D.: Actualisatie van de bodemkaart in degemeente Vijfheerenlanden: herkartering van de verbreiding van veen, WOt-rapport 151, WOT Natuur & Milieu, Wageningen, (last access: 25 January 2024), 2023. a, b

Brus, D., Hengl, T., Heuvelink, G., Kempen, B., Mulder, T. V., Olmedo, G. F., Poggio, L., Ribeiro, E., and Omuto, C. T.: Soil Organic Carbon Mapping Cookbook, edited by: Yigini, Y., Baritz, R., and Vargas, R. R., FAO, Rome, 1st edn., ISBN 978-92-5-130440-2, 2017. a

Brus, D. J.: Sampling for digital soil mapping: A tutorial supported by R scripts, Geoderma, 338, 464–480,, 2019. a

Brus, D. J.: Spatial sampling with R, The R Series, CRC Press, (last access: 16 December 2023), 2022. a, b

Brus, D. J. and Heuvelink, G. B. M.: Towards a Soil Information System with quantified accuracy: Three approaches for stochastic simulation of soil maps, Statutory Research Tasks Unit for Nature and the Environment 58, Alterra, Wageningen, 2007. a, b

Brus, D. J., Vašát, R., Heuvelink, G. B. M., Knotters, M., Vries, F. d., and Walvoort, D. J. J.: Towards a Soil Information System with quantified accuracy. A prototype for mapping continuous soil properties, Statutory Research Tasks Unit for Nature and the Environment 197, Alterra, Wageningen, 2009. a, b, c, d, e, f, g, h

Brus, D. J., Kempen, B., and Heuvelink, G. B. M.: Sampling for validation of digital soil maps, Eur. J. Soil Sci., 62, 394–407,, 2011. a, b, c, d

Buchanan, S., Triantafilis, J., Odeh, I. O. A., and Subansinghe, R.: Digital soil mapping of compositional particle-size fractions using proximal and remotely sensed ancillary data, Geophysics, 77, WB201–WB211,, 2012. a

Buringh, P., Stuer, G. G. L., and Vink, P.: Some techniques and methods of soil survey in the Netherlands, Netherlands Journal of Agricultural Science, 10, 17, 1962. a

BZK: Ministerie van Binnenlandse Zaken en Koninkrijksrelaties (BZK): Uitvoeringsregeling Meststoffenwet, (last access: 17 January 2024), 2022. a

Carré, F., McBratney, A. B., and Minasny, B.: Estimation and potential improvement of the quality of legacy soil samples for digital soil mapping, Geoderma, 141, 1–14,, 2007. a

Chagas, C. D. S., de Carvalho Junior, W., Bhering, S. B., and Calderano Filho, B.: Spatial prediction of soil surface texture in a semiarid region using random forest and multiple linear regressions, CATENA, 139, 232–240,, 2016. a

Chen, S., Arrouays, D., Leatitia Mulder, V., Poggio, L., Minasny, B., Roudier, P., Libohova, Z., Lagacherie, P., Shi, Z., Hannam, J., Meersmans, J., Richer-de Forges, A. C., and Walter, C.: Digital mapping of GlobalSoilMap soil properties at a broad scale: A review, Geoderma, 409, 115567,, 2022. a, b, c, d, e, f, g, h, i

Chen, S., Saby, N. P. A., Martin, M. P., Barthès, B. G., Gomez, C., Shi, Z., and Arrouays, D.: Integrating additional spectroscopically inferred soil data improves the accuracy of digital soil mapping, Geoderma, 433, 116467,, 2023. a

Clement, J.: GIS Vierde Bosstatistiek: Gebruikersdocumentatie, Documentatie van bestanden, Tech. rep., Research Instituut voor de Groene Ruimte, Alterra, Wageningen, 2001. a

Cochran, W. G.: Sampling techniques, John Wiley & Sons, New York, 3rd edn., 1977. a

Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L., Wehberg, J., Wichmann, V., and Böhner, J.: System for Automated Geoscientific Analyses (SAGA) v. 2.1.4, Geosci. Model Dev., 8, 1991–2007,, 2015. a

de Bakker, H. and Schelling, J.: Systeem van bodemclassificatie voor Nederland: De hogere niveaus, Centrum voor Landbouwpublikaties en Landbouwdocumentatie, Wageningen, the Netherlands, 1st edn., first version, (last access: 16 December 2023), 1966. a, b, c, d

de Bakker, H. and Schelling, J.: Systeem van bodemclassificatie voor Nederland: de hogere niveaus: With Engl. summary: A system of soil classification for the Netherlands, Centrum voor Landbouwpublikaties en Landbouwdocumentatie, Wageningen, the Netherlands, second, revised edn., ISBN 978-90-220-0997-0, 1989. a, b, c, d, e

de Bruin, S., Bregt, A., and Ven, M. V. D.: Assessing fitness for use: the expected value of spatial data sets, Int. J. Geogr. Inf., 15, 457–471,, 2001. a

de Bruin, S., Brus, D. J., Heuvelink, G. B. M., van Ebbenhorst Tengbergen, T., and Wadoux, A. M. J.-C.: Dealing with clustered samples for assessing map accuracy by cross-validation, Ecol. Inform., 69, 101665,, 2022. a, b

de Gruijter, J. J., Walvoort, D. J. J., and van Gams, P. F. M.: Continuous soil maps – a fuzzy set approach to bridge the gap between aggregation levels of process and distribution models, Geoderma, 77, 169–195,, 1997. a

de Gruijter, J. J., van der Horst, J. B. F., Heuvelink, G. B. M., Knotters, M., and Hoogland, T.: Grondwater opnieuw op de kaart; methodiek voor de actualisering van grondwaterstandsinformatie en perceelsclassificatie naar uitspoelingsgevoeligheid voor nitraat, Alterra-rapport 915, Alterra, Wageningen, (last access: 16 December 2023), 2004. a, b

de Gruijter, J. J., Brus, D., Bierkens, M., and Knotters, M.: Sampling for Natural Resource Monitoring, Springer, The Netherlands, ISBN 978-3-540-33161-2, 2006. a

de Vries, F. and Al, E. J.: De groeiplaatsgeschiktheid voor bosdoeltypen in beeld met ALBOS, Tech. Rep. 234, DLO-Staring Centrum, (last access: 16 December 2023), 1992. a

de Vries, F., de Groot, W., Hoogland, T., and Denneboom, J.: De Bodemkaart van Nederland digitaal;. Toelichting bij inhoud, actualiteit en methodiek en korte beschrijving van additionel informatie, Alterra-rapport 811, Alterra, Research Instituut voor de Groene Ruimte, Wageningen, the Netherlands, 2003. a, b, c, d, e, f, g, h, i, j

de Vries, F., Brus, D. J., Kempen, B., Brouwer, F., and Heidema, A. H.: Actualisatie bodemkaart veengebieden: deelgebied en 2 in Noord Nederland, Alterra-rapport 2556, Alterra, Wageningen, (last access: 16 December 2023), 2014. a

de Vries, F., Walvoort, D., and Brouwer, F.: Basisregistratie Ondergrond (BRO) Actualisatie bodemkaart; Herkartering van de eenheden met slappe kleilagen, Wageningen Environmental Research Rapport 2834, Wageningen Environmental Research, Wageningen,, 2017. a

de Vries, F., Brouwer, F., and Walvoort, D.: Basisregistratie Ondergrond (BRO) actualisatie bodemkaart: Herkatering westelijk veengebied Waterschap Drents Overijsselse Delta, Wageningen Environmental Research Rapport 2887, Wageningen Environmental Research, Wageningen,, 2018. a

Delmas, M., Saby, N., Arrouays, D., Dupas, R., Lemercier, B., Pellerin, S., and Gascuel-Odoux, C.: Explaining and mapping total phosphorus content in French topsoils, Soil Use Manage., 31, 259–269,, 2015. a

Delta Programme: National Delta Programme 2024: Now for the Future, Tech. rep., Ministry of Infrastructure and Water Management, the Ministry of Agriculture, Nature and Food Quality, and the Ministry of the Interior and Kingdom Relations, Den Haag, (last access: 16 December 2023), 2023. a

Dokuchaev, V.: Report to the Transcaucasian Statistical Committee on Land Evaluation in General and Especially for the Transcaucasia. Horizontal and Vertical Soil Zones, Off. Press Civ, Affairs Commander-in-Chief Cacasus, Tiflis, Russia, 1899 (in Russian). a, b, c

Edelmann, C.: Soils of the Netherlands, vol. 53, North-Holland Publishing Company, Amsterdam,, 1950. a, b

Erisman, J. W.: Setting ambitious goals for agriculture to meet environmental targets, One Earth, 4, 15–18,, 2021. a

Eurofins Agro: Bemesting Wetgeving, (last access: 16 December 2023), 2024a. a

Eurofins Agro: BemestingsWijzer, (last access: 16 December 2023), 2024b. a

European Commission: A Soil Deal for Europe: 100 living labs and lighthouses to lead the transition towards healthy soils by 2030, Implementation Plan, European Commission, (last access: 16 December 2023), 2021. a, b, c

European Commission: Proposal for a Directive on Soil Monitoring and Resilience, Tech. rep., European Commission, (last access: 16 December 2023), 2023. a, b

Eurostat: Key figures on the European food chain, (last access: 16 December 2023), 2022. a

EZK: Fysisch Geografische Regio's 2013; Ministerie van Economische Zaken en Klimaat (EZK; Ministry of Economic Affairs and Climate), (last access: 16 December 2023), 2013. a

EZK: Basisregistratie Gewaspercelen (BRP): 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019. Ministerie van Economische Zaken en Klimaat (EZK; Ministry of Economic Affairs and Climate), Agrarische Areaal Nederland, (last access: 16 December 2023), 2019. a

Felix, R.: Bodemkartering voor 1943 – het geologisch perspectief, in: Van bodemkaart tot informatiesysteem, edited by: Buurman, P. and Sevink, J., 1–17, Wageningen Pers, Wageningen, 1995. a

Finke, P. A.: On digital soil assessment with models and the Pedometrics agenda, Geoderma, 171–172, 3–15,, 2012. a

Finke, P. A., de Gruijter, J. J., and Visschers, R.: Status 2001 Landelijke Steekproef Kaarteenheden en toepassingen; Gestructureerde bemonstering en karakterisering Nederlandse bodems, Alterra-rapport 389, Alterra, Research Instituut voor de Groene Ruimte, Wageningen, 2001. a

GDAL/OGR contributors: GDAL/OGR geospatial data abstraction software library, (last access: 16 December 2023), 2023. a

Gregoire, T. G. and Valentine, H. T.: Sampling Strategies for Natural Resources and the Environment, CRC Press, Boca Raton, USA, ISBN 978-0-203-49888-0, 2007. a

Grimm, R. and Behrens, T.: Uncertainty analysis of sample locations within digital soil mapping approaches, Geoderma, 155, 154–163,, 2010. a

Guyon, I., Weston, J., Barnhill, S., and Vapnik, V.: Gene Selection for Cancer Classification using Support Vector Machines, Mach. Learn., 46, 389–422,, 2002. a

Hack-ten Broeke, M., van Beek, C. L., Hoogland, T., Knotters, M., Mol-Dijkstra, J. P., Schils, R., Smit, A., and de Vries, F.: Kaderrichtlijn Bodem; Basismateriaal voor eventuele prioritaire gebieden, Alterra-rapport 2007, Alterra, Wageningen, the Netherlands, 2009. a

Hack-ten Broeke, M. J. D., Mulder, H. M., Bartholomeus, R. P., van Dam, J. C., Holshof, G., Hoving, I. E., Walvoort, D. J. J., Heinen, M., Kroes, J. G., van Bakel, P. J. T., Supit, I., de Wit, A. J. W., and Ruijtenberg, R.: Quantitative land evaluation implemented in Dutch water management, Geoderma, 338, 536–545,, 2019. a

Hartemink, A. E. and Sonneveld, M. P. W.: Soil maps of The Netherlands, Geoderma, 204–205, 1–9,, 2013. a

Hastie, T., Tibshirani, R., and Friedman, J. H.: The elements of statistical learning: data mining, inference, and prediction, Springer series in statistics, Springer, New York, NY, 2nd edn., ISBN 978-0-387-84857-0, 978-0-387-84858-7, 2009. a

Hazeu, G., Schuiling, R., Thomas, D., Vittek, M., Storm, M., and Bulens, J. D.: Landelijk Grondgebruiksbestand Nederland 2021 (LGN2021): achtergronden, methodiek en validatie, Rapport 3235, Wageningen Environmental Research, Wageningen, (last access: 16 December 2023), 2023. a

Hazeu, G. W., Vittek, M., Schuiling, R., Bulens, J. D., Storm, M. H., Roerink, G. J., and Meijninger, W. M. L.: LGN2018: een nieuwe weergave van het grondgebruik in Nederland, Tech. Rep. 3010, Wageningen Environmental Research, Wageningen, (last access: 16 December 2023), 2020. a, b, c

Heinen, M., Mulder, H. M., Bakker, G., Wösten, J. H. M., Brouwer, F., Teuling, K., and Walvoort, D. J. J.: The Dutch soil physical units map: BOFEK, Geoderma, 427, 116123,, 2022. a, b, c, d, e, f, g, h

Helfenstein, A.: 4 Dimensional Information About the Skin of the Earth, Wageningen University & Research [video],, last access: 21 June 2024a. a

Helfenstein, A.: BIS-4D. In Earth System Science Data, Zenodo [code],, 2024b. a

Helfenstein, A., Baumann, P., Viscarra Rossel, R., Gubler, A., Oechslin, S., and Six, J.: Quantifying soil carbon in temperate peatlands using a mid-IR soil spectral library, SOIL, 7, 193–215,, 2021. a

Helfenstein, A., Mulder, V. L., Heuvelink, G. B., and Okx, J. P.: Tier 4 maps of soil pH at 25 m resolution for the Netherlands, Geoderma, 410, 115659,, 2022. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

Helfenstein, A., Mulder, V. L., Hack-ten Broeke, M. J., van Doorn, M., Teuling, K., Walvoort, D. J., and Heuvelink, G. B.: BIS-4D: Maps of soil properties and their uncertainties at 25 m resolution in the Netherlands, 4TU.ResearchData [data set],, 2024a. a, b, c

Helfenstein, A., Mulder, V. L., Hack-ten Broeke, M. J., van Doorn, M., Teuling, K., Walvoort, D. J., and Heuvelink, G. B.: Spatially explicit environmental variables at 25 m resolution for spatial modelling in the Netherlands, 4TU.ResearchData [data set],, 2024b. a, b

Helfenstein, A., Mulder, V. L., Heuvelink, G. B. M., and Hack-ten Broeke, M. J. D.: Three-dimensional space and time mapping reveals soil organic matter decreases across anthropogenic landscapes in the Netherlands, Commun. Earth Environ., 5, 1–16,, 2024c. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s

Helfenstein, A., Teuling, K., Walvoort, D. J., Hack-ten Broeke, M. J., Mulder, V. L., van Doorn, M., and Heuvelink, G. B.: Georeferenced point data of soil properties in the Netherlands, 4TU.ResearchData [data set],, 2024d. a

Hengl, T. and MacMillan, R. A.: Predictive Soil Mapping with R, OpenGeoHub foundation, Wageningen, the Netherlands, ISBN 978-0-359-30635-0, (last access: 17 January 2024), 2019. a

Hengl, T., Heuvelink, G. B. M., Kempen, B., Leenaars, J. G. B., Walsh, M. G., Shepherd, K. D., Sila, A., MacMillan, R. A., Mendes de Jesus, J., Tamene, L., and Tondoh, J. E.: Mapping Soil Properties of Africa at 250 m Resolution: Random Forests Significantly Improve Current Predictions, PLOS ONE, 10, e0125814,, 2015. a

Hengl, T., Nussbaum, M., Wright, M. N., Heuvelink, G. B. M., and Gräler, B.: Random forest as a generic framework for predictive modeling of spatial and spatio-temporal variables, PeerJ, 6, e5518,, 2018. a

Hessel, R., Stolte, J., and Riksen, M. J. P. M.: Huidige maatregelen tegen water- en winderosie in Nederland, Bodem, 21, 11–12, (last access: 17 January 2024), 2011. a

Heuvelink, G., Brus, D., De Vries, F., Kempen, B., Knotters, M., Vasat, R., and Walvoort, D.: Implications of digital soil mapping for soil information systems, in: 4th Global Workshop on Digital Soil Mapping, p. 6, Rome, Italy, (last access: 17 January 2024), 2010. a

Heuvelink, G. B. M.: Uncertainty quantification of GlobalSoilMap products, in: GlobalSoilMap: Basis of the global spatial soil information system, 335–340, CRC Press, ISBN 978-1-315-77558-6, 2014. a

Heuvelink, G. B. M.: Uncertainty and Uncertainty Propagation in Soil Mapping and Modelling, in: Pedometrics, edited by: McBratney, A. B., Minasny, B., and Stockmann, U., Progress in Soil Science, 439–461, Springer International Publishing, Cham, ISBN 978-3-319-63439-5,, 2018. a, b, c

Heuvelink, G. B. M., Angelini, M. E., Poggio, L., Bai, Z., Batjes, N. H., Bosch, R. v. d., Bossio, D., Estella, S., Lehmann, J., Olmedo, G. F., and Sanderman, J.: Machine learning in space and time for modelling soil organic carbon change, Eur. J. Soil Sci., 72, 1607–1623,, 2020. a

Hoogland, T., Knotters, M., Pleijter, M., and Walvoort, D. J. J.: Actualisatie van de grondwatertrappenkaart van holoceen Nederland: resultaten van het veldonderzoek, Alterra-rapport 2612, Alterra Wageningen UR, Wageningen, (last access: 17 January 2024), 2014. a, b

Jenny, H.: Factors of Soil Formation: A System of Quantitative Pedology, McGraw- Hill, New York, ISBN 978-0-486-68128-3, (last access: 17 January 2024), 1941. a, b, c

Jukema, G., Ramaekers, P., and Berkhout, P.: De Nederlandse agrarische sector in internationaal verband: Editie 2023, Wageningen Economic Research, ISBN 978-94-6447-546-3,, 2023. a

Keller, A., Franzen, J., Knüsel, P., Papritz, A., and Zürrer, M.: Bodeninformations-Plattform Schweiz (BIP-CH), Thematische Synthese TS4 des Nationalen Forschungsprogramms “Nachhaltige Nutzung der Ressource Boden” (NFP 68), Schweizerischer Nationalfonds zur Förderung der wissenschaftlichen Forschung (SNF), Bern, 2018. a

Kempen, B., Brus, D. J., Heuvelink, G. B. M., and Stoorvogel, J. J.: Updating the 1:50 000 Dutch soil map using legacy soil data: A multinomial logistic regression approach, Geoderma, 151, 311–326,, 2009. a, b

Kempen, B., Brus, D. J., and Stoorvogel, J. J.: Three-dimensional mapping of soil organic matter content using soil type–specific depth functions, Geoderma, 162, 107–123,, 2011. a

Kempen, B., Brus, D. J., and Heuvelink, G. B. M.: Soil type mapping using the generalised linear geostatistical model: A case study in a Dutch cultivated peatland, Geoderma, 189–190, 540–553,, 2012a. a

Kempen, B., Brus, D. J., Stoorvogel, J. J., Heuvelink, G. B. M., and de Vries, F.: Efficiency Comparison of Conventional and Digital Soil Mapping for Updating Soil Maps, Soil Sci. Soc. Am. J., 76, 2097–2115,, 2012b. a

Kempen, B., Heuvelink, G. B. M., Brus, D., and Walvoort, D.: Towards products for The Netherlands, GlobalSoilMap: Basis of the Global Spatial Soil Information System – Proceedings of the 1st GlobalSoilMap Conference, 85–90,, 2014. a

Kempen, B., Brus, D. J., and de Vries, F.: Operationalizing digital soil mapping for nationwide updating of the 1:50,000 soil map of the Netherlands, Geoderma, 241–242, 313–329,, 2015. a

Keskin, H., Grunwald, S., and Harris, W. G.: Digital mapping of soil carbon fractions with machine learning, Geoderma, 339, 40–58,, 2019. a

Khaledian, Y. and Miller, B. A.: Selecting appropriate machine learning methods for digital soil mapping, Appl. Math. Model., 81, 401–418,, 2020. a

KNMI: Koninklijk Nederlands Meteorologisch Instituut (KNMI) Dataplatform, Koninklijk Nederlands Meteorologisch Instituut, (last access: 12 March 2021), 2020. a, b

Knotters, M. and Vroon, H. R. J.: The economic value of detailed soil survey in a drinking water collection area in the Netherlands, Geoderma Regional, 5, 44–53,, 2015. a, b

Knotters, M. and Walvoort, D.: Hoge resolutie, nauwkeurige kaarten?, (last access: 16 October 2020), 2020. a

Knotters, M., Broeke, M. J. D. H.-t., Hinssen, P. J. W., Kolk, J. W. H. v. d., and Okx, J. P.: Betekenis van BRO/BIS Nederland voor WOT Natuur & Milieu: Een risicoanalyse, WOT Natuur & Milieu, (last access: 17 January 2024), 2015a. a

Knotters, M., Okx, J., Hack-ten Broeke, M., and de Vries, F.: Bodem in beweging: BIS NEderland informeert, Bodem, 25, 11–13, 2015b. a

Knotters, M., Walvoort, D. J. J., Brouwer, F., Stuyt, L. C. P. M., and Okx, J. P.: Landsdekkende, actuele informatie over grondwatertrappen digitaal beschikbaar, H2O online, (last access: 17 January 2024), 2018. a, b, c, d

Knotters, M., Teuling, K., Reijneveld, A., Lesschen, J. P., and Kuikman, P.: Changes in organic matter contents and carbon stocks in Dutch soils, 1998–2018, Geoderma, 414, 115751,, 2022. a, b, c, d, e, f, g

Kombrink, H., Doornenbal, J. C., Duin, E. J. T., Dulk, M. d., Veen, J. H. t., and Witmans, N.: New insights into the geological structure of the Netherlands; results of a detailed mapping project, Netherlands Journal of Geosciences, 91, 419–446,, 2012. a

Koomen, A. and Maas, G.: Geomorfologische Kaart Nederland (GKN); Achtergronddocument bij het landsdekkende digitale bestand, Altera-rapport 1039, Alterra, Wageningen, 2004. a

Kortleve, A. J., Mogollón, J. M., Heimovaara, T. J., and Gebert, J.: Topsoil Carbon Stocks in Urban Greenspaces of The Hague, the Netherlands, Urban Ecosystems, 26, 725–742,, 2023. a

Kroes, J., Van Dam, J., Bartholomeus, R., Groenendijk, P., Heinen, M., Hendriks, R., Mulder, H., Supit, I., and Van Walsum, P.: SWAP version 4; Theory description and user manual, Report 2780, Wageningen Environmental Research, Wageningen,, 2017. a

Kuhn, M.: Building Predictive Models in R Using the Caret Package, J. Stat. Softw., 28, 1–26,, 2008. a

Kuhn, M.: The caret Package, GitHub [code], (last access: 17 January 2024), 2019. a, b

Kuhn, M.: Classification and Regression Training: package “caret”, GitHub [code],, 2022. a

Kull, A., Kikas, T., Penu, P., and Kull, A.: Modeling Topsoil Phosphorus – From Observation-Based Statistical Approach to Land-Use and Soil-Based High-Resolution Mapping, Agronomy, 13, 1183,, 2023. a

Lamigueiro, O. P. and Hijmans, R. J.: Visualization Methods for Raster Data: package “rasterVis”, GitHub [code], (last access: 17 January 2024), 2023. a

Lark, R. M. and Bishop, T. F. A.: Cokriging particle size fractions of the soil, Eur. J. Soil Sci., 58, 763–774,, 2007. a

Lehmann, J., Bossio, D. A., Kögel-Knabner, I., and Rillig, M. C.: The concept and future prospects of soil health, Nat. Rev. Earth Environ., 1, 544–553,, 2020. a

Lemercier, B., Lagacherie, P., Amelin, J., Sauter, J., Pichelin, P., Richer-de Forges, A. C., and Arrouays, D.: Multiscale evaluations of global, national and regional digital soil mapping products in France, Geoderma, 425, 116052,, 2022. a

Li, J., Heap, A. D., Potter, A., and Daniell, J. J.: Application of machine learning methods to spatial interpolation of environmental variables, Environ. Model. Softw., 26, 1647–1659,, 2011. a

Liu, L., Ji, M., and Buchroithner, M.: Transfer Learning for Soil Spectroscopy Based on Convolutional Neural Networks and Its Application in Soil Clay Content Mapping Using Hyperspectral Imagery, Sensors, 18, 3169,, 2018. a

Loiseau, T., Chen, S., Mulder, V. L., Román Dobarco, M., Richer-de Forges, A. C., Lehmann, S., Bourennane, H., Saby, N. P. A., Martin, M. P., Vaudour, E., Gomez, C., Lagacherie, P., and Arrouays, D.: Satellite data integration for soil clay content modelling at a national scale, Int. J. Appl. Earth Ob., 82, 101905,, 2019. a

Lookman, R., Vandeweert, N., Merckx, R., and Vlassak, K.: Geostatistical assessment of the regional distribution of phosphate sorption capacity parameters (FeOX and AlOX) in northern Belgium, Geoderma, 66, 285–296,, 1995. a, b

Ma, Y., Minasny, B., McBratney, A., Poggio, L., and Fajardo, M.: Predicting soil properties in 3D: Should depth be a covariate?, Geoderma, 383, 114794,, 2021. a, b

Maas, G., van der Meij, M., Delft, S., and Heidema, A.: Toelichting bij de legenda Geomorfologische kaart van Nederland 1:50000 (2019), Wageningen Environmental Research, Wageningen, (last access: 17 January 2024), 2019. a, b

Malone, B., Searle, R., Malone, B., and Searle, R.: Updating the Australian digital soil texture mapping (Part 1*): re-calibration of field soil texture class centroids and description of a field soil texture conversion algorithm, Soil Res., 59, 419–434,, 2021. a

Malone, B. P., McBratney, A. B., and Minasny, B.: Spatial Scaling for Digital Soil Mapping, Soil Sci. Soc. Am. J., 77, 890–902,, 2013. a

Malone, B. P., Minasny, B., and McBratney, A. B.: Using R for Digital Soil Mapping, Progress in Soil Science, Springer International Publishing, Cham, ISBN 978-3-319-44325-6 978-3-319-44327-0,, 2017. a

Maring, L., Vries, F. d., Brouwer, F., Groot, H., Kiden, P., Leeters, E. E. J. M., and Mol, G.: IMBOD deelactiviteit 5: inhoudelijke afstemming, Alterra-rapport 1817, Alterra, Wageningen, (last access: 17 January 2024), 2009. a, b, c, d

Matos-Moreira, M., Lemercier, B., Dupas, R., Michot, D., Viaud, V., Akkal-Corfini, N., Louis, B., and Gascuel-Odoux, C.: High-resolution mapping of soil phosphorus concentration in agricultural landscapes with readily available or detailed survey data, Eur. J. Soil Sci., 68, 281–294,, 2017. a

McBratney, A., Mendonça Santos, M., and Minasny, B.: On digital soil mapping, Geoderma, 117, 3–52,, 2003. a, b

McKeague, J. A.: An evaluation of 0.1 m pyrophosphate and pyrophosphate-dithionite in comparison with oxalate as extractants of the accumulation products in podzols and some other soils, Ca. J. Soil Sci., 47, 95–99,, 1967. a

McKeague, J. A., Brydon, J. E., and Miles, N. M.: Differentiation of Forms of Extractable Iron and Aluminum in Soils, Soil Sci. Soc. Am. J., 35, 33–38,, 1971. a

Meinshausen, N.: Quantile Regression Forests, J. Mach. Learn. Res., 7, 983–999, 2006. a

Meyer, H.: “caret” Applications for Spatial-Temporal Models: package “CAST”, GitHub [code], (last access: 13 January 2024), 2023. a

Mulder, M., Walvoort, D., Brouwer, F., Tol-Leenders, D. V., and Verzandvoort, S.: Bodemgeschiktheidskaarten voor landbouw in de provincie Noord-Brabant: Een toepassing van Waterwijzer Landbouw, Rapport 3206, Wageningen Environmental Research, Wageningen, the Netherlands, (last access: 18 January 2024), 2022. a

Mulder, V. L., Lacoste, M., Richer-de Forges, A. C., and Arrouays, D.: GlobalSoilMap France: High-resolution spatial modelling the soils of France up to two meter depth, Sci. Total Environ., 573, 1352–1369,, 2016. a

Møller, A. B., Beucher, A. M., Pouladi, N., and Greve, M. H.: Oblique geographic coordinates as covariates for digital soil mapping, SOIL, 6, 269–289,, 2020. a, b

NEN 5753: Soil – Determination of clay content and particle size distribution in soil and sediment by sieve and pipet, Standards ICS codes: 13.080.20, Anorganische parameters (Milieukwaliteit), (last access: 18 January 2024), 2018 edn., 2020 (in Dutch). a, b, c

Neyroud, J.-A. and Lischer, P.: Do different methods used to estimate soil phosphorus availability across Europe give comparable results?, J. Plant Nutr. Soil Sci., 166, 422–431,, 2003. a

NHI: Nederlands Hydrologisch Instrumentarium (NHI), (last access: 18 January 2024), 2023. a

Nussbaum, M., Walthert, L., Fraefel, M., Greiner, L., and Papritz, A.: Mapping of soil properties at high resolution in Switzerland using boosted geoadditive models, SOIL, 3, 191–210,, 2017. a

Nussbaum, M., Zimmermann, S., Walthert, L., and Baltensweiler, A.: Benefits of hierarchical predictions for digital soil mapping – An approach to map bimodal soil pH, Geoderma, 437, 116579,, 2023. a

Odeh, I. O. A., Todd, A. J., and Triantafilis, J.: Spatial prediction of soil particle-size fractions as compositional data, Soil Sci,, 168, 501,, 2003. a

Padarian, J., Minasny, B., and McBratney, A. B.: Transfer learning to localise a continental soil vis-NIR calibration model, Geoderma, 340, 279–288,, 2019a. a

Padarian, J., Minasny, B., and McBratney, A. B.: Using deep learning for digital soil mapping, SOIL, 5, 79–89,, 2019b. a, b

Pahlavan-Rad, M. R. and Akbarimoghaddam, A.: Spatial variability of soil texture fractions and pH in a flood plain (case study from eastern Iran), CATENA, 160, 275–281,, 2018. a

Pan, S. J. and Yang, Q.: A Survey on Transfer Learning, IEEE T. Knowl. Data Eng., 22, 1345–1359,, 2010. a

Panagos, P., Montanarella, L., Barbero, M., Schneegans, A., Aguglia, L., and Jones, A.: Soil priorities in the European Union, Geoderma Regional, 29, e00510,, 2022. a, b

Papadopoulos, G., Edwards, P., and Murray, A.: Confidence estimation methods for neural networks: a practical comparison, IEEE T. Neur. Netw., 12, 1278–1287,, 2001. a

Pawlowsky-Glahn, V. and Buccianti, A.: Compositional Data Analysis: Theory and Applications, John Wiley & Sons, Ltd, West Sussex, ISBN 978-1-119-97646-2,, 2011. a

Pawlowsky-Glahn, V., Egozcue, J. J., and Tolosana-Delgado, R.: Modeling and Analysis of Compositional Data, John Wiley & Sons, West Sussex, ISBN 978-1-118-44306-4, 2015. a

Piikki, K., Wetterlind, J., Söderström, M., and Stenberg, B.: Perspectives on validation in digital soil mapping of continuous attributes – A review, Soil Use Manage., 37, 7–21,, 2021. a, b

Poggio, L. and Gimona, A.: 3D mapping of soil texture in Scotland, Geoderma Regional, 9, 5–16,, 2017. a

Poggio, L., de Sousa, L. M., Batjes, N. H., Heuvelink, G. B. M., Kempen, B., Ribeiro, E., and Rossiter, D.: SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty, SOIL, 7, 217–240,, 2021. a, b, c

QGIS Development Team: QGIS: A Free and Open Source Geographic Information System, (last access: 18 January 2024), 2023. a

R Core Team: R: A language and environment for statistical computing, Tech. rep., R Foundation for Statistical Computing, Vienna, Austria, (last access: 18 January 2024), 2023. a

RIVM: Grootschalige Concentratie- en Depositiekaarten Nederland (GCN, GDN), Rijksinstituut voor Volksgezondheid en Milieu (RIVM), (last access: 18 January 2024), 2020. a

Roerink, G. and Mücher, S.: Groenmonitor, (last access: 18 January 2024), 2023. a, b

Román Dobarco, M., Arrouays, D., Lagacherie, P., Ciampalini, R., and Saby, N. P. A.: Prediction of topsoil texture for Region Centre (France) applying model ensemble methods, Geoderma, 298, 67–77,, 2017. a

Römkens, P. F. A. M. and Oenema, O.: Quick scan soils in the Netherlands; overview of the soil status with reference to the forthcoming EU soil strategy, Alterra-rapport 948, Alterra, Wageningen, (last access: 18 January 2024), 2004. a, b, c, d

Ros, G. H.: Open Bodemindex (OBI): Eenvoudig en betaalbaar inzicht in bodemkwaliteit en bodemverbetering., (last access: 18 January 2024), 2023. a, b

Ros, G. H., Verweij, S. E., Janssen, S. J. C., De Haan, J., and Fujita, Y.: An Open Soil Health Assessment Framework Facilitating Sustainable Soil Management, Environ. Sci. Technol., 56, 17375–17384,, 2022. a, b

Ros, G. H., Haan, J. J. d., Fuchs, L. M., and Molendijk, L.: Bodembeoordeling van landbouwgronden voor diverse ecosysteemdiensten: ontwikkeling van de BLN, versie 2.0, Rapport / Wageningen University & Research, Business unit Open Teelten WPR-OT-1030, Wageningen Plant Research, Wageningen,, 2023. a, b

Sandri, M. and Zuccolotto, P.: A Bias Correction Algorithm for the Gini Variable Importance Measure in Classification Trees, J. Comput. Graph. Stat., 17, 611–628,, 2008. a

Sandri, M. and Zuccolotto, P.: Analysis and correction of bias in Total Decrease in Node Impurity measures for tree-based algorithms, Stat. Comput., 20, 393–407,, 2010. a

Schmidinger, J. and Heuvelink, G. B. M.: Validation of uncertainty predictions in digital soil mapping, Geoderma, 437, 116585,, 2023. a, b, c

Schoumans, O. F.: Description of the Phosphorus Sorption and Desorption Processes in Lowland Peaty Clay Soils, Soil Sci., 178, 291,, 2013. a

Schoumans, O. F. and Chardon, W. J.: Phosphate saturation degree and accumulation of phosphate in various soil types in The Netherlands, Geoderma, 237–238, 325–335,, 2015. a, b, c

Schulte, R. P. O., Bampa, F., Bardy, M., Coyle, C., Creamer, R. E., Fealy, R., Gardi, C., Ghaley, B. B., Jordan, P., Laudon, H., O'Donoghue, C., Ó'hUallacháin, D., O'Sullivan, L., Rutgers, M., Six, J., Toth, G. L., and Vrebos, D.: Making the Most of Our Land: Managing Soil Functions from Local to Continental Scale, Front. Environ. Sci., 3, 81,, 2015. a

Scull, P., Franklin, J., Chadwick, O. A., and McArthur, D.: Predictive soil mapping: a review, Prog. Phys. Geogr. Earth and Environment, 27, 171–197,, 2003. a, b

Seidel, M., Hutengs, C., Ludwig, B., Thiele-Bruhn, S., and Vohland, M.: Strategies for the efficient estimation of soil organic carbon at the field scale with vis-NIR spectroscopy: Spectral libraries and spiking vs. local calibrations, Geoderma, 354, 113856,, 2019. a

Sekulić, A., Kilibarda, M., Heuvelink, G. B. M., Nikolić, M., and Bajat, B.: Random Forest Spatial Interpolation, Remote Sens.-Basel, 12, 1687,, 2020. a

Slier, T., Mi-Gegotek, Y., and Lesschen, J.: Potentie voor koolstofvastlegging in minerale landbouwbodems in de Nederlandse provincies, Slim Landgebruik, Wageningen Environmental Research, Wageningen, (last access: 18 January 2024), 2023. a, b

Stettler, M., Keller, T., Weisskopf, P., Lamandé, M., Lassen, P., and Schjønning, P.: Terranimo – a web-based tool for evaluating soil compaction, Landtechnik, 69, 132–138, 2014. a

Stokstad, E.: Nitrogen crisis threatens Dutch environment – and economy, Science, 366, 1180–1181,, 2019. a

Szatmári, G., Pásztor, L., and Heuvelink, G. B. M.: Estimating soil organic carbon stock change at multiple scales using machine learning and multivariate geostatistics, Geoderma, 403, 115356,, 2021. a, b

Taghizadeh-mehrjardi, R., Toomanian, N., Khavaninzadeh, A. R., Jafari, A., and Triantafilis, J.: Predicting and mapping of soil particle-size fractions with adaptive neuro-fuzzy inference and ant colony optimization in central Iran, Eur. J. Soil Sci., 67, 707–725,, 2016. a

ten Cate, J., van Holst, A., Kleijer, H., and Stolp, J.: Handleiding bodemgeografisch onderzoek: richtlijnen en voorschriften. Deel A: Bodem, Technisch Document 19A, DLO-Staring Centrum, Wageningen, (last access: 19 January 2024), 1995. a, b, c, d

Teuling, K., Knotters, M., van Tol-Leenders, T. P., Lesschen, J. P., and Reijneveld, J. A.: Nieuwe steekproefomvang voor landelijke monitoring koolstof en bodemkwaliteit, Vervolg op rapportages CC-NL en De staat van de Nederlandse landbouwbodems in 2018, Tech. rep., Wageningen Environmental Research, Wageningen, 2021. a

Tobler, W. R.: A Computer Movie Simulating Urban Growth in the Detroit Region, Economic Geography, 46, 234–240,, 1970. a, b

United Nations: Transforming our World: The 2030 Agenda for Sustainable Development, Tech. rep., United Nations, New York, NY, 2015. a

van Asselen, S., Erkens, G., Stouthamer, E., Woolderink, H. A. G., Geeraert, R. E. E., and Hefting, M. M.: The relative contribution of peat compaction and oxidation to subsidence in built-up areas in the Rhine-Meuse delta, The Netherlands, Sci. Total Environ., 636, 177–191,, 2018. a

van Dam, J. C., Huygen, J., Wesseling, J. G., Feddes, R. A., and Kabat, P.: Theory of SWAP version 2.0; Simulation of waterflow, solute transport and plant growth in the Soil-Water-Atmosphere-Plant environment, Technical document 45, Wageningen Agricultural University and DLO Winand Staring Centre, Wageningen, the Netherlands, (last access: 19 January 2024), 1997. a

van Delft, B. and Maas, G.: Landschappelijke Bodemkartering (LBK): Achtergronden, toepassingen en technische documentatie, WOt-technical report 248, Wettelijke Onderzoekstaken Natuur & Milieu, Wageningen,, 2023. a, b

van Delft, S. P. J. and Maas, G. J.: De Landschappelijke Bodemkaart van Nederland; versie 2022, (last access: 19 January 2024), 2022. a, b

van den Akker, J. J. H. and Hoogland, T.: Comparison of risk assessment methods to determine the subsoil compaction risk of agricultural soils in The Netherlands, Soil Till. Res., 114, 146–154,, 2011. a

van den Akker, J. J. H., de Vries, F., Vermeulen, G., Hack-ten Broeke, M., and Schouten, T.: Risico op ondergrondverdichting in het landelijk gebied in kaart, Alterra-rapport 2409, Alterra, Wageningen, the Netherlands, 2012. a

van den Berg, F., Tiktak, A., Hoogland, T., Poot, A., Boesten, J., van der Linden, A. M. A., and Pol, J. W.: An improved soil organic matter map for GeoPEARL_NL: Model description of version 4.4.4 and consequence for the Dutch decision tree on leaching to groundwater, Tech. rep., Wageningen Environmental Research (Alterra), Wageningen,, 2017. a, b, c, d, e, f

van den Elsen, E., van Tol-Leenders, D., Teuling, K., Römkens, P., de Haan, J., Korthals, G., and Reijneveld, A.: De staat van de Nederlandse landbouwbodems in 2018: Op basis van beschikbare landsdekkende dataset (CC-NL) en bodem-indicatorenlijst (BLN), Tech. Rep. 3048, Wageningen Environmental Research, Wageningen, (last access: 19 January 2024), 2020. a, b

van der Meulen, M., Doornenbal, J., Gunnink, J., Stafleu, J., Schokker, J., Vernes, R., van Geer, F., van Gessel, S., van Heteren, S., van Leeuwen, R., Bakker, M., Bogaard, P., Busschers, F., Griffioen, J., Gruijters, S., Kiden, P., Schroot, B., Simmelink, H., van Berkel, W., van der Krogt, R., Westerhoff, W., and van Daalen, T.: 3D geology in a 2D country: perspectives for geological surveying in the Netherlands, Netherlands J. Geosci., 92, 217–241,, 2013. a

van der Westhuizen, S., Heuvelink, G. B. M., Hofmeyr, D. P., and Poggio, L.: Measurement error-filtered machine learning in digital soil mapping, Spatial Statistics, 47, 100572,, 2022. a, b

van der Zee, S., van Riemsdijk, W., and de Haan, F.: Het protokol fosfaatverzadigde gronden I: toelichting, verslagen en mededelingen 1990-1A, Landbouwuniversiteit Wageningen, Landbouwuniversiteit Wageningen, Wageningen, (last access: 19 January 2024), 1990. a

van Doorn, M., Helfenstein, A., Ros, G. H., Heuvelink, G. B. M., van Rotterdam-Los, D. A. M. D., Verweij, S. E., and de Vries, W.: High-resolution digital soil mapping of amorphous iron- and aluminium-(hydr)oxides to guide sustainable phosphorus and carbon management, Geoderma, 443, 116838,, 2024. a

van Leeuwen, C. C. E., Mulder, V. L., Batjes, N. H., and Heuvelink, G. B. M.: Statistical modelling of measurement error in wet chemistry soil data, Eur. J. Soil Sci., 73,e13137,, 2021. a

van Tol-Leenders, D., Knotters, M., Groot, W. d., Gerritsen, P., Reijneveld, A., van Egmond, F., Wösten, H., and Kuikman, P.: Koolstofvoorraad in de bodem van Nederland (1998–2018): CC-NL, Rapport 2974, Wageningen Environmental Research, Wageningen,, 2019. a, b

Varón-Ramírez, V. M., Araujo-Carrillo, G. A., and Guevara Santamaría, M. A.: Colombian soil texture: building a spatial ensemble model, Earth Syst. Sci. Data, 14, 4719–4741,, 2022. a

Vasenev, V. I., Stoorvogel, J. J., Vasenev, I. I., and Valentini, R.: How to map soil organic carbon stocks in highly urbanized regions?, Geoderma, 226–227, 103–115,, 2014. a

Vasenev, V. I., Varentsov, M., Konstantinov, P., Romzaykina, O., Kanareykina, I., Dvornikov, Y., and Manukyan, V.: Projecting urban heat island effect on the spatial-temporal variation of microbial respiration in urban soils of Moscow megalopolis, Sci. Total Environ., 786, 147457,, 2021. a

Vaysse, K. and Lagacherie, P.: Using quantile regression forest to estimate uncertainty of digital soil mapping products, Geoderma, 291, 55–64,, 2017. a

Viscarra Rossel, R. A., Chen, C., Grundy, M. J., Searle, R., Clifford, D., and Campbell, P. H.: The Australian three-dimensional soil grid: Australia’s contribution to the GlobalSoilMap project, Soil Res., 53, 845,, 2015. a

Visschers, R., Finke, P. A., and de Gruijter, J. J.: A soil sampling program for the Netherlands, Geoderma, 139, 60–72,, 2007. a

von Hippel, P. T.: Mean, Median, and Skew: Correcting a Textbook Rule, J. Stat. Educ., 13, 2,, 2005. a

Vos, P.: Origin of the Dutch coastal landscape: long-term landscape evolution of the Netherlands during the Holocene, described and visualized in national, regional and local palaeogeographical map series, Barkhuis, Groningen, ISBN 978-94-91431-82-1, 2015. a

Vos, P., Meulen, M. V. D., Weerts, H., and Bazelmans, J.: Atlas of the Holocene Netherlands, landscape and habitation since the last ice age, Amsterdam University Press, Amsterdam, (last access: 17 March 2020), 2020. a

Wadoux, A. M. J. C.: Using deep learning for multivariate mapping of soil with quantified uncertainty, Geoderma, 351, 59–70,, 2019. a, b, c

Wadoux, A. M. J.-C. and Heuvelink, G. B. M.: Uncertainty of spatial averages and totals of natural resource maps, Methods Ecol. Evolut., 14, 1320–1332,, 2023. a

Wadoux, A. M. J.-C., Padarian, J., and Minasny, B.: Multi-source data integration for soil mapping using deep learning, SOIL, 5, 107–119,, 2019. a

Wadoux, A. M. J. C., Heuvelink, G. B. M., de Bruin, S., and Brus, D. J.: Spatial cross-validation is not the right way to evaluate map accuracy, Ecol. Model., 457, 109692,, 2021a. a, b

Wadoux, A. M. J. C., Heuvelink, G. B. M., Lark, R. M., Lagacherie, P., Bouma, J., Mulder, V. L., Libohova, Z., Yang, L., and McBratney, A. B.: Ten challenges for the future of pedometrics, Geoderma, 401, 115155,, 2021b. a, b, c, d

Walvoort, D. J. J. and de Gruijter, J. J.: Compositional Kriging: A Spatial Interpolation Method for Compositional Data, Math. Geol., 33, 951–966, 2001. a

Wamelink, G. W. W., Walvoort, D. J. J., Sanders, M. E., Meeuwsen, H. A. M., Wegman, R. M. A., Pouwels, R., and Knotters, M.: Prediction of soil pH patterns in nature areas on a national scale, Appl. Veg. Sci., 22, 189–199,, 2019. a

Wang, Z. and Shi, W.: Mapping soil particle-size fractions: A comparison of compositional kriging and log-ratio kriging, J. Hydrol., 546, 526–541,, 2017. a

WENR: Landelijk Grondgebruik Nederland (LGN), (last access: 18 January 2024), 2020. a, b, c

Withers, P. J. A., Sylvester-Bradley, R., Jones, D. L., Healey, J. R., and Talboys, P. J.: Feed the Crop Not the Soil: Rethinking Phosphorus Management in the Food Chain, Environ. Sci. Technol., 48, 6523–6530,, 2014.  a

Wright, M. N. and Ziegler, A.: ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R, J. Stat. Softw., 77, 1–17,, 2017. a

Short summary
Earth system models and decision support systems greatly benefit from high-resolution soil information with quantified accuracy. Here we introduce BIS-4D, a statistical modeling platform that predicts nine essential soil properties and their uncertainties at 25 m resolution in surface 2 m across the Netherlands. Using machine learning informed by up to 856 000 soil observations coupled with 366 spatially explicit environmental variables, prediction accuracy was the highest for clay, sand and pH.
Final-revised paper