###### Abstract
Thus, it becomes far from straightforward to ensure consistency between SOC stored in single layers and the stocks of the entire soil colum at every single site in the resulting maps. Analogous problems arise for the resulting depth gradients of SOC, which might adopt artificial shapes. Furthermore, fitting of individual models for respective soil layers forces the user to deal with several different estimates of variable importance. Except for soil layers with a priori well-known differences regarding the main processes for carbon enrichment, such as the organic layer and mineral soil, it will hardly be possible to logically interpret the influence of covariates in terms of soil forming factors. This may hinder the optimal selection of predictors and detection of spurious covariates to be excluded. Last but not least, individual regionalisation of depth layers may enhance the propagation of potentially large sampling errors caused by incorrect separation. Such errors especially result from mingling of the organic layer and mineral soil [29]. The usage of depth functions to describe the vertical variation of soil properties has a long tradition in soil science (see [30] for a detailed overview). The application of soil attribute depth functions in three-dimensional spatial prediction of soil properties was already proposed by [PERSON] et al. [31]. It has also already been shown that the concept of depth functions is beneficial in mapping SOC distribution throughout the soil profile: Equal-area quadratic smoothing splines were used to derive SOC for homogeneous depth increments from legacy profile data as a data basis for the development of artificial neural networks and regression models [21; 32; 33]. To obtain maps of continuous depth functions, [PERSON] et al. [32] utilised equal-area quadratic splines for interpolation through predicted SOC contents of discrete soil layers after regionalisation. To calculate weighted means of SOC for mapping units from available legacy soil data, [PERSON] et al. [34] applied equal-area quadratic smoothing splines. From a more conceptual perspective, depth functions are used for two reasons within 2.5D soil mapping frameworks [35]: The first reason is the harmonisation of soil observations. Equal-area quadratic smoothing splines are then often fitted to depth interval values after mapping to derive continuous depth functions for each pixel mapped. The equal-area smoothing splines of [PERSON] et al. (2018) can be considered as the state of the art for these two applications for one reason in particular: The interpolation process preserves the masses originally measured or predicted on a depth interval basis (especially for small values of \(\lambda\)). Direct utilisations of depth functions during the mapping process itself are given by [PERSON] et al. (2018) and [PERSON] et al. (2018). [PERSON] et al. (2018) used ordinary kriging to interpolate parameters of exponential functions fitted to single-point observations beforehand. [PERSON] et al. (2018) predicted parameters of soil-type-specific depth functions of SOC using cokriging. By using relative depth gradients, we try to strengthen the interpretability of SCORPAN factors in three-dimensional mapping of SOC in forest soils. Our further objectives were to ensure the consistency between solum stocks and stocks within single depth layers, to avoid the construction of artificial depth gradients, and to lessen the propagation of sample separation errors. The basic conceptual hypothesis of the method is that distinguishing between covariates that control horizontal and vertical distribution of SOC simplifies variable selection during model development. Even though increasing interpretability may decrease model performance measures, we believe that a good understanding of covariates in terms of soil forming factors is valuable for the mapping results.