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-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-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-24]. Hence, models that are based on terrain parameters reproduce displacement and reallocation of soil (i.e. mass movements and soil erosion) and are of particular interest when modelling SOC at catchment scale [25]. Furthermore, terrain can not only be used to estimate or model soil displacement and reallocation, but also as a proxy for environmental covariates, which are not used as predictors, or inaccessible scorpan-factors. For instance, slope and aspect can serve as proxy for microclimate through its influence on local solar insulation [24]. The catchment area can serve as a proxy for soil fertility because of terrain driven water and SOC accumulation [19] and elevation, slope and aspect can act as proxy for parent material, tectonics and periglacial climate through strike and dip of the geological sediments and down-cutting processes [22,23,26].

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 three-dimensional mapping is to consider the vertical dimension as multiple two-dimensional predictions, which can be interpreted in a three-dimensional way [17,27-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-37].

Besides geostatistical frameworks [38,39], different depth functions have been applied for 3D modelling: power, logarithmic [32,40], exponential decay [32,33], polynomial [34,36] and equal-area spline functions [31,41].

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 2nd and 3rd 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.

Material and methods

Study area and sampling design

The BEF-China study sites are artificial biodiversity experiments on property leased and managed by the project partner Institute of Botany, Chinese Academy of Sciences, 20 Nanxincun, Xiangshan, Bejing, 100093, PR China. Field studies did not involve endangered or protected species and no specific permissions for field research were required.

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].

Fig 1

Study area in mainland China with BEF-China plot scheme and indication of sampled plots.

Upper right panel with permission by R. Hijmans; https://gadm.org/.

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 clear-cut 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.

Fig 2

Datasets for SOC and BD used in this study summarized in boxplots.

The boxplots show the variation of the SOC and BD values for each depth increment. SOC and BD samples were taken in five depth increments and 9 cores per plot were bulked (Note that depth increments do not increase linearly). The grey lines show model depth functions (3rd degree polynomial for SOC and natural logarithmic function for BD; see subsection “3D mapping with soil depth functions”).

Bulk density samples (n = 55) were taken in April 2015 with soil sample rings (100 cm3) 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.

Fig 3

Empirical cumulative distribution functions (ECDF) for SOC and BD datasets.

The ECDFs show the locations of the sampling sites in the state space of the elevation (DEM) in metres above sea level (m a.s.l.). The aim is to show the coverage of the DEM feature space by the samples. It can be seen that most samples are located in the mid-range of the elevation values. Therefore, predictions at grid locations which are only sparsely covered by the samples (i.e. locations close to the minimum and maximum values of the DEM) may be less accurate. The minimum, median and maximum values of both datasets (DEM and sampling locations) are shown with vertical lines (dashed grey: DEM, dashed black: sampling locations) to compare the full range of the respective feature spaces.

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].

Table 1

Terrain attributes used for SOC and bulk density modelling.

	Covariates	Method	Author(s)
Local	Slope and aspect	Fitted 2nd degree polynomial	[54]
		Fitted 3rd degree polynomial	[55]
		Least squares fitted plane	[56]
		Maximum triangle slope	[57]
		Fitted 2nd degree polynomial	[58]
	Plan, profile, longitudinal, tangential and flowline curvature	Fitted 2nd degree polynomial	[54]
		Fitted 3rd degree polynomial	[55]
		Fitted 2nd degree polynomial	[58]
	Vertical distance to channel network		[53]
	Sky visibility, sky view factor, direct and diffusive insolation		[59]
Regional	Catchment area	Top-down	[60]
		Recursive
Combined	Topographic Wetness Index (TWI)	Any combination of slope and catchment area	[61]
	Slope length and steepness factor (LS-Factor)	Any combination of slope and catchment area	[61,62]

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 n