Exploiting Soil and Remote Sensing Data Archives for 3D Mapping of Multiple Soil Properties at the Swiss National Scale

Abstract

Soils play a central role in ecosystem functioning, and thus, mapped soil property information is indispensable to supporting sustainable land management. Digital Soil Mapping (DSM) provides a framework to spatially estimate soil properties. However, broad-scale DSM remains challenging because of non-purposively sampled soil data, large data volumes for processing extensive soil covariates, and high model complexities due to spatially varying soil–landscape relationships. This study presents a three-dimensional DSM framework for Switzerland, targeting the soil properties of clay content (Clay), organic carbon content (SOC), pH value (pH), and potential cation exchange capacity (CEC_pot). The DSM approach is based on machine learning and a comprehensive exploitation of soil and remote sensing data archives. Quantile Regression Forest was applied to link the soil sample data from a national soil data base with covariates derived from a LiDAR-based elevation model, from climate raster data, and from multispectral raster time series based on satellite imagery. The covariate set comprises spatially multiscale terrain attributes, climate patterns and their temporal variation, temporarily multiscale land use features, and spectral bare soil signatures. Soil data and predictions were evaluated with respect to different landcovers and depth intervals. All reference soil data sets were found to be spatially clustered towards croplands, showing an increasing sample density from lower to upper depth intervals. According to the R² value derived from independent data, the overall model accuracy amounts to 0.69 for Clay, 0.64 for SOC, 0.76 for pH, and 0.72 for CEC_pot. Reduced model accuracies were found to be accompanied by soil data sets showing limited sample sizes (e.g., CEC_pot), uneven statistical distributions (e.g., SOC), and low spatial sample densities (e.g., woodland subsoils). Multiscale terrain covariates were highly influential for all models; climate covariates were particularly important for the Clay model; multiscale land use covariates showed enhanced importance for modeling pH; and bare soil reflectance was a major driver in the SOC and CEC_pot models.

1. Introduction

1.1. Soil Information Needs

Soils are an intrinsic component of terrestrial ecosystems, being formed and exposed at the intersections of the lithosphere, biosphere, atmosphere, and hydrosphere [1]. Therefore, soils determine the provision and regulation of numerous ecosystem services and functions related to the quality of air and water, food production, climate change, and biodiversity, among others [2,3]. The central role of soils for ecosystem functioning entails the need for extensive soil property information to support sustainable land management [4,5,6].

1.2. Digital Soil Mapping

Digital Soil Mapping (DSM) emerged over the last two decades as a labor- and cost-efficient, quantitative, reproducible, and updateable method to map soil properties [7]. DSM describes the creation of spatial soil information systems based on numerical models, linking soil observations with covariates to infer soil properties and associated uncertainties at unsampled locations [8,9,10]. Conceptually, DSM is rooted in Jenny’s clorpt model for soil formation, where soil at a specific location is described by climate, vegetation, relief, parent material, and time [11]. The DSM principle acknowledges the spatial variation of soils, where a soil property is described as a function of the so-called scorpan factors, thus covariates that are related to soil information (s), climate (c), vegetation (o), relief (r), parent material (p), time (a), and the geographic position (n) [12].

In the past decade, DSM has been increasingly applied at broader scales to address extensive and soil-related environmental issues, as well as to support local and regional soil sampling campaigns [13,14,15,16]. However, DSM at broader scales often requires coping with the spatially varying complexity and non-stationarity of the soil-landscape relationship, large data volumes, as well as sparse, unevenly distributed, and spatially clustered soil sample data [10,17].

Broad-scale DSM most commonly utilizes machine learning, a technique in data science that applies statistical and mathematical methods to correlate data and build models with potentially high computational efficiency [5,18,19]. Among the large group of machine learning algorithms, tree-based ensemble models, such as Quantile Regression Forest (QRF), a variant of Random Forest (RF), are most popular [20,21,22,23,24]. Both RF and QRF aggregate the results of numerous randomized regression trees for the final predictions. The randomization is based on a bootstrap sample for each tree and a random selection of covariates for each tree node. RF records the mean of the conditional distribution of the predicted target variable for all end nodes of a fitted tree. Instead, QRF preserves the entire conditional distribution for all end nodes, allowing us to infer a lower and upper quantile and, thus, a prediction interval (PI) for any location. In DSM, PIs are usually leveraged to describe the prediction uncertainty, indicated by the range in which the true value is expected to be related to a given probability. Several studies have approved QRF to produce plausible patterns of uncertainty variations [10,25]. Moreover, the success of the tree-based ensemble models can be attributed to further characteristics of the algorithms. The models are reported to be relatively resistant to overfitting, varying sample sizes, multicollinearity, and high data dimensionality [26,27,28]. Furthermore, no assumptions on data distributions are required; only a small number of hyperparameters need to be tuned, while the algorithms capture nonlinear relationships [29,30,31]. The combination of these benefits enables a robust representation of the soil–landscape relationship based on the scorpan concept, as well as efficient workflows with respect to data preprocessing, model calibration, and prediction [18]. Recently, tree-based ensemble models have been increasingly applied for broad-scale DSM [32,33,34]. Moreover, numerous DSM studies reported the performance superiority of these models compared to other machine learning algorithms, such as Support Vector Machines, Multiple Linear Regression, K-Nearest Neighbor, or Artificial Neural Networks [35,36,37].

According to a recent review by Chen et al. [5], broad-scale DSM studies of the past two decades predominantly applied 2D model designs, thus calibrating a model referring to a single soil depth layer only [38,39,40]. The recently prevalent 2.5D design acknowledges the depth dimension of soil by combining independent models that are fitted for different depth layers [41,42,43]. Most recently, an increasing number of DSM studies have applied 3D model designs that incorporate the depth dimension more explicitly by spatially fitting depth functions or integrating depth as a covariate [44,45,46,47]. The 3D design, which applies depth as a covariate, is associated with a straightforward and flexible workflow, allowing predictions at any location and depth based on a single calibrated model. Moreover, the calibration sample size is maximized, while complex assumptions about depth trends are avoided [33,34,48]. However, Nauman and Duniway (2019) [49] stated that these 3D designs might lead to increased prediction uncertainties for soil properties with an enhanced depth gradient. Ma et al. [50] found unrealistically stepped depth functions. Nevertheless, model performances between 2.5D and 3D designs are reported to be generally comparable, while any design preference might be case-specific, depending on the soil property and sample size [5,32,50].

Besides the advent of machine learning DSM, the increasing availability of air- and spaceborne remote sensing data as well as high-performance computing facilities has initiated major progress for broad-scale DSM [51,52,53]. Particularly, airborne LiDAR-based digital elevation models (DEM) in high spatial resolution and multispectral satellite imagery in high spatio-temporal resolution, combined with methodological advancements, allow for the extensive and more accurate processing of scorpan covariates [54,55,56].

1.3. Terrain Covariates for DSM

Generally, the scorpan factors determine soil formation by describing soil-forming processes that interact and operate across multiple spatial scales [57,58,59]. In this context, Behrens et al. [60] emphasized that spatially multiscale terrain attributes from high-resolution DEMs represent proxies for multiple scorpan factors related to gravimetry, hydrology, climate, land use, and lithology. This has been acknowledged by several DSM studies, proposing methods to process multiscale terrain covariates based on resampling techniques, expanding convolution kernels, or wavelet transformations [60,61,62,63]. Although the aforementioned studies revealed improved model performances, the majority of DSM studies make use of terrain covariates at a single spatial scale only, and thus, soil-forming processes operating on a different scale might be inadequately described [64,65,66].

