Comparison of catchment scale 3D and 2.5D modelling of soil organic carbon stocks in Jiangxi Province, PR China

As limited resources, soils are the largest terrestrial sinks of organic carbon. In this respect, 3D modelling of soil organic carbon (SOC) offers substantial improvements in the understanding and assessment of the spatial distribution of SOC stocks. Previous three-dimensional SOC modelling approaches usually averaged each depth increment for multi-layer two-dimensional predictions. Therefore, these models are limited in their vertical resolution and thus in the interpretability of the soil as a volume as well as in the accuracy of the SOC stock predictions. So far, only few approaches used spatially modelled depth functions for SOC predictions. This study implemented and evaluated an approach that compared polynomial, logarithmic and exponential depth functions using non-linear machine learning techniques, i.e. multivariate adaptive regression splines, random forests and support vector machines to quantify SOC stocks spatially and depth-related in the context of biodiversity and ecosystem functioning research. The legacy datasets used for modelling include profile data for SOC and bulk density (BD), sampled at five depth increments (0-5, 5-10, 10-20, 20-30, 30-50 cm). The samples were taken in an experimental forest in the Chinese subtropics as part of the biodiversity and ecosystem functioning (BEF) China experiment. Here we compared the depth functions by means of the results of the different machine learning approaches obtained based on multi-layer 2D models as well as 3D models. The main findings were (i) that 3rd degree polynomials provided the best results for SOC and BD (R2 = 0.99 and R2 = 0.98; RMSE = 0.36% and 0.07 g cm-3). However, they did not adequately describe the general asymptotic trend of SOC and BD. In this respect the exponential (SOC: R2 = 0.94; RMSE = 0.56%) and logarithmic (BD: R2 = 84; RMSE = 0.21 g cm-3) functions provided more reliable estimates. (ii) random forests with the exponential function for SOC correlated better with the corresponding 2.5D predictions (R2: 0.96 to 0.75), compared to the 3rd degree polynomials (R2: 0.89 to 0.15) which support vector machines fitted best. We recommend not to use polynomial functions with sparsely sampled profiles, as they have many turning points and tend to overfit the data on a given profile. This may limit the spatial prediction capacities. Instead, less adaptive functions with a higher degree of generalisation such as exponential and logarithmic functions should be used to spatially map sparse vertical soil profile datasets. We conclude that spatial prediction of SOC using exponential depth functions, in conjunction with random forests is well suited for 3D SOC stock modelling, and provides much finer vertical resolutions compared to 2.5D approaches.