1.4. Multispectral Covariates for DSM

Recently, long-term multispectral raster time series, derived from imagery archives of the Landsat (LA) and Sentinel-2 (S2) satellite missions, have been leveraged to process scorpan covariates related to soil surface and land use characteristics [67,68,69].

Referring to soil surface characteristics, the spectral signature of bare soil surfaces correlates with mineralogical, physical, and chemical soil properties [70,71,72]. These bare soil covariates are approved to be expedient for DSM, particularly concerning topsoil clay content, pH, organic carbon, and cation exchange capacity [40,56,73]. The processing of bare soil covariates comprises a filtering of the multispectral time series data for pixels that represent bare soil surfaces and a subsequent aggregation of the filtered raster stack across time [38,74,75].

Referring to land use characteristics, the spectral signature of vegetation, as a land use indicator, correlates with biological, chemical, physical, and hydraulic soil properties [20,76,77]. Related scorpan covariates can be processed as averages across time based on vegetation indices derived from the multispectral raster time series [28,78,79]. These land use covariates describe temporarily static vegetation conditions for one defined period, are commonly applied for DSM, and are being approved to be particularly relevant for mapping organic carbon, pH, carbonates, texture, and cation exchange capacity [19,80,81]. However, vegetation characteristics tend to change over time and across multiple temporal scales [82,83,84]. For example, the temporal variability of vegetation differs across seasonal, inter-annual, or inter-decadal time scales, and thus, it is related to different and interacting soil-forming processes [54,85]. Similar to multiscale terrain analysis, the vast majority of DSM studies disregard this temporal multiscale character of land use by using temporarily static land use covariates based on a single time scale only. However, Heuvelink et al. [82] applied DSM to map soil organic carbon in space and time, using temporarily multiscale land use covariates based on a seasonal vegetation index processed for each year and retrospectively across multiple time scales, defined by an exponential decay function.

1.5. Study Outline and Objectives

This study provides a broad-scale 3D DSM framework for Switzerland, addressing all crop-, grass-, and woodland areas (‘Cropland’, ‘Grassland’, and ‘Woodland’) and three depth intervals (‘Topsoil’: 0–30 cm; ‘Subsoil’: 30–60 cm; and ‘Deep Subsoil’: 60–120 cm). The target soil properties are clay content (Clay, %), soil organic carbon concentration (SOC, %), pH value (pH, CaCl₂), and potential cation exchange capacity (CEC_pot, mm_ckg⁻¹), which were derived from archived soil sample data sets. Scorpan covariates were derived by exploiting remote sensing data archives, thus a LiDAR-based DEM, two multispectral raster time series based on LA and S2 imagery, respectively, as well as climate raster data. The covariate set comprises spatially multiscale terrain attributes, temporarily multiscale land use features, spectral bare soil signatures, as well as climate patterns and their temporal variation. We applied the machine learning algorithm QRF for a 3D model design, linking each soil sample set with the scorpan covariates and using soil depth as an additional covariate. Hence, the introduced framework for broad-scale DSM combines recent advancements related to the engineering of scorpan covariates with a state-of-the-art 3D model design. The latter allows us to efficiently derive map predictions and associated uncertainties at any location and depth based on a single calibrated model. Moreover, the evaluation scheme of the applied 3D model design allows for the derivation of model accuracies across landcovers and depth intervals based on independent soil data. The specific objectives of the study are defined by (i) analyzing soil sample data sets with respect to spatial sample density and statistical distributions’ (ii) evaluating predicted maps for Clay, SOC, pH, and CEC_pot, referring to landcover and depth intervals; and (iii) assessing the influence of covariates related to terrain, bare soil, land use, and climate on model performances.

2. Materials and Methods

2.1. Study Area

Switzerland is a landlocked country in the center of Europe, covering 41,291 km². Altitudes range from 196 m to 4634 m (a.s.l.). Annual precipitation ranges from 530 mm to 3300 mm, while the mean temperature in January amounts to 1 °C in January and 17 °C in July [86,87,88]. With an areal proportion of 45%, the largest proportion of the Swiss territory is devoted to agriculture. This encompasses croplands (9%; ‘Cropland’) and grasslands (36%; ‘Grassland’), which include permanent grasslands (14%), and extensively managed alpine grasslands (22%). Forest areas (‘Woodland’) and settlements cover 28% and 8% of the study area, respectively. Approximately 19% of the Swiss territory is unproductive due to rock, water, and ice surfaces. The study area comprises ‘Cropland’, ‘Grassland’, and ‘Woodland’, covering approximately 30,394 km² (Figure 1) [89].

2.2. Soil Model Design

We used QRF to link soil sample data for Clay, SOC, pH, and CEC_pot with scorpan covariates. The final prediction was derived based on the median of the conditional distribution recorded for each end node of a fitted tree. The prediction uncertainty is based on the QRF prediction intervals, defined by the 5%- and 95%-quantiles of the conditional distribution. The final prediction uncertainty (RPI) was then defined by the ratio between the PI and the median of the conditional distribution (cf. 1.2).

The subset of sample data, which is excluded for the calibration of a single QRF tree, allows us to derive the model accuracy as the average across all trees (OOB_error). The hyperparameters ‘mtry‘ and ‘numtrees‘ define the size of the covariate subset at each node and the number of trees in the ensemble, respectively. These hyperparameters were determined by grid tuning based on the OOB_error and using the Coefficient of Determination (R²; Equation (1)) as an accuracy indicator [90]. The R² is a measure of covariation between observed and predicted values, as follows:

$$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{{(y}_i^{observed}-y_i^{predicted})}^2}{{{\sum }_{i=1}^n\left(y_i^{observed}-\overline{y^{observed}}\right)}^2}}$$

(1)

The scorpan covariates represent extensively available raster data sets related to terrain, climate, land use, and the spectral signature of bare soil. All covariates were resampled to a common raster cell size of 30 × 30 m and assigned to the soil sample data sets. Covariates related to terrain, climate, and bare soil were assumed to be static in time and assigned to the soil data based on the geographic position only. Covariates for land use describe temporarily dynamic phenomena and are annually available. These covariates were assigned to the soil data based on the geographic position and the year of soil data acquisition.

The soil depth served as an additional covariate and was primarily harmonized to the standard depth intervals of ‘Topsoil’ (0–30 cm), ‘Subsoil’ (30–60 cm), and ‘Deep Subsoil’ (60–120 cm). The depth harmonization is based on fitting an equal-area quadratic spline down the soil profile of each sample location and soil property. The soil property value within each depth interval was defined by the mean of the fitted spline [91,92]. Since the fitting of a spline requires more than one data value, sample values at locations with only one sampled depth were transferred to their associated standard depth interval [93].

The accuracy evaluation is based on five-fold stratified random data splitting. Thus, the soil data sets were stratified by the depth intervals and landcover classes ‘Cropland’, ‘Grassland’, and ‘Woodland’. Within each stratum, the soil data sets were randomly split by a ratio of 80/20 to form a calibration set and an independent validation set not used in the model (IV) for all strata. The R² (Equation (1)) and the Root Mean Square Error (RMSE; Equation (2)), absolute measures of the average error, were used as accuracy indicators.

$$RMSE=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^n{(y_i^{predicted}-y_i^{observed})}^2}{n}}}$$

(2)

The R² and RMSE were derived using the OOB_error, a 10-fold cross-validation approach (CV10), and IV. Moreover, both accuracy indicators and, additionally, the uncertainty indicator RPI were derived for each stratum of the IV. Based on different data splits, the modeling was repeated fivefold to generate independent submodels. The overall model accuracy and uncertainty were defined by the mean of the respective indicators across all submodels. Additionally, we evaluated the model robustness by analyzing the standard deviation for each accuracy indicator across all submodels (SD_R2; SD_RMSE), as well as by comparing the results of the OOB_error, CV10, and IV.

The best submodel was defined by the closest distance to the mean across all submodels, referring to R² and RMSE based on the IV. This model was then used to derive spatial predictions of the soil properties and their respective RPI for all depth intervals using the most recent temporarily dynamic covariates.

A covariate importance analysis was conducted to assess the influence of the covariates on the model performance. Covariates related to terrain, climate, land use, and bare soil were each defined as a group. In the next step, the performance of additional calibrated models, for which all covariates of a particular group were randomly permuted, was compared to the performance of the original models, which were calibrated with the original data [94]. The decline in model performance therefore indicates the importance of the particular covariate group. For this specific analysis, the model performance was assessed by R² based on the OOB_error. The covariate importance was derived for each covariate group and submodel, while the final covariate importance was defined by the mean across all submodels.

2.3. Soil Sample Data

The soil sample data for Clay (%), SOC (%), pH (CaCl₂), and CEC_pot (mm_ckg⁻¹) stem from two sources. First, the ‘National Soil Information System NABODAT’ collects, merges, and maintains soil profile data from various soil surveys at regional scales across Switzerland since the 1960s. We used version 6 of the ‘National Soil Dataset’ [95]. Second, the ‘Swiss Biodiversity Monitoring’ (BDM), which includes a soil survey conducted between 2011 and 2015 based on a regular grid network of 6 by 4 km spanning across Switzerland. At around 1200 sampling plots, 2–4 topsoil samples (0–20 cm) were collected from the outer rim in divergent directions from a center point and finally averaged [96,97,98].

The soil sample data represent a conglomerate of the most diverse soil surveys of the past six decades and show increased heterogeneity with respect to underlying methodological and technical standards. This requires several data harmonization and cleansing steps. Hence, soil property labels and units were standardized, while data from surveys prior to 1985 were excluded. Profile data with overlapping soil horizons, duplicates, and implausible values were identified and removed. Data associated with sealed surfaces, according to current land cover models and referring to the time of data acquisition, indicate uncertain georeferencing and were removed. With respect to all soil properties, the soil sample data comprise 17,283 unique sample locations (n_xy).

Clay sample data were analyzed according to the sedimentation method [99]. The majority of the SOC data (81%) were analyzed using the Swiss standard method based on wet oxidation and retitration of potassium dichromate [100]. Moreover, 19% of the SOC sample data were analyzed by dry combustion using a CN-analyzer [101]. Since SOC values determined by wet oxidation are generally lower compared to an analysis by dry combustion, we adapted these samples by applying a conversion factor of 1.18 [102]. Moreover, SOC values from organic soils are inherently higher compared to SOC values from mineral soils, while available data prevent a definite distinction. Thus, we removed SOC data with values exceeding three times the interquartile range from the mean [79,103]. Sample data for pH were measured potentiometrically in a 0.01 M CaCl₂ solution [104]. The majority of the CEC_pot data (89%) were analyzed by extracting exchangeable cations (Ca, Mg, K, and Na⁺) in a 0.1 M BaCl₂ or a 0.1 M HCI + H₂SO₄ solution. Elemental concentrations were measured by atomic absorption spectroscopy (AAS) or a flame photometer. The summarized exchangeable charges of Ca, Mg, K, Na, and H correspond to CEC_pot [104]. Exchangeable cations of fewer CEC_pot data (11%) were extracted in a 1 M NH₄Cl solution, including for Al and Fe, while atomic emission spectroscopy (ICP-AES) was applied to measure elemental concentrations [105,106,107].

Finally, the soil sample data were analyzed according to their spatial sample density (SSD; n_xy*km⁻²) and statistical distributions with respect to landcover classes and depth intervals.

2.4. Terrain Covariates

The data source for terrain covariates is the digital elevation model (DEM) swissALTI^3D with a spatial resolution of 2 × 2 m. For areas below 2000 m a.s.l., swissALTI^3D is based on airborne LiDAR measurements and shows an accuracy of 0.5 m with respect to position and elevation. For areas above 2000 m a.s.l., the swissALTI^3D is predominantly based on stereo orthophoto correlation of overlapping aerial images, showing an accuracy of 1–3 m [88].

We conducted a multiscale terrain analysis, where terrain attributes were derived at different spatial scales using mixed scaling [108]. Therefore, the DEM was rescaled to multiple versions of increasingly coarser resolutions. Then, terrain attributes were derived from each of the rescaled DEMs. In a final step, the terrain attributes were rescaled to the initial resolution of the DEM (Figure S1) [60,108,109].

We spatially extended the swissALTI^3D data by 30 km across the Swiss external border to reduce artefacts that would arise from processing terrain attributes close to the border. For that purpose, we used elevation data from the Shuttle Topography Radar Mission (SRTM), available at a spatial resolution of 30 × 30 m [110]. The swissALTI^3D and the SRTM were resampled to a consistent spatial resolution of 8 × 8 m using bilinear interpolation and merged to form the pivotal database for the multiscale terrain analysis.

In total, we derived 11 terrain attributes (

Table 1. List of terrain attributes used for the multiscale terrain analysis.

Terrain Attribute	Reference
Elevation (a.s.l.) [m]	[88]
Slope Inclination [°]	[112]
Terrain Ruggedness Index [-]	[113]
Terrain Roughness Index [-]	[114]
Topographic Position Index [-]	[115]
Profile Curvature [-]	[112]
Planform Curvature [-]	[112]
Flow Accumulation, log10 [-]	[116]
Convergence Index [-]	[117]
Eastness [-], sin (aspect in rad)	[114]
Northness [-], cos (aspect in rad)	[114]

) at 14 different spatial scales between 8 m and 4880 m. The selected terrain attributes were approved to be relevant for the target soil properties [15,50,111]. The spatial scales represent the initial resolution of the DEM factorized by the Fibonacci numbers n (n = 1, 2, …, 610). In total, the multiscale terrain analysis resulted in 154 terrain covariates.

2.5. Climate Covariates

Climate covariates are based on annual raster data for precipitation totals, mean relative sunshine duration, as well as means of daily average, maximum, and minimum temperatures. These data sets are available at a spatial resolution of 1000 × 1000 m for the period from 1980 to 2021 (

Table 2. List of annually available (1980–2021) climate grid data used to derive long-term means and mean interannual changes.

Climate Raster	Reference
Precipitation totals [mm]	[118]
Temperature daily average [°C]	[119]
Temperature daily maximum [°C]	[119]
Temperature daily minimum [°C]	[119]
Relative sunshine duration [%]	[120]

) [87]. According to Frei et al. [118,119,120], the precipitation and temperature grids were interpolated based on a network of 430 rain gauges and 125 temperature stations across Switzerland. The relative sunshine duration represents the ratio between the effective sunshine and the maximally possible sunshine duration, assuming that no clouds are present. The grids for relative sunshine duration are based on 70 measurement stations and a clearness index derived from Meteosat satellite imagery. We derived climate covariates as the mean and the mean interannual change in each climate parameter across 1980 to 2021 (Figure S2). In total, the analysis results in 10 individual climate covariates.

2.6. Spectral Covariates

2.6.1. Multispectral Raster Time Series