Introduction Soils are a fundamental part of ecosystem functioning and services [1]. As finite resources, soils contribute to food production, nutrient cycling, biodiversity and freshwater quality [2]. Furthermore, they are interconnected with other ecosystem functions and services, such as local and global climate alteration; and therefore, contribute indirectly to human well-being [3]. Among soil properties, soil organic carbon (SOC) plays an important role in this context. SOC increases the water-holding capacity (e.g. important for agriculture, forest and flood management), improves the physical properties of soils, such as nutrient availability for plants in agriculture and forestry, and accounts for carbon sequestration to mitigate climate change [4][5][6]. In forestry, there is strong interest in the effects of tree species and tree diversity on soil carbon input and mineralization as well as the net effects of these processes [7]. Knowledge about the interconnection between SOC, forests and the diversity of tree species as well as SOC stock degradation by soil erosion [8,9] and land cover change [10,11] can also help to implement countermeasures to reduce global warming [7]. Consequently, the implementation of a credible soil carbon auditing and monitoring to verify changes in SOC is crucial regarding soil security and carbon sequestration [7,12].
To preserve the functions and services provided by soils, a good quantitative understanding of the SOC stocks is required-both in the vertical domain of a soil profile as well as in the spatial domain over landscapes [13,14]. However, conventional soil maps use soil classes in horizontal dimension and soil horizons in vertical dimension. This categorical setup is often not precise enough and not well suited for interpreting soil functions and processes as well as for decision-making, since soil properties mostly vary continuous in space and time [15,16].
For the spatial prediction of continuous soil properties, such as SOC, methods of digital soil mapping (DSM) are suitable [17][18][19]. DSM is based on the soil forming factor concept [20] and the scorpan model introduced by McBratney et al. [21]. Both approaches illustrate soil information as a function of environmental covariates, influencing the process of soil formation. Terrain parameters, describing the shape of the land surface, are used widely as an environmental covariate in DSM. Terrain is an essential factor of soil formation and controls the effects of gravity, climate, lithology, water and biota [22][23][24]. Hence, models that are based on For spatially modelling soil properties, different approaches have been established to derive relationships between soil properties and environmental covariates. However, for a reliable estimation of SOC stocks, the vertical dimension is crucial [13]. A common way of threedimensional mapping is to consider the vertical dimension as multiple two-dimensional predictions, which can be interpreted in a three-dimensional way [17,[27][28][29]. Because, multi-layered predictions do not provide full 3D soil information, since they are limited to the mapped depth increments. Information of the space between the mapped depth increments has to be derived on an interpretative and subjective basis. One approach is to vertically interpolate the single layers to construct a volumetric model, which is computationally intensive [30,31].
Therefore, multi-layered models are referred to as pseudo-3D mapping or 2.5D mapping [32]. To overcome these drawbacks, it is favourable to map soil properties as continuous depth function in the spatial domain [13,18], where the vertical distribution of soil properties is represented by depth functions, that are predicted spatially. These predictions allow the calculation of SOC stocks over the integral of the functions [33] as well as the calculation of fully three-dimensional maps at any vertical resolution [32,[34][35][36][37].
While with 2.5D mapping soil properties are directly predicted at specific depth levels using the environmental covariates [17,29], 3D approaches use environmental covariates to predict parameters of the depth functions [34], which are abstract soil properties. According to the scorpan model, soil properties can be spatially mapped with neighbourhood relations solely [21], which also have been used for 3D modelling [36,40,42,43]. Over the past years, machine learning techniques have become a standard technique in DSM due to several advantages like dealing with non-linearity or the handling of large datasets. Aldana Jague et al. [33] used multiple linear regression (MLR) to model SOC incorporating terrain covariates, while Gasch et al. [43] compared spatial and terrain covariates using random forests (RF) and regression kriging for mapping SOC at different depth layers. Piikki et al. [27] used multivariate adaptive regression splines (MARS) to model clay and sand fractions as well as organic matter based on proximal soil sensing data. Several other studies also suggest that machine learning techniques, such as artificial neural networks (ANN; [41,44]), random forests (RF; [17]) and support vector machines (SVM; [45]), can be applied successfully in DSM.
The objectives of this study were to test the spatial prediction of four soil profile depth functions for modelling SOC content and bulk density with different machine learning methods based on multi-scale terrain covariates. The tested soil profile depth functions are polynomials of 2 nd and 3 rd degree, natural logarithmic and exponential functions. The machine learning methods used to model the depth functions spatially were multivariate adaptive regression splines (MARS), random forests (RF) and support vector machines (SVM) with radial basis functions. We validated the machine learning models with 10-fold cross-validation and evaluated the results of the 3D mapping approach by comparing it with the predictions of the more common multi-layered 2.5D modelling approach based on five layers.

Study area and sampling design
The biodiversity and ecosystem functioning (BEF) China project [46] is located near Xingangshan, Jiangxi Province, PR China (UTM/WGS84: 50R 588000 3222000), about 400 km south-west of Shanghai (Fig 1). The study site is a topographically heterogeneous environment in a small catchment of 26.7 ha leased by the Institute of Botany of the Chinese Academy of Sciences (CAS). It features an elevation ranging from 105 to 275 m a.s.l., slopes inclined 29˚in average and a maximum slope inclination of 45˚, which are typically convex [19]. Non-calcareous slates with varying sand and silt content and grey-green sandstone constitute the bedrock. Predominant soil types are Endoleptic Cambisols with Anthrosols at the hillsides and Gleysols at the valley bottom. The mean soil depth is 0.6 m with underlaying isomorphic weathered slate (saprolite; [19]). Soil texture ranges from silt loam to silty clay loam [47]. The climate is typically subtropical with monsoons in summer, a mean annual temperature of about 17˚C and long-term average annual rainfall of about 1800 mm [48] but with a drier period from 2009 to 2012 [49].
About 18 ha were covered with 271 experimental plots. In total 8.7 ha at the valley bottom were not part of the experimental design due to paths and rivulets. Plots had a size of 25.8 m × 25.8 m (traditional Chinese unit of 1 mu, 1/15 ha) and were replanted in 2008 after clearcut of a commercial Chinese fir plantation. One plot comprised 400 (20 × 20) trees in monocultures and mixtures of 2, 4, 8, 16 and 24 species. Species composition of the plots was based on random as well as non-random (plant trait-oriented) extinction scenarios, where all species were represented equally (broken-stick design). The datasets used in this study comprised soil samples from random subsets of all species and species richness levels referred to as VIPs (Very Intensively Studied Plots). For details on the experimental design, see Bruelheide et al. [46] and Trogisch et al. [50].

Datasets
All described datasets are part of the legacy database of BEF-China. Soil sampling was conducted in 2014. Nine cores on a regular grid basis (3 cm in diameter) were taken at each of the 67 VIPs according to the BEF-China experimental design (Fig 1; [46]). The samples were bulked for each depth increment (0-5 cm, 5-10 cm, 10-20 cm, 20-30 cm and 30-50 cm) and were referred to as dataset SOC (n = 67; Fig 2). Fine roots and charcoal were sorted out manually. For dry combustion CNS-analysis, a Vario EL III (Elementar, Hanau, Germany) was used. Due to acidic soil conditions there was no detectable carbonate fraction, and thus total carbon represented SOC [19]. SOC content ranged from 5.06 to 0.35% decreasing with depth.
Bulk density samples (n = 55) were taken in April 2015 with soil sample rings (100 cm 3 ) and five replicates for each depth increment at the VIPs. Bulk density was determined gravimetrically and was referred to as dataset BD (Fig 2). Bulk density ranged from 0.75 to 1.84 g cm -3 increasing with depth.
Since some plots with SOC samples did not have BD data (Fig 1), both soil properties were modelled individually instead of calculating and modelling the SOC stocks directly. This 'model-then-calculate' approach is a useful alternative to the 'calculate-then-model' approach. Both were compared by Orton et al. [51]. The digital elevation model (DEM) had a resolution of 5 m and was generated by ordinary kriging [52] based on differential global positioning system data (DGPS) with 1956 points (73 points per ha; [19]). The distribution of datasets SOC and BD over the DEM is shown in Fig 3. Dataset SOC covered the elevation data more comprehensively compared to the dataset BD.

Digital terrain analysis
Environmental covariates that describe the morphometry of a landscape are grouped in four major classes of terrain attributes: local, regional, combined (i.e. combinations of local and regional) and solar morphometric variables. Given that many terrain attributes can be calculated based on different equations or modelling approaches and because it is unknown which version would be most suitable for modelling SOC and BD within the study area, we used multiple established methods to derive single terrain attributes, if available. Given the circular nature of aspect, we used sine and cosine transformations to derive eastness and northness. Overall, we calculated 58 terrain attributes (Table 1) with SAGA GIS 2.3.1 [53].
Terrain attributes derived from a DEM with a given resolution may not be suitable for landscape characterization and for digital soil mapping due to a non-representative DEM resolution [63], since the terrain attributes are not derived on the most relevant scale [64,65]. To examine the influence of scale, [65] applied simple smoothing (mean) filters with different neighbourhood sizes. This approach was applied on every terrain attribute used in this study with five circular neighbourhoods (radii of 1, 2, 4, 6 and 8 pixels), resulting in 290 terrain attributes in total. The maximum radius was set to 8 pixels to represent the local catena scale of 90 m.

Machine learning techniques
We compared three data mining methods to test the 3D prediction of soil profile depth functions for SOC and BD based on terrain covariates. Given the large number of 290 covariates (instances) and sample sizes of n = 67 and n = 55, not all available techniques could be applied. For example, the interpretable multiple linear regression (MLR) analysis used for spatial modelling of polynomial depth functions by Aldana Jague [34] requires more samples (n) than instances (p; [66]). Furthermore, we have to account for multi-collinearity. Many terrain covariates in this study are calculated by different algorithms for the same terrain attribute and on different spatial scales with the same algorithm, which is often seen as a constraint in machine learning [66]. To reduce the covariate space to either enable MLR or handle the 'curse of dimensionality', principal component analysis (PCA) is often applied. However, feature reduction with PCA can have negative effects on model accuracy with multi-scale terrain data and models with the full set of covariates have higher accuracies [65]. Other feature reduction methods increase accuracy only marginally [65]. In this study, we applied multivariate adaptive regression splines (MARS), random forests (RF) and support vector machine (SVM). These machine learning methods are robust against multi-collinearity, can handle  Table 1. Terrain attributes used for SOC and bulk density modelling.

Covariates
Method Author(s)
Multivariate adaptive regression splines (MARS). MARS was introduced by Friedman [69] and is a generalisation of recursive partitioning regression approaches using piecewise linear models. With its linear basis functions, it overcomes the discontinuous response of other recursive partitioning models like Classification and Regression Trees (CART; [70]) and can generate continuous surfaces. Therefore, prediction accuracy of MARS is expected to be higher [69]. MARS is a partial linear function, where each new part is added with an exhaustive search for best fit and models a finite quantity of the regression. Thus, the model measures variable importance by its nature and is insensitive to non-informative instances. MARS require very little pre-processing and are non-affected by collinearity, since the predictor selection is random during iteration and redundant features are used equally [66]. This may affect measurement of variable importance and interpretation, which, however, is out of scope in this study. For modelling using MARS, the earth package version 4.4.6 [71] was used.
Random forests (RF). RF is a widely used machine learning technique in digital soil mapping [17,22,64,72]. It was introduced by Breiman [73] and is an ensemble technique with CART [70] as a base learner. The single decision tree uses binary splits to create more homogenous groups in respect to the response. To grow an ensemble of trees, different random subsets of covariates (bootstrap sampling) and features (random set of features for every split) are used to build a single tree. The final prediction is created by averaging all individual tree outputs. Breiman [73] has proven that random forests with a large number of trees is robust against overfitting. Moreover, it is robust against noise, non-informative and correlated features. RF also returns feature importance measures (affected by correlation as MARS; [66]) and there is little need for fine-tuning [74]. The randomForest package version 4.6-12 [75] was used for modelling with RF. Support vector machine (SVM). Originally, SVM has been developed for classification problems [76]. It is a kernel method and uses hyperplanes to linearly separate classes of objects. For regression problems, Drucker et al. [77] developed support vector regression machines (SVR), which are an extension of SVM. Therefore, the term SVM is often used in both cases. The kernel function defines a transformation of the input data into a high dimensional feature space. In this feature space, it is possible to derive a linear regression hyperplane for non-linear relationships. Afterwards, it is back-transformed to non-linear space. Smola and Schölköpf [78] provide a comprehensive and detailed insight into SVR. The kernel used in this study is a radial basis function, where the scaling parameter σ is estimated by caret after a method by Caputo et al. [79]. In contrast to MARS, Drucker et al. [77] suggest that SVM should be used when the number of features is larger than the number of instances, since its optimisation does not depend on the dimensionality of feature space. Furthermore, SVM is partially insensitive to outliers (depending on cost factor) and does not require feature reduction to reduce multi-collinearity [66]. The kernlab package version 0.9-25 [80] was used for radial support vector regression modelling.
Data pre-processing. Some algorithms are sensitive to the scale and the range of the covariate space (e.g. SVM). To reduce effects of small values and little variance, SVM needs centred and scaled covariates [66], which was computed using the scale and centre-option in caret. To make all models comparable, this was also done for MARS and RF.

Spatial 2.5D and 3D models
Differences between 2.5D models and spatial prediction of depth functions. The environmental covariates were used to train regression models (MARS, RF and SVM) to predict SOC and BD. For 2.5D predictions this was done for each sampled depth increment individually, were we assigned the mid-depth of the sampled increments as depth of the respective layer. This method to obtain volumetric soil information has several advantages. For modelling of each standard depth individually, there are no further requirements to abstract soil information in terms of vertical variability, i.e. a soil profile function. Furthermore, there is no error propagation through secondary models that describe depth functions. On the other hand, in contrast to 3D modelling, 2.5D modelling has the disadvantage that the individual model outcomes are purely two-dimensional. Soil properties of the depth increments between the standard depths are not used in the model and have to be derived on an interpretative basis [15] or through further processing [30] after spatial prediction. However, this is a well-established and well-documented approach. Therefore, we compare the results of the 3D approach described below directly with the 2.5D results.
3D mapping with soil depth functions. For the spatial modelling of depth functions, which we handled similar to the soil properties in terms of modelling, we applied 3 rd degree polynomial functions proposed by Aldana Jague [34] and less flexible 2 nd degree polynomials as well as logarithmic and exponential functions [32]. The workflow of the 3D mapping (Fig 4) of this study involved five main steps: i). Mathematical approximation of depth functions to the five depth increments with a linear least squares approach. These were [34] [32] cf. [32] where f 1,2,3,4 (x) is SOC and BD at a specific depth x (depth of the lower corner of a voxel in cm), c 0 is the intercept that equals SOC and BD at depth 0 (cm) and the function coefficients c 1 , c 2 and c 3 are dimensionless. This altogether described the vertical distribution of SOC in respect to depth x at a certain location. ii). Evaluation of model error for all equations in (i).
iii). Spatial modelling of the function coefficients c 1 , c 2 , c 3 and c 0 (analogous to two-dimensional modelling of SOC and BD) of the depth function with the lowest error (ii) with MARS, RF and SVM. The depth function parameters were treated and evaluated similar to a soil property.
iv). Evaluation of the cross-validation results for MARS, RF and SVM models of the depth function coefficients. v). Solving the depth functions with spatially modelled coefficients (iii) at each grid location to generate a three-dimensional model. The depth functions were solved for depths from 0 cm to 50 cm in 5 cm increments. The resulting 11 depth layers (matrices) were stacked to two three-dimensional models (one for Catchment scale 3D and 2.5D modelling of soil organic carbon stocks SOC and BD each), where individual values are represented by voxels, which are the volumetric 3D analogue of 2D pixels. Due to the nature of the polynomial depth functions, negative SOC predictions in the profiles are possible. Consequently, the values of these voxels had to be set to zero. This is not required for logarithmic and exponential functions.
Compared to the standard depth method, the main advantages of spatially modelled depth functions are a higher vertical resolution and the fact that the result can be interpreted as volumetric structure. Instead of pixels with SOC and BD information in multiple layers, volumetric elements-so called voxels-in a three-dimensionally georeferenced stack of matrices with userdefined vertical resolution are obtained. Since the depth functions are secondary models, the error which is propagated by the depth function model to the spatial model depends on the chosen function. Due to the limited number of samples per profile, cross-validation of the depth functions was omitted.
The final models for SOC and BD were validated internally against the measured values of the input datasets.

Validation and evaluation
The evaluation consists of two independent steps for the 2.5D multi-layered model predictions and the volumetric 3D model predictions of SOC and BD, where we treat the depth function parameters as soil properties.
In a first step, we evaluated each model of the soil properties SOC and BD as well as the spatial models of the depth function parameters, by using a 10-fold cross-validation with the coefficient of determination (R 2 ) and the root mean square error (RMSE) as quality criterion. In this step, the models were tuned over the default grid-or hyper-learning sequence of parameters [81] using the tune grid function of caret to identify the most suitable combination of tuning parameters with the lowest RMSE and to reduce the model error, while preserving the models ability to generalise. The tuning parameters are degree and nprune for MARS, mtry for RF and cost for SVM. For RF ntree was set to the default value and σ for SVM was calculated by a method after Caputo et al. [79]. All models used the same set of folds to make cross-validation results comparable. The final models were selected from this sequence by the lowest RMSE.
To estimate the effect of overfitting of the depth function models based on grid learning, we evaluated the 3D model results with the datasets SOC and BD by R 2 and RMSE (observed-predicted-evaluation). Overfitting is indicated by large differences in the prediction error between the training and the validation sets [81].
Further, we compared the 3D models against the 2.5D predictions of the same datasets to evaluate the performance of the 3D models. We chose this approach, because the legacy datasets are too small to hold out a larger subset for independent validation. The model results should be similar, if the spatial prediction of depth function parameters is reproducing the spatial distribution of the soil properties. This means that independently from the modelling framework (modelling of SOC and BD or modelling depth function as soil property) the results of the 3D model are reasonable, if both models are similar.
We see this comparison as a valid method for the evaluation of the 3D models, since Brus et al. [38] report strong correspondence between 2.5D and 3D geostatistical models and MARS, RF and SVM are well established for 2D and 2.5D soil mapping and in data science [17,27,66]. Therefore, we use the 2.5D layered predictions at the specific mid-depth of the increments as reference predictions. For the comparison between the 2.5D models and the corresponding depths of the 3D models, we used the coefficient of determination R 2 , Lin's concordance correlation coefficient (ρ c ; [82]), which validates the models against the 1:1 line, and the RMSE.