The data source for covariates related to land use consists of two raster time series based on LA and S2 data, respectively. The data source for covariates related to the spectral signature of bare soil is based only on the LA raster time series [121,122].

The LA raster time series was compiled from all available scenes of all LA sensors (TM, ETM+, and OLI) for the period from 1985 to 2021. The LA sensors record a blue (BLUE), green (GREEN), red (RED), and near-infrared (NIR) band, as well as two shortwave-infrared (SWIR1, SWIR2) bands at a spatial resolution of 30 m × 30 m. LA scenes of the sensors TM and ETM+ were readily corrected to surface reflectance according to the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [123]. LA scenes of the OLI sensor were corrected according to the Land Surface Reflectance Code (LaSRC) [124]. All pixels covered by clouds and cloud shadows were masked using the CFmask algorithm [125]. Since there are differences in the reflective wavelengths of the LA sensors TM/ETM+ and OLI, we harmonized the data using the RMA (Reduced Major Axis) regression coefficients provided by Roy et al. [126].

The spectral vegetation indexes Normalized Difference Vegetation Index (NDVI; Equation (3)) and Normalized Burn Ratio 2 (NBR2; Equation (4)) were derived as additional bands for each scene of the LA raster time series [127,128]. The NDVI describes a proxy for plant condition and scales between −1 and 1. Negative NDVI values correspond to water, low positive values to nearly non-vegetated areas, and increasing values to enhanced biomass and healthier vegetation.

$$NDVI={\displaystyle \frac{NIR-RED}{NIR+RED}}$$

(3)

The NBR2 is based on the ratio of two shortwave-infrared bands (SWIR1 and SWIR2) and serves as a proxy for litter deposit, soil moisture, and soil texture, among others [67]. The NBR2 index also scales between −1 and 1, while positive values correspond to vegetation and negative values indicate bare surfaces [129].

$$NBR2={\displaystyle \frac{SWIR1-SWIR2}{SWIR1+SWIR2}}$$

(4)

The S2 raster time series was compiled from all available scenes of the S2 sensors for the period from 2017 to 2021. The S2 sensors record a blue, green, red, and near-infrared band at a spatial resolution of 10 m × 10 m. S2 scenes were corrected to surface reflectance using the sen2cor algorithm [130,131]. All pixels covered by clouds and cloud shadows were masked using the S2cloudless algorithm [132,133]. Moreover, we derived the NDVI as an additional band for each scene of the S2 raster time series (Equation (3)).

Both the LA and S2 raster time series were temporarily restricted to scenes acquired from the 1st of March to the 31st of October of each year. Also, the raster time series were spatially limited to areas of existing soil surfaces, and thus, sealed surfaces as well as areas covered by rock, snow, ice, and water were masked. The raster time series were obtained from the Google Earth Engine Data Catalogue and processed using the Google Earth Engine API [134,135].

2.6.2. Land Use Covariates

The land use covariates are based on multiple annual aggregations of the LA and S2 NDVI raster time series using the median, maximum, minimum, and standard deviation (

Table 3. List of NDVI raster time series based on annually aggregated data used for the land use trend analysis (NDVI: Normalized Difference Vegetation Index; S2: Sentinel-2; LA: Landsat).

Data	Annual Aggregation	Timespan	Resolution
S2	Median	2017–2022	10 m
LA		1985–2022	30 m
S2	Maximum	2017–2022	10 m
LA		1985–2022	30 m
S2	Minimum	2017–2022	10 m
LA		1985–2022	30 m
S2	Standard deviation	2017–2022	10 m
LA		1985–2022	30 m

). For each aggregation, we computed multiple statistical trend functions retrospectively for each year of the time series and across multiple temporal scales. We defined the trend functions by the mean, standard deviation, and the Theil-Sen trend estimator [136,137]. The temporal scales were defined by five-year intervals for the LA and by annual intervals for the S2 raster time series. Since the LA and S2 raster time series span across 37 and 5 years, respectively, the number of temporal scales depends on the referred year of the time series. In total, we derived 1558 covariates related to land use (Figures S3 and S4).

2.6.3. Spectral Bare Soil Covariates

The spectral bare soil covariates represent band-wise composites based on a median aggregation of the LA raster time series, using pixels that flag bare surfaces only (

Table 4. List of bare soil composites and their spectral range used to derive bare soil covariates.

Bare Soil Composite	Spectral Range [µm]
BS_BLUE	0.45–0.52
BS_GREEN	0.52–0.60
BS_RED	0.63–0.69
BS_NIR	0.76–0.90
BS_SWIR1	1.55–1.75
BS_SWIR2	2.08–2.35

). Pixels showing NDVI values > 0.25, which implies green vegetation and growing crops, were removed [40,73]. Pixels with NBR2 values > 0.075 were also excluded since these pixels indicate straw, crop residues, or biased spectral signatures due to high soil moisture content [38,75]. Pixel stacks with a bare soil pixel frequency < 3 were rejected to prevent biased median aggregations due to the limited data population. Moreover, pixel stacks showing an elevation > 1200 m (a.s.l.) and a slope > 25° were excluded since topographic roughness combined with low solar angles causes biased spectral signatures due to shading effects [138,139].

The bare soil composites cover the study area only partially due to the existence of permanently vegetated areas (e.g., forests). These data gaps were filled for each bare soil composite using the spatial modeling approach of the aforementioned soil model design (cf. 2.2) with covariates related to terrain, climate, and land use (cf. 2.4–2.6). For this purpose, we extracted calibration sample sets from the bare soil composites, representing their geographical and spectral space evenly. Therefore, we derived spatial clusters using unsupervised K-means based on the geographic position and the spectral bare soil signature. Within these clusters, we randomly selected subsamples of size n = 5000, which we combined to form the calibration sample set. The optimal number of clusters was determined using the Silhouette Method (Figure S5) [140].

3. Results

3.1. Soil Sample Data

The soil sample data were mainly acquired in the 1990s and generally cover the entire study area (Figure S6). The spatial sample density (SSD, n_xy*km⁻²) amounts to 0.34 for Clay, 0.37 for SOC, 0.56 for pH, and 0.06 for CEC_pot. All sample data sets show a generally decreasing SSD from ‘Cropland’ to ‘Grassland’ and ‘Woodland’, as well as from the ‘Topsoil’ to ‘Subsoil’ and ‘Deep Subsoil’. The SSDs of ‘Cropland’ are distinctively larger by approximately 80% compared to ‘Grassland’ and ‘Woodland’. The declines of SSD from the upper to lower depth intervals are more continuous, ranging between 23% and 48% (

Table 5. Spatial sample density SSD (n_xy*km⁻²; n_xy: Count of unique sample locations) for the soil sample data (Clay, SOC, pH, and CEC_pot). Provided values refer to the entire study area (‘Total area’) across all depth intervals, as well as to specific depth intervals (‘Topsoil’, ‘Subsoil’, and ‘Deep Subsoil’) and landcover classes (‘Cropland’, ‘Grassland’, and ‘Woodland’).

Soil Property	Depth	Total Area	Cropland	Grassland	Woodland
Clay	0–120 cm	0.34	1.20	0.22	0.20
	0–30 cm	0.34	1.19	0.22	0.20
	30–60 cm	0.26	0.96	0.16	0.15
	60–120 cm	0.20	0.73	0.12	0.12
SOC	0–120 cm	0.37	1.25	0.23	0.25
	0–30 cm	0.37	1.24	0.22	0.24
	30–60 cm	0.20	0.62	0.11	0.17
	60–120 cm	0.12	0.33	0.06	0.12
pH	0–120 cm	0.56	1.84	0.36	0.37
	0–30 cm	0.56	1.83	0.35	0.37
	30–60 cm	0.29	0.97	0.17	0.19
	60–120 cm	0.21	0.72	0.12	0.15
CECpot	0–120 cm	0.06	0.21	0.04	0.04
	0–30 cm	0.06	0.20	0.04	0.04
	30–60 cm	0.04	0.15	0.02	0.02
	60–120 cm	0.03	0.13	0.02	0.01

).

The count of unique sample locations (n_xy) for the soil sample data amount to 10,467 for Clay, 11,329 for SOC, 17,052 for pH, and 1879 for CEC_pot. The descriptive statistics of the Clay (%) data show a marginally right-skewed distribution with a median of 22.04 and a standard deviation of 10.78. The SOC (%) data show an explicitly right-skewed distribution with a median and standard deviation of 1.75 and 1.76, respectively. The pH (CaCl₂) data set is left-skewed, with a median of 6.46 and a standard deviation of 1.14. The CEC_pot (mm_ckg⁻¹) data set is right-skewed with a median and standard deviation of 150.48 and 88.43, respectively (

Table 6. Descriptive statistics of the soil sample data (Clay; SOC; pH; CEC_pot; Min: Minimum; Max: Maximum; SD: Standard deviation; n_xy: Count of unique sample locations).

Soil Property	Min	Max	Mean	Median	SD	nxy
Clay [%]	0.20	69.17	24.16	22.04	10.78	10,476
SOC [%]	0.01	10.19	2.14	1.75	1.76	11,329
pH [-]	2.50	9.00	6.24	6.46	1.14	17,052
CECpot [mmckg−1]	1.04	554.35	170.46	150.48	88.43	1879

).

The frequency distributions of the soil sample data for Clay (%), SOC (%), and CEC_pot (mm_ckg⁻¹), stratified by landcover classes and depth intervals, are generally similar across the strata. This similarity is confirmed for the pH (CaCl₂) sample data across depth intervals, as well as across ‘Cropland’ and ‘Grassland’, while showing a left-skewed distribution. However, the pH (CaCl₂) data set for ‘Woodland’ shows a bimodal distribution across all depth intervals. Apart from that, the distributions of the stratified data sets are comparable to those of the entire data population (Figure 2).

3.2. Clay Model

The Clay model shows an R² of 0.69, an RMSE of 6.07, and an RPI of 0.89 as the mean across five independent submodels and based on the IV. With respect to depth intervals, the model achieved the highest accuracy and lowest uncertainty for ‘Subsoil’ (R² = 0.76; RMSE = 5.39; RPI = 0.84), followed by ‘Topsoil’ (R² = 0.69; RMSE = 5.89; RPI = 0.91) and ‘Deep Subsoil’ (R² = 0.62; RMSE = 7.13; RPI = 0.92). With respect to landcover classes and according to the R², the best results were obtained for ‘Grassland’ (R² = 0.76), followed by a similar R² for ‘Cropland’ and ‘Woodland’ (R² = 0.68). According to the RMSE and RPI, the results for ‘Cropland’ (RMSE = 5.98; RPI = 0.77) outperform ‘Grassland’ (RMSE = 6.03; RPI = 0.86), with a further performance decrease for ‘Woodland’ (RMSE = 6.32; RPI = 1.03;

Table 7. Performance of the Clay model according to accuracy (R² and RMSE) and uncertainty (RPI) based on independent validation data. Provided values are averaged across five independent submodels and refer to the entire study area (‘Total area’) across all depth intervals, as well as to specific depth intervals (‘Topsoil’, ‘Subsoil’, and ‘Deep Subsoil’) and landcover classes (‘Cropland’, ‘Grassland’, and ‘Woodland’).

Indicator	Depth	Total Area	Cropland	Grassland	Woodland
R2	0–120 cm	0.69	0.68	0.76	0.68
	0–30 cm	0.69	0.68	0.73	0.62
	30–60 cm	0.76	0.75	0.77	0.77
	60–120 cm	0.62	0.61	0.63	0.64
RMSE	0–120 cm	6.07	5.98	6.03	6.32
	0–30 cm	5.89	5.58	5.91	6.51
	30–60 cm	5.39	5.35	5.40	5.44
	60–120 cm	7.13	7.27	6.97	6.95
RPI	0–120 cm	0.89	0.77	0.86	1.03
	0–30 cm	0.91	0.80	0.86	1.06
	30–60 cm	0.84	0.71	0.84	0.97
	60–120 cm	0.92	0.81	0.89	1.07

; Figure 3 and Figure S7).

The Clay model’s performance is robust, showing low variability across the submodels with a SD_R2 of 0.01 and a SD_RMSE of 0.02. The accuracy indications based on the IV are confirmed, referring to the OOB_error (R² = 0.69; SD_R2 = 0.01; RMSE = 5.99; SD_RMSE = 0.02) and the CV10 (R² = 0.68; SD_R2 = 0.01; RMSE = 6.13; SD_RMSE = 0.02).

The covariate importance analysis revealed a high influence of covariates related to terrain and climate. The impact of spectral bare soil and land use covariates on the model performance was intermediate and low, respectively (Figure 3).

The Clay mapping shows consistent patterns across landcover classes and depth intervals. The mean Clay values range between 20% and 25% for all strata. However, referring to ‘Grassland’ and ‘Woodland’, Clay values are marginally larger for the ‘Subsoil’ (Figure 3, Figure 4 and Figure S8).

3.3. SOC Model

The SOC model shows an R² of 0.64, an RMSE of 1.06, and an RPI of 3.04 as the mean across five independent submodels and based on the IV. With respect to depth intervals as well as according to R² and RPI, the model performs best for ‘Topsoil’ (R² = 0.50; RPI = 1.57), followed by ‘Subsoil’ (R² = 0.45; RPI = 3.71) and ‘Deep Subsoil’ (R² = 0.37; RPI = 3.85). According to RMSE, the best model performance is indicated for ‘Subsoil’ and ‘Deep Subsoil’ (RMSE = 0.92). With respect to landcover classes and according to R², the model shows the highest accuracy for ‘Grassland’ (R² = 0.66), followed by ‘Woodland’ (R² = 0.62) and ‘Cropland’ (R² = 0.61). According to RMSE and RPI, the model performs best on ‘Cropland’ (RMSE = 0.94; RPI = 2.83), followed by ‘Grassland’ (RMSE = 1.09; RPI = 2.87) and ‘Woodland’ (RMSE = 1.21; RPI = 3.43;

Table 8. Performance of the SOC model according to accuracy (R² and RMSE) and uncertainty (RPI) based on independent validation data. Provided values are averaged across five independent submodels and refer to the entire study area (‘Total area’) across all depth intervals, as well as to specific depth intervals (‘Topsoil’, ‘Subsoil’, and ‘Deep Subsoil’) and landcover classes (‘Cropland’, ‘Grassland’, and ‘Woodland’).

Indicator	Depth	Total Area	Cropland	Grassland	Woodland
R2	0–120 cm	0.64	0.61	0.66	0.62
	0–30 cm	0.50	0.51	0.48	0.37
	30–60 cm	0.45	0.46	0.42	0.47
	60–120 cm	0.37	0.37	0.38	0.41
RMSE	0–120 cm	1.06	0.94	1.09	1.21
	0–30 cm	1.16	0.96	1.13	1.5
	30–60 cm	0.92	0.91	1.01	0.83
	60–120 cm	0.92	0.84	1.01	0.90
RPI	0–120 cm	3.04	2.83	2.87	3.43
	0–30 cm	1.57	1.49	1.47	1.75
	30–60 cm	3.71	3.33	3.67	4.12
	60–120 cm	3.85	3.67	3.47	4.42

; Figure 5 and Figure S7).