Estimation of SOC stocks
The three-dimensional array of SOC stocks was calculated by where SOC stocks (g voxel -1 ) is the soil organic carbon storage, SOC is SOC content (%), BD is bulk density (g cm -3 ), 500 2 is the base area of a voxel (cm 2 ) related to the DEM resolution of 500 cm and 5 is the vertical resolution in cm. Consequently, 1 voxel represented 1.25 m 3 of soil. Adjustment with the fraction of coarse material (> 2 mm) was omitted, since the coarse fraction was negligible low (< 5 vol.-%) at the VIPs and cannot be determined precisely by coring. According to Orton et al. [51] calculating the SOC stocks from two models of SOC and BD is an useful alternative when the samples are not taken at the same locations.

2.5D predictions of standard depths as reference
For the models of SOC, the mean cross-validation R 2 of MARS was 0.33 with a root mean square error of 0.39%, compared to RF with an R 2 of 0.41 (RMSE 0.34%) and SVM with an R 2 of 0.39 (RMSE 0.35%; cf. Table 2). Models for BD showed a mean R 2 of 0.43 (MARS), 0.39 (RF) and 0.39 (SVM) and mean RMSE values of 0.09 g cm -3 (MARS), 0.08 g cm -3 (RF) and 0.08 g cm -3 (SVM). In addition to the mean values, Table 2 shows the prediction accuracies and the RMSE's for each depth increment and all three machine learning techniques of both SOC and BD.