The SOC model performance is robust, as indicated by the low variability across the submodels with a SD_R2 and a SD_RMSE of 0.01. The model accuracy indications based on the IV are comparable to those derived from OOB_error (R² = 0.62; SD_R2 = 0.01; RMSE = 1.08; SD_RMSE = 0.01) and CV10 (R² = 0.63; SD_R2 = 0.01; RMSE = 1.10; SD_RMSE = 0.01).

The covariate importance analysis indicated that covariates related to spectral bare soil and terrain are highly influential for the SOC model. Covariates related to climate and land use are likewise less important (Figure 5).

With respect to predicted SOC values for landcover classes and depth intervals, the mean SOC decreases significantly from ‘Topsoil’ to ‘Subsoil’ by 89% and is less pronounced (59%) from ‘Subsoil’ to ‘Deep Subsoil’. Moreover, the mean SOC decreased for ‘Cropland’ (1.25%) compared to ‘Grassland’ (1.88%) and ‘Woodland’ (1.81%; Figure 5, Figure 6 and Figure S9).

3.4. pH Model

The pH model shows an R² of 0.76, an RMSE of 0.56, and an RPI of 0.39 as the mean across five independent submodels and based on the IV. With respect to depth intervals as well as according to R² and RMSE, the highest accuracy was found for ‘Subsoil’ (R² = 0.86; RMSE = 0.45), followed by ‘Deep Subsoil’ (R² = 0.8; RMSE = 0.57) and ‘Topsoil’ (R² = 0.70; RMSE = 0.60). The lowest uncertainty was indicated for ‘Subsoil’ (RPI = 0.36), while ‘Topsoil’ (RPI = 0.39) revealed a decreased uncertainty compared to ‘Deep Subsoil’ (RPI = 0.43). With respect to landcover classes, R² is comparable between ‘Cropland’ (R² = 0.72), ‘Grassland’ (R² = 0.72), and ‘Woodland’ (R² = 0.71). However, according to RMSE and RPI, the model performance is best for ‘Cropland’ (RMSE = 0.43; RPI = 0.25), followed by ‘Grassland’ (RMSE = 0.51; RPI = 0.39) and ‘Woodland’ (RMSE = 0.79; RPI = 0.53;

Table 9. Performance of the pH model according to accuracy (R² and RMSE) and uncertainty (RPI) based on independent validation data. Provided values are averaged across five independent submodels and refer to the entire study area (‘Total area’) across all depth intervals, as well as to specific depth intervals (‘Topsoil’, ‘Subsoil’, and ‘Deep Subsoil’) and landcover classes (‘Cropland’, ‘Grassland’, and ‘Woodland’).

Indicator	Depth	Total Area	Cropland	Grassland	Woodland
R2	0–120 cm	0.76	0.72	0.72	0.71
	0–30 cm	0.70	0.61	0.66	0.63
	30–60 cm	0.86	0.85	0.83	0.84
	60–120 cm	0.80	0.79	0.75	0.81
RMSE	0–120 cm	0.56	0.43	0.51	0.79
	0–30 cm	0.60	0.47	0.53	0.88
	30–60 cm	0.45	0.34	0.43	0.61
	60–120 cm	0.57	0.44	0.56	0.77
RPI	0–120 cm	0.39	0.25	0.39	0.53
	0–30 cm	0.39	0.26	0.34	0.56
	30–60 cm	0.36	0.24	0.32	0.51
	60–120 cm	0.43	0.25	0.50	0.53

; Figure 7 and Figure S7).

The pH model performance is robust with low variability across five independent submodels, as indicated by a SD_R2 and a SD_RMSE of 0.01. Moreover, accuracy levels based on the IV are similar to those derived from OOB_error (R² = 0.75; SD_R2 = 0.01; RMSE = 0.56; SD_RMSE = 0.01) and CV10 (R² = 0.74; SD_R2 = 0.01; RMSE = 0.58; SD_RMSE = 0.01).

The covariates related to terrain and land use are indicated as highly influential for the performance of the pH model. The impact of covariates related to climate and spectral bare soil is less pronounced (Figure 7).

With respect to map values, the mean pH increases from ‘Topsoil’ to ‘Subsoil’ and ‘Deep Subsoil’ between 1% and 8% across landcover classes. Generally, the mean pH is larger for ‘Cropland’ (6.57) compared to ‘Grassland’ (5.87) and ‘Woodland’ (5.71; Figure 7, Figure 8 and Figure S10).

3.5. CEC_pot Model

The CEC_pot model shows an R² of 0.72, an RMSE of 48.61, and an RPI of 1.31 as the mean across five independent submodels and based on the IV. With respect to depth intervals and according to R², the model performs best for ‘Topsoil’ (R² = 0.71), followed by ‘Subsoil’ (R² = 0.67) and ‘Deep Subsoil’ (R² = 0.59). According to RMSE, the highest accuracy was found for ‘Subsoil’ (RMSE = 44.61), followed by ‘Deep Subsoil’ (RMSE = 45.08) and ‘Topsoil’ (RMSE = 50.65). The RPI is lowest for ‘Topsoil’ (RPI = 1.18), followed by ‘Deep Subsoil’ (RPI = 1.19) and ‘Subsoil’ (RPI = 1.55). With respect to landcover classes, as well as according to R² and RMSE, the model performs best for ‘Cropland’ (R² = 0.71; RMSE = 40.45), followed by ‘Grassland’ (R² = 0.69; RMSE = 54.17) and ‘Woodland’ (R² = 0.60; RMSE = 58.83). The uncertainty increases from ‘Cropland’ (RPI = 1.16) to ‘Woodland’ (RPI = 1.24) and to ‘Grassland’ (RPI = 1.52;

Table 10. Performance of the CEC_pot model according to accuracy (R² and RMSE) and uncertainty (RPI) based on independent validation data. Provided values are averaged across five independent submodels and refer to the entire study area (‘Total area’) across all depth intervals, as well as to specific depth intervals (‘Topsoil’, ‘Subsoil’, and ‘Deep Subsoil’) and landcover classes (‘Cropland’, ‘Grassland’, and ‘Woodland’).

Indicator	Depth	Total Area	Cropland	Grassland	Woodland
R2	0–120 cm	0.72	0.71	0.69	0.60
	0–30 cm	0.71	0.71	0.68	0.48
	30–60 cm	0.67	0.76	0.62	0.60
	60–120 cm	0.59	0.66	0.57	0.31
RMSE	0–120 cm	48.61	40.45	54.17	58.83
	0–30 cm	50.65	42.56	52.84	65.65
	30–60 cm	44.61	34.27	59.14	38.87
	60–120 cm	45.08	40.2	52.16	45.71
RPI	0–120 cm	1.31	1.16	1.52	1.24
	0–30 cm	1.18	1.04	1.33	1.17
	30–60 cm	1.55	1.38	1.94	1.32
	60–120 cm	1.19	1.06	1.29	1.22

; Figure 9 and Figure S7).

The CEC_pot model performance is robust, referring to R², showing low variability across independent submodels (SD_R2 = 0.02). However, the SD_RMSE is 1.96, which is relatively large. Accuracy indications based on the IV are comparable to those derived from OOB_error (R² = 0.69; SD_R2 = 0.01; RMSE = 48.87; SD_RMSE = 0.55) and CV10 (R² = 0.69; SD_R2 = 0.01; RMSE = 49.71; SD_RMSE = 0.54).

The covariate importance analysis revealed a high influence of covariates related to terrain and spectral bare soil. Covariates related to climate and land use were less important for the CEC_pot model (Figure 9).

With respect to map values, the mean CEC_pot decreases from ‘Topsoil’ to ‘Subsoil’ and ‘Deep Subsoil’ between 3% and 34%. The mean CEC_pot is lower for ‘Cropland’ (137.76 mm_ckg⁻¹) compared to ‘Grassland’ (179.21 mm_ckg⁻¹) and ‘Woodland’ (184.79 mm_ckg⁻¹; Figure 9, Figure 10 and Figure S11).

4. Discussion

4.1. Evaluation Strategy

We applied stratified random data splitting to determine the prediction accuracy, indicated by averaging the R² and RMSE of multifold independent data splits and models, respectively. Data splitting, in combination with the aforementioned accuracy indicators, belongs to the prevalent accuracy estimation strategies for DSM [5,22]. The data splits were stratified, referring to depth intervals and landcover classes. This allows us to derive depth- and landcover-related accuracy estimations while using a single model calibrated across all strata. In comparison, recent DSM studies using a 3D model design as presented here disregarded stratification in the validation data and were limited to overall accuracy estimations across depth intervals or geographic regions [32,33,34].

The presented study provides a spatially explicit measure of the prediction uncertainty, framed by machine learning DSM for continuous soil property data. Generally, the uncertainty measure describes the limited ability of any DSM model to reproduce the soil–landscape relationship due to deficiencies in the soil sample and covariate data among others [15,141,142]. Uncertainty maps are particularly required for assessing broad-scale DSM products since the model is transferred to a variety of spatial domains, which may be inconsistently represented by the model [143,144,145]. Moreover, uncertainty maps, derived a priori based on archived soil data, can be leveraged as an information base for planning upcoming soil surveys. Accordingly, several studies have proposed sampling designs that are optimized by targeting areas of high uncertainty for improved model performances combined with reduced field costs [13,15,146].

Since we applied QRF, PIs and thus our uncertainty measure (RPI) harness the built-in mechanism of QRF to predict conditional quantiles [141]. While the QRF-based uncertainty estimation has been generally approved [10,147], the RPI measure (RPI) is further approved by our results, showing corresponding trends between RPI and the accuracy indicators (R², RMSE) across depth intervals and landcover classes.

Contrary to this, Kasraei et al. [25] proposed an uncertainty approach, applying a QRF-like quantile regression as a post-processing mechanism and thus being independent from the prediction model. Moreover, bootstrapping is another technique frequently applied in DSM that allows to derive uncertainty estimations independently from the type of prediction model. According to the bootstrapping method, multiple map realizations and thus PIs are derived based on randomized subsets of the training data [15,148]. Alternatively, Malone et al. [149] introduced an empirical method to derive spatially explicit PIs independently. The method clusters the covariate space based on a fuzzy K-means routine, then PIs are estimated based on the model errors within each cluster and finally regionalized according to the grade of membership to each cluster. Recently, Meyer and Pebesma [143] proposed a method to derive spatial uncertainties based on comparing the covariate spaces between a specific location and the calibration data rather than using PIs. This method acknowledges that prediction uncertainties increase when the covariate space of specific locations becomes increasingly dissimilar to the covariate space of the calibration data. The aforementioned approaches might be appealing since they are independent from the prediction model and applicable to a wide range of DSM designs. However, this also might include a more complex workflow, computational demand, and increased requirements related to data distributions and sample density, which are impractical in the context of broad-scale DSM. Instead, the use of QRF provides a straightforward workflow that integrates adequate model performance with reliable uncertainty estimations based on a single calibrated model [32].

4.2. Model Performance and Soil Sample Data

Numerous DSM studies identified the spatial sample density (SSD, n_xy*km⁻²) as a major driver for the model performance [150,151,152]. Accordingly, accuracy levels increase with an increasing SSD up to a certain threshold, which depends on the target soil property, model design, and environmental setting [64,153,154]. Chen et al. [50] reviewed 244 recent broad-scale DSM studies and reported an SSD variation between 0.0001 and 1, while only 5% of these studies outline an SSD larger than 0.1. Most broad-scale DSM studies in Europe showed an SSD between 0.01 and 0.1. A Czech DSM study by Žížala et al. [42] reported a relatively high SSD of 0.83 for Clay and 0.07 for SOC. A French DSM study by Mulder et al. [155] reported an SSD of 0.03. Zhou et al. [156] mapped SOC for Switzerland using a soil sample set with an SSD of 0.005 only. Generally, the SSDs of our study are considerably higher, ranging between 0.06 and 0.56 across the target soil properties.

Our results for CEC_pot revealed increasing SSDs combined with improved R² values from deeper to lower soil depths, as well as from ‘Woodland’ to ‘Grassland’ and to ‘Cropland’ (Table 5 and Table 10). Since the SSD of CEC_pot is relatively low, the consistent trends indicate SSD as one driving factor for the DSM model’s performance. However, referring to Clay, SOC, and pH, which show relatively high SSDs, the accuracy trends across depth intervals and landcover classes are ambiguous, most likely since SSDs are above a certain threshold to influence accuracy levels (Table 5, Table 7, Table 8 and Table 9).

Besides a sufficiently high SSD, adequate DSM model performance depends on a uniform spatial distribution of the sample data to optimally represent the true soil variation within the study area [31,151]. However, similar to more than 50% of the broad-scale DSM studies, the present study is based on legacy soil data, which are not acquired purposefully for countrywide DSM [5]. Accordingly, our sample data show partly uneven frequency distributions and are spatially clustered towards ‘Cropland’ (Table 5; Figure 1, Figure 2 and Figure S6). In consequence, similar to other broad-scale DSM studies, our model performances show spatial inconsistencies (Table 7, Table 8 and Table 9; Figures S8–S11). The deficiencies of our sample data, indicated by a depth-related decrease in SSD and spatial clustering, can be addressed by the purposeful and uncertainty-guided acquisition of additional soil samples [13,15,146].

4.3. Model Performance in Comparison

The review by Chen et al. [5] reported median R² values of broad-scale DSM studies, predominantly published between 2014 and 2021. In comparison, the R² values of the presented study are 14% to 33% higher across the target soil properties, referring to 3D DSM studies that disregard a depth-specific model evaluation (cf. 4.1). With respect to ‘Topsoil’ only, the R² values of the presented study are 2% to 95% higher (

Table 11. Comparison of R² values between the presented study and the median of broad-scale DSM studies, reviewed by Chen et al. [5]. ‘Total’ refers to 3D DSM models that are not evaluated for specific depth intervals.

Soil Property	Depth	R2—Chen et al. [5]	R2—Presented Study	Change
Clay	0–30 cm	0.44	0.69	+57%
	Total	0.52	0.69	+33%
SOC	0–30 cm	0.49	0.50	+2%
	Total	0.56	0.64	+14%
pH	0–30 cm	0.60	0.70	+17%
	Total	0.60	0.76	+27%
CECpot	0–30 cm	0.37	0.71	+95%
	Total	0.59	0.72	+22%

).