Soil depth functions
For SOC, all equations showed R 2 values higher than 0.9 (0.99 for f 1 , 0.96 for f 2 , 0.96 for f 3 and 0.94 for f 4 ) with a RMSE ranging from 0.36 (f 1 ) to 0.7% (f 2 ). For BD, the performance in terms of R 2 was similar (RMSE = 0.07 g cm -3 ), except for f 3 with R 2 = 0.84 (RMSE = 0.22 g cm -3 ), which is the natural logarithmic function. The 3 rd degree polynomial (f 1 ) resulted in the best fits for SOC and BD. However, the general trend of SOC in the profiles was exponential ( Fig  2). Hence, both the 3 rd degree polynomial and the exponential functions were chosen for further spatial modelling and comparison in this study. With higher errors and without being able to reproduce the general trend in the profiles profile the 2 nd order polynomial (f 2 ) was omitted in the following steps.

Spatial modelling of soil depth functions
The cross-validation results for the machine learning methods applied on the depth functions (c.f. Table 3) showed, that the polynomial depth functions for MARS, RF and SVM for SOC were comparable in their goodness of fit with marginal differences (mean R 2 from 0.3 to 0.32). R 2 of the exponential depth functions ranged from 0.3 for MARS to 0.44 for RF. The models of the function coefficients could not be compared directly because c 0 represented the SOC in % and BD in g cm -3 , whereas c 1 , c 2 and c 3 were dimensionless. Hence, we compared these models by the normalised RMSE (nRMSE), which is the RMSE divided by the coefficients range (Table 3). The nRMSE showed little variation of around 0.18 for all coefficient predictions of the 3 rd polynomial depth function of SOC. RF had the lowest mean of nRMSE over all coefficients (0.17). The lowest nRMSE (0.09) for SOC was achieved by the exponential depth functions (RF and SVM).
The models based on the 3 rd degree polynomial depth functions of BD had a mean R 2 of about 0.23-0.4, while the mean nRMSE was about 3×10 4 , due to the low performance of models with c 3 . Given such high errors, none of the models could reasonably predict the 3 rd degree polynomial depth function for bulk density. The exponential function was not able to reproduce the vertical trend of BD. Thus, we used the logarithmic depth function, although it fitted the five depth increments least. However, these spatial depth function models performed better (mean R 2 from 0.36 to 0.45; nRMSE of about 0.16 for SVM).

Evaluation of 3D predictions
For the comparison of 3D models against the 2.5D reference predictions, we used the coefficient of determination R 2 , Lin's concordance correlation coefficient ρ c and the RMSE in corresponding depths (Table 4). The three-dimensional MARS prediction for SOC with the 3 rd degree polynomial depth function showed the largest difference to its counterpart. The prediction at 2.5 cm ranged from close to zero to 15% SOC compared to 1.5 to 4% SOC in the two-dimensional prediction (Fig 5). The other depth increments showed a similar pattern with values down to -15% SOC. For the 2.5 cm increment the performance of RF was slightly better than that of SVM, but subsequently dropped with increasing depth. Especially at 40 cm, but also at 25 cm and 15 cm, the three-dimensional prediction of RF differed more from the two-dimensional predictions than the three-dimensional predictions of SVM differed from their counterparts. There was no distinct over-or underestimation of RF, but random scattering between -4 and 4% SOC for 40 cm (Fig 5). SVM showed lower deviation at 15 cm, 25 cm and even 40 cm. There were less predictions with negative values and less scattering. The predicted depth intersections of spatially modelled depth functions corresponded to the two-dimensional predictions by SVM largely by R 2 and ρ c , while RMSE is low (Table 4).

Table 3. Performance of a 10-fold cross-validation for MARS, RF and SVM applied on function coefficients of a 3 rd degree polynomial (f 1 for SOC and BD with four coefficients) and natural logarithmic function (f 3 for BD with two coefficients
In contrast, the 3D predictions of RF and SVM based on the exponential function showed good correspondence for all five depth increments ( Table 4). The 3D predictions overestimated SOC for the 0-5 and 5-10 cm increments and underestimated it for 20-30 and 30-50 cm slightly due to the exponential nature of the equation, but there was no wide scattering as it was the case with the polynomial prediction for RF.
The results of the internal validation showed high correspondence between the chosen models (RF with exponential function for SOC and RF with logarithmic function for BD) and respective input data at all five sampled depth increments ( Table 5). The R 2 and RMSE values of the internal validation were similar to the validation results of the model comparison, indicating that model overfitting of both models is similar (Table 4). This partly accounts to the propagation error of the profile depth function. The spatial prediction of the exponential function for SOC had an average R 2 of 0.79 with an average RMSE of 0.33% and Table 4   Catchment scale 3D and 2.5D modelling of soil organic carbon stocks the prediction of the logarithmic function used for BD had a R 2 of 0.77 with an average RMSE of 0.14 g cm -3 .

SOC stocks
The 2.5D models showed SOC stocks of 61.9 Mg ha -1 from 0 to 40 cm, with 19, 14.7, 12, 8.9 and 7.3 Mg ha -1 in the individual depth increments (from surface downwards). The 3D model predicted 78.3 Mg ha -1 over the whole interval. The upper 20 cm of soil contained about 46.4 Mg ha -1 . This depth is often designated as topsoil [83,84] and is also a critical soil depth for modelling plant productivity and community assembly [85]. 31.9 Mg ha -1 SOC are stored in the subsoil from 20 to 40 cm. Considering that the rooting depth varies, depending on the species and individual age, a static discrimination between topsoil and subsoil may be not appropriate. The model showed that plants with shallow roots down to 5 cm mainly interacted with a carbon pool of 10.9 Mg ha -1 , whereas plants with roots in 25 cm depth interacted with a pool of 54.5 Mg ha -1 . Fig 6 shows

2.5D predictions of standard depths as reference
As RF returned the lowest error for the 2.5D models, this was the best choice for modelling SOC. SVM ranked slightly below. Compared to the results presented by Lacoste et al. [30], who used Cubist for 2.5D SOC stock mapping, the accuracy of our results was similar and reasonable.
However, the sampled VIPs do not represent the terrain of the study site adequately, since they were chosen based on species richness levels, which were distributed randomly, and not representative for the study site. For example, a representative sampling design could be achieved with Conditioned Latin Hypercube Sampling (cLHS) [72,86].
For bulk density SVM and RF performed equal by means of R 2 and RMSE and showed a similar pattern, especially at 15 cm and 40 cm. MARS performed least for BD. In general, RF resulted in the most stable predictions and is therefore recommended over SVM.

Evaluation of 3D predictions
The negative values in the prediction results and the pronounced difference between the 3D models, with predictions up to 15% SOC, and the 2.5D models indicated that MARS is not capable of adequately predicting the depth functions in space, although the cross-validation showed similar results as for RF and SVM models. The latter showed better correspondence between the 3D and the 2.5D models (Fig 5, Table 4). According to the results of the direct comparison between the multi-layered prediction and the corresponding depths in the 3D model, RF with exponential functions was most suitable for SOC modelling. RF and SVM with polynomials performed well at upper depth increments and less in lower increments. MARS models were not suitable of reproducing the 2.5D predictions. Lower performance of all techniques with polynomials in the lower depth increments may be referred to lower influence of the terrain as a driving factor to explain SOC accumulation and redistribution (e.g. by erosion). Other factors accounting for SOC redistribution in deeper soil horizons may be bioturbation or vertical transport in the liquid soil phase. Additionally, it is possible, that accumulation layers in the solum, that would reflect the lateral distribution, were not fully covered by the legacy dataset and, therefore, the interpretation remains difficult. All these processes and others relevant for SOC concentration as well as SOC stocks cannot be fully covered by a distinct set of terrain parameter and lead to a dilution effect by predicting the deeper horizons. Lower accordance of the models also may be referred to uncertain models of function coefficients (c 3 ) and (c 2 ), which have significant influence at greater depths (cubic and squared) and exponentiate up this error. Based on the results, we chose RF with exponential depth functions for three-dimensional mapping of SOC and the logarithmic depth function for BD.

SOC stocks
Compared to other studies in this area, the estimated SOC stocks were well in line. Scholten et al. [19] calculated mean SOC stocks of 70 Mg ha -1 for the upper 50 cm with the same data but a different approach. Chen et al. [83] compared five plantations with different species in five age groups and calculated SOC stocks for the upper 20 cm. Especially the age of the trees and shrubs and their biomass have a strong impact on SOC stocks. Very young forest communities showed SOC stocks ranging from 20 to 25 Mg ha -1 and plantations with older trees of 7 to 10 years 30-40 Mg ha -1 . The latter were slightly older than the trees of BEF-China, where 42 Mg ha -1 were predicted. Diverse species pools in these studies may explain differences. Tang et al. [87] found SOC stocks in the top 60 cm in bamboo forests ranging from 60 to 200 Mg ha -1 .
The introduced approach is capable of summing SOC stocks at any depth interval. Since topsoil depth varies spatially, conventional static assumptions of topsoil thickness can result in inaccurate SOC stock calculations for individual horizons. Incorporating spatial models of topsoil depth into 3D SOC stock mapping can overcome this drawback and help to improve ecological and biodiversity models as conducted in the BEF-China experiment. In particular, consideration of biotic predictors like forest biomass, tree species richness and functional plant diversity might further improve model fit and accuracy of estimated SOC stocks [14]. This would allow one to quantify terrain-specific effects of changes in forest cover and composition on SOC stocks. The developed models could also help to identify areas that are especially prone to loss of SOC stocks (e.g. by soil erosion or land cover change).
Furthermore, continuous three-dimensional SOC mapping can support models of a national SOC inventory. Yang et al. [37] applied depth functions to categorical soil types and estimated SOC stocks for mainland China. Combining models with high vertical resolution by Yang et al. [37] and continuous spatial modelling like in this study can improve accuracy of SOC mapping compared to the categorical mapping approach. This combination can also help to estimate and understand carbon fluxes between topsoil and subsoil [88] as well as between soil and the atmosphere [5]. Both objectives play major roles in inventory estimation, SOC auditing and decision making in respect to ecosystem services and carbon sequestration [1,5,7,12,89].

Conclusion
This study comprises the spatial prediction of soil depth functions for three-dimensional modelling of SOC and bulk density. The spatial prediction of the function coefficients enabled the calculation of two three-dimensional arrays by solving the depth functions for depths from 0 to 50 cm by 5 cm increments. This was used to estimate the SOC stocks in high spatial (5 m) and vertical (5 cm) resolution. The main conclusions of this study are: • The general trend of SOC as visualised by the boxplots (Fig 2) was exponential. However, polynomial depth functions described the soil profiles for SOC with higher accuracy and the logarithmic functions for BD showed better results in spatial modelling. Therefore, we conclude that functions resulting in high accuracies based on the soil profile data may not be the most suitable for spatial modelling, as they may overfit the vertical trend of SOC content.
• The 3D RF models correspond best with the 2.5D counterparts (R 2 up to 0.96). Thus, RF is recommended to predict SOC based on exponential depth functions and bulk density with logarithmic depth functions in high vertical resolution. The 2.5D and 3D predictions of SOC with RF correlated much better, especially when using exponential functions, and lacked accuracy in deeper layers for SOC when modelled based on polynomial functions.
• Comparisons between conventional 2D and 2.5D predictions at the sampled depth and the corresponding depth of the three-dimensional predictions showed that MARS is not suitable for modelling corresponding 2.5D and 3D models, although cross-validation of the individual models showed similar performance in R 2 .
Minor conclusions are: polynomial functions may be an option, when the problem of propagated errors and the ability to generalise in the horizontal domain is investigated further, however, polynomials of any degree have to be used carefully. To overcome these shortcomings, a higher sampling density in the vertical and horizontal domain and in combination with other depth functions, such as equal-area splines [90], should be considered, since exponential functions are not suitable for soil properties that do not increase or decrease continuously.
The 3D approach presented in this study is promising for SOC auditing in various disciplines and especially for decision making regarding climate and land use policies. Future work should focus on sampling design to cover valley positions outside the established plots at site A of BEF-China project. Given the dynamics of SOC stocks, we recommend the analyses of time series data and the expansion of the current database for four-dimensional models.