The accuracy levels of the presented study are widely comparable with the most recent broad-scale DSM studies. Lu et al. [19] applied an extremely randomized RF model for mapping ‘Topsoil’ pH at the European scale and a spatial resolution of 250 m × 250 m, reporting an R² of 0.70. Žížala et al. [42] mapped Clay, SOC, and pH at the Czech national scale with a spatial resolution of 20 m × 20 m using a model design similar to the presented study. In comparison to Žížala et al. [42], our study shows higher R² values for modeling Clay, with +13% for ‘Topsoil’ and +33% for ‘Subsoil’. Referring to SOC, the R² values of our study are 25% higher for ‘Topsoil’ and 20% lower for ‘Subsoil’. The R² for our pH model is 59% higher. Furthermore, Liu et al. [157] mapped ‘Topsoil’ pH, SOC, and CEC_pot at a spatial resolution of 90 m × 90 m for the entire territory of China based on QRF. Referring to the RMSE, the accuracy levels of our study are higher, with +21% for pH, +13% for SOC, and +37% for CEC_pot.

In the Swiss context, the results of the presented study outperformed existing national-scale DSM studies. Baltensweiler et al. [35] mapped pH, SOC, and Clay for the Swiss ‘Woodland’ at a spatial resolution of 25 m × 25 m using RF. In comparison, our R² values for ‘Woodland’ pH predictions are higher, with +31% for ‘Topsoil’ and +68% for ‘Subsoil’. The SOC predictions of the presented study are higher, with +32% for ‘Topsoil’ and +104% for ‘Subsoil’. The R² values of our Clay predictions are higher, with +35% for ‘Topsoil’ and by 75% for ‘Subsoil’. Zhou et al. [156] applied Boosted Regression Trees to map SOC for the entire Swiss territory in a spatial resolution of 100 m × 100 m based on 150 topsoil samples only (0–20 cm). Our study shows a 6% higher R² value referring to ‘Topsoil’ (0–30 cm). Furthermore, Gupta et al. [64] mapped SOC and Clay for the entire territory of Switzerland at a spatial resolution of 30 m × 30 m using QRF. While their accuracy levels were not analyzed for specific depth intervals, our study shows a higher overall R² with +21% for the SOC and +73% for the Clay model.

4.4. Influential Covariates for the Soil Property Distributions

Due to the wide availability of DEMs, combined with the integrative character of the terrain to represent multiple scorpan factors, terrain covariates are widely used for DSM and often drive the model performance [5,60,158]. Numerous DSM studies found a dominant role for terrain covariates in predicting various soil properties in different environments [35,51,159]. Particularly, spatially multiscale terrain covariates are approved to be beneficial for DSM [42,60,61,160,161]. This is confirmed by the results of the presented study, which designate the group of multiscale terrain covariates as highly influential for the model performance across all target soil properties (Figure 3, Figure 5, Figure 7 and Figure 9).

Climate influences soil formation due to interactions with hydrological cycles, whereby these interactions are related to changes in soil biology, soil physics, and soil hydrology [20,162,163]. Accordingly, the results of numerous DSM studies revealed a high influence of climate covariates in predicting Clay, pH, SOC, and CEC_pot [19,164,165,166]. In the presented study, climate covariates are a dominant factor for the Clay model but are medium-important for pH and less-important for predicting SOC and CEC_pot. However, since our model performances for pH, SOC, and CEC_pot are reasonable (Table 8, Table 9 and Table 10), the underlying soil-forming processes seem to be adequately described by the interactions of other covariates than those directly related to climate (Figure 3, Figure 5, Figure 7 and Figure 9).

Land use and its change influence soil carbon and nutrient cycles, thereby changing soil biology and soil physics [77,167,168]. Several DSM studies reported the relatively high importance of land use covariates for predicting Clay, pH, SOC, and CEC_pot [19,80,81,169]. In the presented study, the covariate importance analysis indicates a lower influence of land use covariates for predicting Clay, SOC, and CEC_pot, while only the pH model is widely driven by covariates related to land use. However, land use patterns are directly influenced by specific climatic and terrain conditions, and thus, land use covariates are likely superimposed by those covariates (Figure 3, Figure 5, Figure 7 and Figure 9) [20].

Covariates related to the spectral signature of bare soil correlate with soil properties such as Clay, pH, SOC, and CEC_pot in varying environments [38,40,56,73,170,171,172]. Thus, as suggested by Demattê et al. [67], spectral bare soil covariates are increasingly used as an integrative part of DSM covariate sets. For example, the Czech DSM study by Žížala et al. [42] reported a high dominance of spectral bare soil covariates for predicting SOC and a medium impact targeting Clay and pH. Also, a French DSM study by Urbina-Salazar et al. [75] identified covariates related to the spectral reflectance of bare soil as a dominant factor for mapping SOC. These findings are widely confirmed by the presented study, where spectral bare soil covariates are identified as a dominant factor for predicting SOC and CEC_pot, as well as medium influential for predicting Clay and pH (Figure 3, Figure 5, Figure 7 and Figure 9). Since the use of spectral bare soil covariates for DSM emerged quite recently, the diversity of available techniques for filtering and compositing bare soil pixels is high [38,68,70,73,138,173]. In this context, Heiden et al. [174] emphasized that the selection of adequate spectral indices to filter bare soil pixels is of particular importance to challenge disturbing factors related to, e.g., varying soil moisture conditions and soil structures, as well as crop residue cover [67]. A wide range of spectral indices, including NDVI, NBR2, Bare Soil Index (BSI), or adaptations thereof, have been applied [38,40,138,174,175]. In the present study, we used a combination of NDVI and NBR2 as two of the most commonly used indices for bare soil detection with respect to LA data and various environments [38,40,69,70,176,177]. However, although our spectral bare soil covariates outlined a major importance for the soil model performances, a more detailed analysis of the spectral indices for bare soil detection is indicated with respect to future applications [174,175]. Furthermore, several studies reported promising results from leveraging hyperspectral imagery for bare soil detection and associated soil mapping [177,178,179,180]. Accordingly, the current advent of available hyperspectral data based on the satellite platforms PRISMA and EnMAP has high potential for future DSM applications at broad scales [177,181,182].

5. Conclusions

This study provides a DSM framework to map multiple soil properties in three dimensions for the entire territory of Switzerland based on machine learning as well as soil and remote sensing data archives. The DSM approach targets the soil properties of clay content, organic carbon concentrations, pH, and the potential cation exchange capacity for three depth intervals. A Quantile Regression Forest model was applied to link archived soil data with a comprehensive set of covariates derived from remote sensing data archives, a digital elevation model, and climate data, describing scorpan factors multiscale in space and time. Accordingly, the covariate set describes the spatially multiscale character of terrain, the spectral signature of bare soil, temporarily dynamic and multiscale patterns of land use, as well as climate and its associated temporal trends. The applied model design is particularly suitable for complex broad-scale DSM, allowing to efficiently derive predictions and uncertainties for multiple depth intervals based on a single calibrated model. Further, the model design allows us to evaluate the model’s accuracy with respect to crop-, grass-, and woodland areas, as well as refer to specific depth intervals. The results show robust model performances, which are superior compared to previous DSM approaches at the national scale in Switzerland and in relation to the majority of recent broad-scale DSM approaches in other countries. Moreover, the results indicate that a spatially uniform and high soil sample density is a decisive factor for adequate DSM performance. The results suggest that in Switzerland, additional soil samples are particularly required for woodland and grassland areas, as well as for the subsoils. The multiscale remote sensing covariates were approved to be beneficial for DSM. Particularly, covariates related to the spectral signature of bare soil were found to be a major driving factor for DSM targeting organic carbon and cation exchange capacity. Generally, this study reveals the potential to update existing broad-scale DSM approaches by combining machine learning with a comprehensive exploitation of soil and remote sensing data archives.