Predictive mapping of soil organic carbon in Northeast Algeria

Catena 190 (2020) 104539

Contents lists available at ScienceDirect

Catena journal homepage: www.elsevier.com/locate/catena

Predictive mapping of soil organic carbon in Northeast Algeria a,⁎

b

a

c

Sana Boubehziz , Kamel Khanchoul , Mohamed Benslama , Abdelraouf Benslama , Alessandro Marchettid, Rosa Francavigliad, Chiara Piccinid

T

a

Laboratory of Soil and Sustainable Development, Department of Biology, Badji Mokhtar University-Annaba, P.O. Box 12, 23000 Annaba, Algeria Laboratory of Soil and Sustainable Development, Department of Geology, Badji Mokhtar University-Annaba, P.O. Box 12, 23000 Annaba, Algeria Laboratory of Mathematics and Applied Sciences, University of Ghardaia, P.O. Box 455, 47000 Ghardaia, Algeria d Council for Agricultural Research and Economics, Research Centre for Agriculture and Environment, via della Navicella 2–4, 00184 Rome, Italy b c

ARTICLE INFO

ABSTRACT

Keywords: Ordinary kriging Regression kriging Semiarid soils Spatial variability Land use

Soils of semi-arid climates undergo organic carbon loss, which in turn affects their agricultural potential. Geostatistics is often used as an interpolation tool to thoroughly describe SOC spatial distribution. To focus on soil organic carbon (SOC) depletion, the Tiffech watershed (Northeast of Algeria), an economically important agricultural area, was chosen due to intensive agricultural practices, decline of forests and occurrence of erosion. The present study aimed to predict the spatial variation of SOC in Tiffech watershed using geostatistics and a Geographical Information System (GIS) software, comparing the performance of two geostatistical methods—Ordinary Kriging (OK) and Regression Kriging (RK)—also assessing the role of auxiliary variables in improving the prediction accuracy and highlighting the role of land cover in SOC storage. The SOC content in Tiffech soils was determined on 42 soil samples from the surface layer (0–10 cm) collected all over the study area and the results were used to estimate SOC density in non-sampled locations. The prediction efficiency of the two methods was evaluated by calculating the Mean Error (ME), the Root Mean Square Error (RMSE) and the Root Mean Square Standardized Error (RMSSE). The interpolation results showed that SOC distribution in the study area was correlated to the topography, the clay index, and general landscape features. SOC content increased northwards in the area, ranging from 0.53 to 6.9 kg·m−2 in relation to land use. As expected, maps figured a good conservation status of SOC stocks in areas with dense vegetation; conversely poor SOC contents were estimated where land degradation factors take place. Cross-validation results showed an outperformance in the interpolation accuracy of RK on OK after the introduction of environmental variables, with an RMSE value of 0.02 versus 0.81. This suggests a higher efficiency of RK in predicting SOC content across the Tiffech area in comparison with OK, confirming that introducing some auxiliary data correlated to the target variable in SOC estimation, considerably improved the interpolation accuracy.

1. Introduction Despite carbon is widespread through all environmental compartments, i.e. living organisms, soil, water, atmosphere, and minerals (Hartemink and McSweeney, 2014), soil carbon in the first meter depth scores the highest amount, with 1500 Pg estimated for the whole earth (Kutsch et al., 2009), making soil carbon the major C pool and top soil the main C sink. Due to the key role of organic carbon in regulating soil chemical, physical and biological properties, soil organic carbon (SOC) is adopted as an indicator of land degradation (Marchetti et al, 2010). Moreover, in agricultural soils it represents an indicator of soil fertility and agricultural production potential (Roose et al., 2005; Bleuler et al.,

⁎

2017). Organic carbon sequestration in soils has gained growing attention as a means to reduce atmospheric CO2 and, in perspective, global warming (Lal, 1999; Pan et al., 2003; Johnson et al., 2006). Previous studies, e.g.Tondoh et al. (2016), showed that sound land management of agricultural soils in semi-arid climate has a potential for organic carbon sequestration. The SOC originates from plant residues and organic fertilizers, being SOC stocks depending on soil decomposition rate, which in turn is driven by microorganisms and environmental factors such as land use, soil type (Piccini et al., 2014), temperature, topography, etc. (Johnson et al., 2006; Cooper et al., 2011; Li et al., 2018). Elevation and temperature also may have an impact on local SOC content (Teng et al., 2017).

Corresponding author. E-mail addresses: [email protected], [email protected] (S. Boubehziz).

https://doi.org/10.1016/j.catena.2020.104539 Received 7 August 2019; Received in revised form 28 February 2020; Accepted 1 March 2020 0341-8162/ © 2020 Elsevier B.V. All rights reserved.

Catena 190 (2020) 104539

S. Boubehziz, et al.

SOC maps, which are a good tool for soil quality evaluation, are often affected by a limited knowledge of the spatial variability of soil types and properties. As a consequence, several methods were developed by researchers to perform the estimation of soil parameters, aiming at reporting the spatial variation of soil properties with acceptable accuracy. The first law of geography given by Tobler in 1970 “Everything is related to everything else, but near things are more related than distant things”, is the basis of spatial interpolation and is still used to estimate the natural variation of soil at unsampled locations (Castrignanò, 2011). In the past 15 years, research has focused on probabilistic (stochastic) methods such as geostatistics and Digital Soil Mapping (DSM) to allow predictive mapping accounting for soil heterogeneity and spatial variability; in DSM soil parameters are interpolated from observed points over an area—using a reduced number of samples—with a known error (Hengl, 2007, Long et al., 2014). Ordinary Kriging (OK) is one of the most commonly used geostatistical methods (Piccini et al., 2014) consisting of a weighted average of the neighborhood data values (Webster and Oliver, 2007). It performs the estimation of a target variable (e.g. a soil property) from point data (Webster and Burgess, 1983) based on their spatial autocorrelation, which determines the coefficient of linear prediction (Odeh et al., 1994). OK is effective in interpolating point data and in analyzing the spatial variability with a known error using a minimum number of samples (Long et al, 2014). However, the accuracy of such method is someway limited as it does not account for other variables—such as environmental variables—that would otherwise improve the interpolation results (Odeh et al., 1995). Therefore, spatial prediction can be improved by introducing attributes like terrain indices, soil texture, etc. (Zhang et al., 2012) as in the case of hybrid methods (McBratney et al., 2000). Regression Kriging (RK), the most used hybrid method (Keskin and Grunwald, 2018), is known as a BLUP (Best Linear Unbiased Prediction) technique (Sun et al., 2012) and is defined by Hengl (2007) as a spatial interpolation technique combining auxiliary variables, derived from Digital Elevation Model (DEM), remote sensing imagery and thematic maps, with simple kriging of regression residuals. According to Odeh et al. (1995), RK provides a cheap and reliable way of predicting land properties, and the spatial prediction is generally improved when soil variables are

estimated into a relatively finer gridded DEM. The introduction of environmental variables having a direct relationship with the target variable could enhance the interpolation results, however the advantage of using RK is debatable. This method in fact could be time consuming in relation to the limited improvement of results and other studies also found the introduction of auxiliary variables in soil organic matter interpolation to not substantially enhance the interpolation results, in comparison with other methods (Piccini et al., 2014). Recent studies about SOC mapping (Arrouays et al., 2017; Song et al., 2017; Wang et al., 2017; Yorulmaz et al., 2017; Costa et al., 2018; Bhunia et al., 2017) and SOC storage in soil (Ma et al., 2015; Chartin et al., 2016; Yang et al., 2016b) across different scales (Kumar, 2013, O’Rourke et al., 2015, Wang et al., 2017) were performed using a variety of statistical methods. Among them, the k-nearest neighbours method known as K-NN, used by Zolfaghari et al. (2016), handles only the information from the neighbours. Boosted regression tree, a simple linear combination, selects a random subset from the complete dataset, requires a large number of samples and is more complicated than kriging methods, needing at least two predictor variables. Conversely, OK uses the measured variable as the unique predictor (Yang et al., 2016a). Many studies have suggested that RK has higher accuracy in SOC estimation than OK (Zhang et al., 2012; Ling et al., 2014) since this method utilizes secondary information that could considerably reduce the Root Mean Square Error (RMSE) (Simbahanet al., 2006). Zhou et al. (2016) suggested that RK is the most accurate method, showing advantages in SOC prediction in comparison with other approaches such as Multiple Linear Regression (MLR) and Inverse Distance Weighting (IDW). On the other hand, Sherpa et al. (2016) illustrated that OK could represent an optimal strategy for SOC stock estimation in several soil agroecosystems. Considered the above, the objectives of this study are: (i) to estimate the SOC stock across the study area and calculate its storage under different land uses, (ii) to assess the spatial variability of SOC at the watershed scale under a semi-arid climate using geostatistics, (iii) to examine the capability of ancillary information to improve SOC estimation accuracy by comparing two different methods (OK and RK), and (iv) to explore land use strategies to mitigate land degradation associated with SOC depletion.

Fig. 1. The study area - southeast of Souk Ahras region, Algeria. 2

Catena 190 (2020) 104539

S. Boubehziz, et al.

• Climate

2. Materials and methods 2.1. Study area

The study area represents a transition from the Mediterranean to desert climate and is classified as a semi-arid area according to the Emberger climatic diagram and in “Csa” class (warm temperate–steppe–hot summer) according to the updated World Map of the Köppen-Geiger Climate Classification (Kottek et al., 2006). The mean annual temperature and rainfall are about 15 °C and 322 mm respectively, with autumn and winter as rainy seasons.

The Tiffech watershed, covering about 89 km2, is located in the Southeast of Souk Ahras region, Algeria (Fig. 1). Mountains surround the region at North, West, and Southeast while the Tiffech river flows along the plain located in the middle and South of the area. The main tributary of this river is Oued Hamimine. The elevation in Tiffech watershed ranges from 831 to 1103 m a.s.l.

2.3. Soil sampling and SOC analysis 2.2. Soil properties

During a soil sampling survey performed in February 2017, a set of soil samples (Fig. 1) was collected from topsoil (about 10 cm depth) with a density of 1/2 sample per km2. The sampling grid was drawn on a Geographic Information System (GIS) spatial database; 275 sampling points were selected. Aiming to reduce the number of samples, the area was split into more homogenous zones by overlapping two maps of land use and soil type; eight homogenous sub-layers were obtained, from which 42 samples were randomly selected. In Table 2 the resulting sublayers with the respective number of samples are reported. Soil samples were air-dried for at least 48 h, crushed and sieved on 2 mm mesh to remove gravel and roots. SOC content was measured according to the modified Walkley-Black method (Nelson and Sommers, 1982). SOC density (kg·m−2) was estimated based on SOC content (g g−1), soil Bulk Density (BD, g cm−3), and soil depth (cm) by using the following equation:

Table 1 shows the main soil properties of the study area used as predictors of the target variable:

• Lithology As stated by the National Office of Geology (1978), the Tiffech basin is made up mainly of sedimentary deposits. The plain consists of deposits and alluvium—gravels, sand, and silt—of late Quaternary age. The eastern hilly area is mainly a complex of conglomerates, sandstones, clays, and Pliocene marls. The western part is characterized by hills of limestone and marls of Cretaceous age, and clays, gypsums, and sandstones of Triassic age.

• Land use

SOC density = SOC × BD × soil depth × 10

The land use map (Fig. 2)—obtained by digitalizing satellite images—indicates that about 87% of the area is farmland, mostly cropland and pasture. Wheat is the most widespread crop, also in rotation with horticultural crops, representing one of the most important products in Algeria (USDA Foreign Agricultural Service, 2018). In the North-Western part of the Tiffech watershed, about 7.8% of the total area, a huge pine forest exploited for timber burned during summer 1956. Another 3.7% of woodland in South-West is a Pinus halepensis forest.

Soil BD was estimated based on soil texture using a soil parameter calculator (Arpa Sardegna, 2018). 2.4. Predictive mapping OK and RK techniques were used to predict the spatial variability of SOC, and to test the ability of auxiliary variables to improve the predictions.

• Soil type

2.4.1. Mapping by OK OK is a method to estimate the values of a target variable in nonmeasured locations Zʹ(x0), using data in the neighborhood Z(xi), (i = 1, 2… n) (Wackernagel, 1994), as follows:

Two main soil types are present in the study area: in the West, soils are developed on limestone rocks, which enriched the soil in calcium carbonates by alteration; in the East, soils contain a significant amount of secondary carbonates in the form of calcic horizons (Durand and Barbut, 1948).

n

Z '(x 0 ) =

Type

Range

DEM

Elevation Slope Aspect Solar Radiation Plan Curvature Profile Curvature Topographic Wetness Index Flow Accumulation Valley Distance Valley Depth NDVI CI NDSI Lithology Soil Type Land Use

[831,1103 m] [0.014,75.04%] [0,360°] [1.3e−0.05,6.3] [−0.009,0.01] [−0.009,0.01] [3.46,22.6] [900,6.92e007] [0,151.3] [0,108] [−0.24,−0.003] [0,0.15] [0,0.5] 10 layers 2 layers 6 layers

Landsat 8 TM Thematic maps

i · Z (x i ) i=1

Table 1 Terrain attributes. Source

(1)

(2)

where i is the kriging weight. According to Webster and Oliver (2007), to ensure that the estimate is unbiased, the weights are made to sum to 1: n i

=1

(3)

i=1

and the expected prediction error is

E [z (x 0 )

Z (x 0)] = 0

(4)

OK uses the variogram as a tool to quantify spatial correlation and similarity of measured data, by measuring the variability between the point Z(xi) and its neighbours Z(xi + h) at a given distance h. To estimate the weights in an objective way, so that the weights reflect the true spatial autocorrelation structure, the variogram is obtained by plotting the half of the average squared difference between points separated by the distance h (lag), according to the following equation:

3

Catena 190 (2020) 104539

S. Boubehziz, et al.

Fig. 2. Land use map.

(h ) =

1 2N (h)

distribution of the target variable (Webster, 2001), the frequency distribution of SOC data was verified by plotting a histogram and by performing a Kolmogorov-Smirnov (K-S) test1 at 0.05 significance level. For SOC estimation in the Tiffech watershed, a cloud variogram was

N (h)

(Z (x i ) i=1

Z (xi + h))2

(5)

where Z (xi ) is the value of a target variable at some sampled location and Z (xi + h) is the value of the neighbour at distance h. Supposing n point observations, there will be n·(n 1) 2 pairs for which a semivariance can be calculated (Hengl, 2009). Since classical geostatistical techniques require a normal

1 K-S test is a nonparametric test of the equality of one-dimensional probability distributions that can be used to compare one or two samples with a reference probability distribution.

4

Catena 190 (2020) 104539

S. Boubehziz, et al.

2.4.2.1. DEM derivatives. Eleven DEM derivatives were used:

Table 2 Homogenous sublayers and sampling density. Soil type

Land use

Samples

Area (km2)

Limestone soil

Cultivated Forest Degraded forest Urban Naked soil Cultivated Degraded forest Meadow

8 4 4 1 1 22 1 1

24.16 3.25 5.95 0.32 1.03 53.29 1.01 0.19

Calcic soil

• Elevation (ELEV), ranging from 831 m a.s.l. in the watershed outlet to 1103 m a.s.l., in the mountains surrounding the study area. • Slope gradient, ranging from 0.014% to 75%. • Aspect (ASP), ranging from 0 to 360°. The eastern area is mainly • •

plotted with all squared differences to describe the spatial variability between data points separated by a lag distance h. By averaging the values to a standard distance, an empirical variogram is obtained. The variogram increases with distance indicating a decrease of correlation with the distance. The distance beyond which there is no more autocorrelation among variables is the range of the variogram. At some distance, the covariance tends to zero and the variogram becomes constant and flat out into a sill. In the case of a smooth curve fitted to the observed semi-variance not passing through the origin and cutting the ordinate at a positive value, this value is the nugget variance (Webster and Burgess, 1983). The nugget/sill ratio can be used as an indicator of spatial dependence: if the ratio is ≤25% the variable is considered to be strongly spatially dependent, if the ratio is between 25% and 75% the variable is considered to be moderately spatially dependent, and if the ratio is ≥75% the variable is considered to be weakly spatially dependent (Cambardella et al., 1994). Point weights for OK are calculated by fitting a theoretical model to the empirical variogram, considering the model with the smallest error by least squares criterion as the best one. Before applying OK, the dataset was checked for anisotropy: the variables showed the same spatial autocorrelation structure in all directions, meaning that the semivariance depends only on the distances among the sampling points (Burgos et al., 2006). Thus, an omnidirectional (isotropic) variogram is plotted considering data pairs separated by a lag distance of 1330 m. Minimizing the sum of the squared deviations among the experimental and theoretical values (Lark, 2000), a suitable theoretical model is fitted to the variogram to create a continuous surface. Several models were tested in order to choose the best fit (exponential, linear, Gaussian, etc.). The spherical model gave the lowest estimation error and was chosen as the best one:

lag Spherical model = nugget + sill 1.5 range

lag 0.5 range

exposed to the Northeast, while the western area is mainly exposed to South and Southwest, following the main river in the watershed. Plan and profile curvature (PLANC and PROFC), obtained from the second derivative of the maximum slope direction and the perpendicular one, respectively, (Olaya, 2009); Topographic Wetness Index (TWI), a parameter correlated to the topography and to the water movement, expressing the spatial distribution of soil moisture, (7)

TWI = ln (As tan )

• • • •

where As is specific catchment area and β is slope (Beven and Kirkby, 1979). Solar radiation. Flow accumulation (FLOWACC). Valley distance and valley depth (VDIST and VDEPTH). Convergence index (CONVI).

2.4.2.2. Landsat 8 TM derivatives. The indices extracted from Landsat 8 TM during the study (March 2017, cloud cover 0%) were:

• Normalized Difference Vegetation Index (NDVI), which could distinguish the presence of green vegetation, or not, calculated according to Colwell (1974):

NIR R NIR + R

NDVI =

(8)

where NIR is Near Infrared (band 5) and R is red (band 4);

• Clay Index (CI), correlated with the soil clay content, calculated according to Hengl (2009):

CI =

•

MIR MIR2

where MIR is Mid Infrared (band 5) and MIR2 is Mid Infrared (band 7); Normalized Difference Soil Index (NDSI), enhancing soil information and largely suppressing impervious areas and vegetation values, calculated according to Deng et al (2015):

NDSI =

3

(6)

(9)

(MIR2 G ) (MIR2 + G )

(10)

where MIR2 is Mid Infrared (band 7), and G is Green (band 3).

• Grain Size Index (GSI), correlated with the fine sand content of the soil:

2.4.2. Mapping by RK Differently from OK, RK was used to perform an estimation of data correlated with environmental variables; terrain attributes are the most commonly used auxiliary variables in SOC mapping (Behrens et al., 2010; Yang et al., 2016b). In this study, terrain attributes were derived from: (i) a DEM (Shuttle Radar Topography Mission ver. 2, SRTM) with 30 m ground resolution, (ii) satellite imagery from Landsat 8 TM, and (iii) thematic maps of geology, lithology and land use of the area. All these variables were used as predictors for SOC estimation (Table 1).

GSI =

(R B) (R + G + B)

(11)

where R is Red (band 4), B is Blue (band 2), and G is Green (band 3) (Xiao et al., 2006). A series of raster files were generated for all extracted indexes. Stepwise regression with a significance level of 0.05 was used to select the smallest set of predictors giving the best linear regression results

Table 3 Descriptive statistics of soil samples. Min. (kg⋅m−2)

Median (kg⋅m−2)

Mean (kg⋅m−2)

Max. (kg⋅m−2)

Skewness

Kurtosis

Std. dev (kg⋅m−2)

CV (%)

Variance

K-S

0.53

2.72

3.17

14.0

2.36

6.60

2.55

80.56

6.51

0.006

5

Catena 190 (2020) 104539

S. Boubehziz, et al.

Table 4 Descriptive statistics without outliers. Min. (kg⋅m−2)

Median (kg⋅m−2)

Mean (kg⋅m−2)

Max. (kg⋅m−2)

Skewness

Kurtosis

Std. dev (kg⋅m−2)

CV (%)

Variance

K-S

0.53

2.69

2.62

7.8

1.04

2.64

1.40

53.72

1.97

0.60

– The Root Mean Square Standardized Error (RMSSE)

Table 5 Parameters of the optimized variogram model. Variable

Nugget

Sill

Range

IGF

SOC

0.812

3.271

27971.73

0.0005

RMSSE =

2.5. Validation The precision of prediction was evaluated by comparing the estimated values z^ (x i ) with actual observations at validation points Z (xi ). The accuracies of the two methods were evaluated by leave-one-out cross-validation, with three statistical indices quantifying the amount of difference in the model simulation deviating from the observation, presented in the unit of the observation:

(Z (x i )

Z (x i ))

(12)

– Root Mean Square Prediction Error (RMSE)

RMSE =

N i=1

Z (x i ))2

(Z (x i ) N

(x i ))]2

N

(14)

Descriptive statistics of sampling points data including mean, minimum and maximum values, and coefficient of variation (CV) are presented in Table 3. The K-S test indicated that data were not normally distributed (P < 0.05); accordingly, the histogram indicated a positive skewness in the distribution of the data (2.36). The high CV value (80.56%) suggested some outliers to likely affect the measures of central tendency. Three outliers were identified on boxplot and then deleted, afterwards the dataset fitted the normal distribution with a K-S Pvalue > 0.25 and was suitable for geostatistical elaboration. Descriptive statistics of the remaining 39 samples, used for geostatistical application, are presented in Table 4. Mean and median values were almost equal, with no more outlier effect, although kurtosis and skewness indicated a high spatial variability in SOC content.

N j=1

Z (x i )

3. Results and discussion

– Mean prediction error (ME)

1 · N

^ (x ) [(Z i

where N is the number of validation points and Z (x i ) is the prediction standard error in location x i . The extraction of ancillary information from DEM and Landsat 8 was performed using SAGA 5.0.0 (SAGA User Group Association, 2017) and QGIS 2.14.14 (QGIS Development Team, 2017); statistical analysis of SOC, distribution models and error estimation for both OK and RK were performed with the statistical software RStudio 1.0.143 (R Development Core Team, 2016) and its packages sp, gstat, Rcmdr.

(Lima et al., 2016). The correlation among the target variable and the 18 auxiliary variables described above was first verified by means of a Principal Component Analysis (PCA).

ME =

N i=1

(13)

Fig. 3. Experimental variograms of SOC: (a) variogram of measured data, (b) variogram of regression residuals. 6

Catena 190 (2020) 104539

ELEV: Elevation; Slope: Slope gradient; ASP: Aspect; SOLAR: Solar radiation; PROFC: Profile Curvature; PLANC: Planc Curvature; TWI; Topographic Wetness Index; FLOWACC: Flow Accumulation; VDIST: Valley Distance; VDEPTH: Valley Depth; CONVI: Convergence Index; NDVI: Normalized Difference Vegetation Index; NDSI: Normalized Difference Soil Index; GSI: Grain Size Index; CI: Clay Index; LU: Land Use; Pedo: soil type; LITHO: Lithology.

**

*

***

23.8 41.03 53.14 63.87 71.29 77.64 83.03 87.15 90.76 93.25 95.37 96.52 97.45 98.31 99.13 99.56 99.87 100.00 Adjusted R2 0.8398 23.8 17.23 12.1 10.73 7.43 6.35 5.39 4.12 3.51 2.59 2.12 1.15 0.92 0.86 0.82 0.42 0.30 0.13 −0.579 0.002 0.485 −0.451 −0.047 0.031 0.043 −0.142 0.150 −0.377 −0.029 −0.102 −0.041 −0.051 −0.106 0.091 −0.031 0.019 −0.51 −0.06 −0.42 0.42 −0.06 −0.12 −0.11 0.19 0.26 −0.13 0.13 −0.28 0.30 0.02 −0.13 0.16 −0.05 −0.03 0.06 0.2 0.2 −0.02 0.20 0.05 0.20 0.06 0.09 0.32 −0.12 −0.66 0.27 −0.33 0.20 −0.15 −0.02 0.14 −0.19 0.23 0.04 0.12 0.11 −0.07 0.09 −0.09 0.57 0.43 −0.12 0.25 −0.34 0.05 0.15 0.24 0.27 0.11 −0.107 0.083 −0.086 0.372 0.031 0.442 0.375 −0.298 −0.236 −0.201 −0.129 −0.068 −0.265 0.008 0.193 0.247 −0.310 0.167 ELEV SLOPE ASP SOLAR PROFC PLANC TWI FLOWACC VDIST VDEPTH CONVI NDVI NDSI GSI CI LU Pedo LITHO

−0.41 −0.31 0.16 0.14 −0.26 −0.19 0.30 0.20 −0.35 0.37 −0.23 0.07 −0.06 0.20 0.10 −0.31 0.01 −0.05

0.05 0.05 0.12 0.09 0.09 0.12 −0.07 −0.02 −0.15 0.14 0.13 −0.49 −0.40 0.36 −0.44 0.10 0.26 −0.30

−0.117 −0.188 −0.542 −0.576 0.269 0.217 0.277 0.280 −0.024 0.104 0.003 −0.046 −0.136 0.017 0.050 0.109 0.012 −0.063

0.3 −0.15 −0.11 −0.07 −0.33 −0.42 0.19 0.09 0.19 −0.20 −0.42 −0.24 −0.27 −0.05 −0.22 0.16 −0.10 0.27

−0.12 0.07 −0.23 −0.18 −0.23 −0.15 −0.34 −0.37 0.05 0.25 0.13 −0.17 −0.32 −0.09 0.17 −0.23 −0.49 −0.17

−0.109 −0.25 −0.231 −0.078 −0.117 0.038 −0.188 −0.519 −0.232 0.090 −0.073 −0.102 0.105 −0.159 −0.009 0.039 0.551 0.357

0.08 −0.2 −0.01 0.04 −0.05 0.26 0.04 0.03 0.38 −0.19 0.24 −0.12 −0.12 0.38 0.13 −0.56 0.07 0.38

−0.063 0.362 0.014 −0.106 −0.429 0.170 0.008 0.297 −0.177 0.232 0.378 0.076 −0.042 −0.109 −0.267 0.121 −0.063 0.464

0.076 −0.572 0.187 0.168 −0.123 0.242 0.007 0.129 0.180 0.080 0.246 0.022 −0.219 −0.573 −0.041 0.087 0.016 −0.169

0.034 −0.019 0.05 0.003 0.079 −0.509 0.258 −0.024 −0.170 −0.157 0.591 −0.113 −0.143 0.026 0.410 0.185 0.150 0.041

−0.007 −0.34 0.137 0.016 0.487 −0.181 −0.031 −0.165 −0.019 0.261 0.133 0.080 0.094 0.142 −0.342 0.125 −0.406 0.385

−0.12 0.25 −0.18 0.12 0.26 −0.18 0.29 −0.12 −0.05 −0.15 0.07 0.16 −0.20 −0.40 −0.41 −0.49 0.10 −0.03

0.14 −0.03 −0.05 −0.10 −0.32 0.05 0.53 −0.41 0.23 0.16 0.20 0.05 0.38 0.15 −0.21 0.05 −0.04 −0.28

PC17 PC16 PC10 PC6 PC5 PC4 PC3 PC2 PC1

Table 6 Principal component coefficients and significance for prediction.

PC7

PC8

PC9

PC11

PC12

PC13

PC14

PC15

PC18

VAR%

CUMVAR%

SOC

S. Boubehziz, et al.

Table 7 Parameters of the optimized variogram model of the residual. Variable

Nugget

Sill

Range

IGF

SOC

0.00

1.329

1239.56

0.001

3.1. Ordinary kriging application The best fitted model parameters are reported in Table 5 with an Indicative Goodness of Fit (IGF) value; being this value very close to zero, the fitted model is adequate. The chosen model had a range of 6731.51 m (Fig. 3a); the presence of a high nugget value refers to measurement error, sampling error, inter-sample error (Burgos et al., 2006). The obtained variogram had a nugget/sill ratio of 24%, indicating that the variable was strongly spatially dependent. 3.2. Regression kriging application Table 6 shows the matrix of transformation coefficients, calculated from the covariance matrix. PC1 component was the most significant for SOC prediction, and was mainly correlated with DEM derived parameters (ELEV, Slope, TWI, VDIST, VDEPTH), and with the land use of the area. The high adjusted R2 value (0.84) for the linear regression indicated a highly significant correlation of the target variable with the principal components. For RK application, elevation, aspect, and CI gave the highest correlation to SOC, indicating that topography, slope direction, and soil texture mostly affect the SOC content and distribution. The high value of R2 (0.97) indicated a strong correlation between predictors and SOC content; these predictors were used as explanatory variables in SOC interpolation and used to develop a SOC prediction model. Fig. 3b shows the experimental variogram for regression residuals, and Table 7 shows the parameters of the fitted model. 3.3. SOC distribution across the Tiffech watershed Figs. 4 and 5 show the SOC distribution maps across the study area according to OK and RK interpolation methods. Although both maps showed similar pattern of SOC distribution, RK map figured a patchy, more detailed, SOC distribution than OK map. The estimated average SOC content in the area was about 3.1 kg·m−2 for OK and 3.65 for RK. The two maps showed that SOC content increased from the South and the South-East to the North-West reaching the highest values (more than 4 kg·m−2) in 2% and 6% of the total area estimated by OK and RK respectively. The two interpolated maps showed a correlation between SOC density and the type of land use in the area, with SOC stock changing on varying land use: cultivated soils were found to hold about 2.50 kg·m−2 of SOC, i.e. almost half of SOC quantity in soils under forest (4.70 kg·m−2). Soils under degraded forest showed higher SOC contents compared to cultivated soils (+25%), with a stock of 3.33 kg·m−2. Results from both interpolated maps confirm land cover to have a remarkable impact on SOC distribution, indicating a similar variation in SOC between croplands and forest areas. We speculate the high SOC content in the North-Western area to be due to the presence of a degraded forest supplying a considerable quantity of organic carbon to soil. Differently, SOC content reached the lowest value (0 to 1.5 kg·m−2) in the watershed outlet, likely due to the semi-arid climate boost of SOC degradation, as well as the highly flooded topography and the intensive agricultural practices. As a consequence of the impact of these factors on SOC content, also the productive capacity of the area (wheat production) and to a greater extent, the country economy are affected in turn. Forestland showed higher SOC content than cultivated 7

Catena 190 (2020) 104539

S. Boubehziz, et al.

Fig. 4. Estimated map of soil organic carbon by OK.

land, probably due to the accumulation on topsoil of conifer debris (cone and leaf litter), which is known to be resistant to decomposition (Soleimani et al., 2017). Besides the fact that forests provide soil protection and dead wood and dead roots accumulate on the soil over time, the canopy reduces runoff and tree roots hinder soil erosion. Conversely, cultivated lands presented a lower SOC content that could be due to the low input of organic materials, which is usually harvested and removed at regular intervals. Farming practices in this area contribute to removing organic materials and harvest residues, that reduce considerably SOC content (Mingjun et al., 2017). Moreover, cultivated lands located at high elevation are more vulnerable to SOC loss compared to those in lowlands, due to the highest exposure to soil degradation factors and non conservative agricultural practices. An intermediate area located between forestland and cropland had a SOC stock lower than forest soils and higher than cultivated soil; this area is a burned forest which still maintains a considerable SOC stock due to its past use. The results showed also a positive correlation of SOC with elevation: the highest SOC density (5 kg·m−2) was measured in the highest mountains in the West (more than 1000 m a.s.l.) and this density decreased with the elevation to reach the lowest value in the watershed outlet (831 m a.s.l). Previous researches have highlighted the strong relationship of SOC content with the topography of the area (Mingjun et al., 2017). Such a relationship can be explained by the low temperature and the high moisture level at high elevations, affecting the rate of SOM decomposition. Moreover, the presence of forests in these areas increases the amount of SOM, preventing soil from erosion. The SOC distribution followed also the variation in clay content, due to its role in soil aggregation and in the formation of humo-mineral complexes (Augustin and Cihacek, 2016; Banerjee et al., 2018).

3.4. Models accuracy A cross-validation on OK and RK estimations was performed to define the best model for predicting SOC distribution in the Tiffech watershed. Measured prediction errors for the two methods are presented in Table 8 where the performance of RK in predicting SOC was higher than OK. The OK map represented the spatial variability with less detail, giving more importance to large distances variability, while the RK map represented more details at short distances with a pattern consistent with the local relief morphology. The MEs of the two interpolation methods were convergent, but the RMSE value was lower using RK. The standard errors showed similar patterns for the two models: as shown in Fig. 6, the error was smaller in the center and increased toward the edges of the watershed, where sampling density was lower. Nonetheless, RK presented lower errors than OK. The results of this study showed that SOC interpolation by OK, based only on measured data, could be not enough accurate, failing in well representing the actual SOC status. The addition of other related variables to the interpolation model, as in RK, could enhance significantly the spatial estimation and hence provide more accurate results. In this work, the results of the comparison of OK and RK performances in predicting SOM distribution, confirm the findings of Hengl et al. (2004) and Piccini et al. (2014). Other studies, e.g. Odeh et al. (1995), found that the performance of OK in interpolating soil properties is poor and the method was significantly outperformed by RK. We speculate that differences in data sources, scale of prediction, types of ancillary variables and their correlation strength with the

8

Catena 190 (2020) 104539

S. Boubehziz, et al.

Fig. 5. Estimated map of soil organic carbon by RK. Table 8 Estimation error values for ordinqry and regression kriging.

SOC

Ordinary kriging ME 0.0007559385

RMSE 0.810934

RMNSE 0.7737199

Regression kriging ME 0.002859975

Fig. 6. Standard error of SOC prediction by a) OK method and b) RK method. 9

RMSE 0.01691507

RMNSE 0.09601438

Catena 190 (2020) 104539

S. Boubehziz, et al.

target variable (Wang et al., 2017; Hengl et al., 2004) may explain the different performance of OK and RK in interpolating SOC content from one study to another. In our study instead, the correlation among ancillary variables and the target variable was as strong as to provide evidence of a reliable interpolation of SOC with a smaller error in comparison to the results by the OK method.

March 2017. Augustin, C., Cihacek, L., 2016. Relationships between soil carbon and soil texture in the northern Great Plains. Soil Sci. 181, 386–392. https://doi.org/10.1097/SS. 0000000000000173. Banerjee, K., Bal, G., Mitra, A., 2018. How soil texture affects the organic carbon load in the mangrove ecosystem? A case study from Bhitarkanika, Odisha. In: Singh V., Yadav S., Yadava R. (Eds.), Environmental Pollution. Water Science and Technology Library, vol. 77. Springer, Singapore, pp. 329–341. https://doi.org/10.1007/978981-10-5792-2_27. Behrens, T., Zhu, A.-X., Schmidt, K., Scholten, T., 2010. Multi-scale digital terrain analysis and feature selection for digital soil mapping. Geoderma 155, 175–185. https://doi. org/10.1016/j.geoderma.2009.07.010. Beven, K., Kirkby, N., 1979. A physically based, variable contributing area model of basin hydrology. Hydrolog. Sci. Bull. 24 (1), 43–69. Bhunia, G.S., Shit, P.K., Pourghasemi, H.R., 2017. Soil organic carbon mapping using remote sensing techniques and multivariate regression model. Geocarto Int. 34 (2), 215–226. https://doi.org/10.1080/10106049.2017.1381179. Bleuler, M., Farina, R., Francaviglia, R., Di Bene, C., Napoli, R., Marchetti, A., 2017. Modelling the impacts of different carbon sources on the soil organic carbon stock and CO2 emissions in the Foggia province (Southern Italy). Agric. Syst. 157, 258–268. https://doi.org/10.1016/j.agsy.2017.07.017. Burgos, P., Madejón, E., Peréz-de-Mora, A., Cabrera, F., 2006. Spatial variability of the chemical characteristics of a trace-element-contaminated soil before and after remediation. Geoderma 130, 157–175. https://doi.org/10.1016/j.geoderma.2005.01. 016. Cambardella, C.A., Moorman, T.B., Novak, J., Parkin, T.B., Karlen, D.L., Turco, R., Konopka, A.E., 1994. Field-scale variability of soil properties in central lowa soils. Soil Sci Soc Am J. 58, 1501–1511. https://doi.org/10.2136/sssaj1994. 03615995005800050033x. Castrignanò, A., 2011. Introduction to Spatial Data Processing. Aracne, Rome, pp. 108. Chartin, C., Stevens, A., Kruger, I., Esther, G., Carnol, M., Wesemael, B.V., 2016. Estimating Soil Organic Carbon stocks and uncertainties for the National inventory Report – a study case in Southern Belgium. Geophys. Res. Abstr. 18 EGU2016-935-1. Colwell, J.E., 1974. Vegetation canopy reflectance. Remote Sens. Environ. 3, 175–183. Cooper, J., Burton, D., Daniell, T., Griffiths, B., Zebarth, J.B., 2011. Carbon mineralization kinetics and soil biological characteristics as influenced by manure addition in soil incubated at a range of temperatures. Eur. J. Soil Biol. 47, 392–399. https://doi.org/ 10.1016/j.ejsobi.2011.07.010. Costa, E., de Souza Tassinari, W., Pinheiro, H., Beutler, S., Anjos, L., 2018. Mapping soil organic carbon and organic matter fractions by geographically weighted regression. J. Environ. Qual. 47 (4), 718–725. https://doi.org/10.2134/jeq2017.04.0178. Deng, Y., Wu, C., Li, M., Chen, R., 2015. RNDSI: A ratio normalized difference soil index for remote sensing of urban/suburban environments. Int. J. Appl. Earth Observ. Geoinform. 39, 40–48. Durand, M.J.-H., Barbut, M.M., 1948. Carte des Sols d'Algérie. Constantine. Service Géographique de l'Armée (in French). Hartemink, A.E., McSweeney, K., 2014. Soil Carbon. Springer, Heidelberg. Hengl, T., 2007. A practical guide to geostatistical mapping of environmental variables. Geoderma 140, 417–427. Hengl, T., 2009. A Practical Guide to Geostatistical Mapping, second ed. University of Amsterdam, Amsterdam, pp. 270. Hengl, T., Heuvelink, G.B.M., Stein, A., 2004. A generic framework for spatial prediction of soil variables based on regression-kriging. Geoderma 120 (1–2), 75–93. Johnson, K.D., 2006. Testing a Geostatistical Regression Model as a Tool for the Management of Soil Organic Matter in a Tropical Watershed. ESE 502. https://www. seas.upenn.edu/~ese502/projects/Johnson_06.pdf (last access 29/07/2019). Keskin, H., Grunwald, S., 2018. Regression kriging as a workhorse in the digital soil mapper's toolbox. Geoderma 326, 22–41. https://doi.org/10.1016/j.geoderma.2018. 04.004. Kottek, M., Grieser, J., Beck, C., Rudolf, B., Rubel, F., 2006. World Map of Köppen-Geiger Climate Classification updated. Meteorol. Z. 15, 259–263. Kumar, S., 2013. Soil organic carbon mapping at field and regional scales using GIS and remote sensing applications. Adv. Crop Sci. Technol. 1 (2). https://doi.org/10.4172/ 2329-8863.1000e105. Kutsch, W.L., Bahn, M., Heinemeyer, A. (Eds.), 2009. Soil Carbon Dynamics: An Integrated Methodology. Cambridge University Press, Cambridge, UK (ISBN 978-0521-86561-6). Lal, R., 1999. World soils and the greenhouse effect. Global Change Newslett. 37, 4–5. Lark, R.M., 2000. A comparison of some robust estimators of the variogram for use in soil survey. Eur. J. Soil Sci. 51, 137–157. Li, J., Wen, Y., Li, X., Li, Y., Yang, X., Lin, Z., Song, Z., Mary Cooper, J., Zhao, B., 2018. Soil labile organic carbon fractions and soil organic carbon stocks as affected by longterm organic and mineral fertilization regimes in the North China Plain. Soil Tillage Res. 175, 281–290. https://doi.org/10.1016/j.still.2017.08.008. Lima, A., Cannon, A.W., Hsieh, W., 2016. Forecasting daily streamflow using online sequential extreme learning machines. J. Hydrol. 537, 431–443. https://doi.org/10. 1016/j.jhydrol.2016.03.017. Ling, L., Wang, H., Dai, W., Lei, X., Yang, X., Li, X., 2014. Spatial variability of soil organic carbon in the forestlands of northeast China. J. For. Res. 25, 867–876. https://doi. org/10.1007/s11676-014-0533-3. Long, X.-H., Zhao, J., Liu, Z.-P., Rengel, Z., Liu, L., Shao, H., Tao, Y., 2014. Applying geostatistics to determine the soil quality improvement by Jerusalem artichoke in coastal saline zone. Ecol. Eng. 70, 319–326. https://doi.org/10.1016/j.ecoleng.2014. 06.024. Ma, K., Liu, J., Zhang, Y., Parry, L., Holden, J., Ciais, P., 2015. Refining soil organic carbon stock estimates for China’s palustrine wetlands. Environ. Res. Lett. 10.

4. Conclusions The Tiffech watershed is mainly occupied by agricultural lands and has an important economic value for the country; this area is prone to SOC depletion due to intensive agricultural practices, in addition to increasing erosion phenomena. The determination of SOC stocks and the spatial distribution of SOC content in the Tiffech area is a tool to evaluate the soil fertility and the vulnerability to SOC loss. The geostatistical analysis revealed that SOC distribution across the Tiffech watershed was highly correlated with the topography of the area: highest values matched with the highest elevations, and the decrease in the middle of the watershed was explained as due to intensive agricultural practices and highly flooded soils. The application of OK and RK to interpolate SOC in unsampled locations provided distribution maps with an acceptable error, but an improvement in the interpolation was noticed using the RK method, i.e. a lower estimation error. The comparison of errors in OK and RK methods supports the idea that SOC interpolation based on measured data only could be effective in SOC mapping, however the introduction of auxiliary information derived from DEM and satellite images improved data interpolation in the Tiffech watershed. Being the accuracy of RK highly related to the spatial relationship between soil variables and environmental variables, and given the strong correlation of soil variables with some of the environmental variables available for the Tiffech watershed, the use of RK improved the accuracy of interpolation. Maps are to be meant as a tool to rate vulnerability to degradation of fertility and for decision making in soil management. The results of this study stress on the crucial role of land use in SOC stocks and soil fertility conservation. Moreover, the increase of the forest area or the combination of crops and forests could help in preventing SOC depletion. Based on SOC distribution maps, and the results of the correlation test, recommendations can be applied to control agricultural practices and to reduce the erosion phenomena in order to prevent depletion of SOC stocks and preserve soil fertility. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This work is part of the PhD thesis of Sana Boubehziz at Badji Mokhtar University – Annaba. The authors wish to thank Dr. Bruno Pennelli from CREA for revising the English language and style of the manuscript, the Editor and the anonymous reviewers for their valuable comments and suggestions. References Arpa Sardegna, 2018. http://www.sar.sardegna.it/servizi/agro/idrosuoli.asp (last access 29/07/2019). Arrouays, D., Minasny, B., McBratney, A., Grundy, M., McKenzie, N., Thompson, J., Gimona, A., Hong, S., Smith, S., Hartemink, A., Chen, S., Martin, M., Mulder, V.L.., Richer de Forges, A., Odeh, I.O.A., Padarian, J., Lelyk, G., Poggio, L., Savin, I., Hempel, J., 2017. Global soil map for soil organic carbon mapping and as a basis for global modeling. In: Global Symposium on Soil Organic Carbon, Rome, Italy, 21–23

10

Catena 190 (2020) 104539

S. Boubehziz, et al. https://doi.org/10.1088/1748-9326/10/12/124016. Marchetti, A., Piccini, C., Francaviglia, R., Santucci, S., Chiuchiarelli, I., 2010. Estimating soil organic matter content by regression kriging. In: Boettinger, J.L., Howell, D.W., Moore, A.C., Hartemink, A.E., Kienast-Brown, S. (Eds.), Digital Soil Mapping. Progress in Soil Science, vol. 2. Springer, Dordrecht, pp. 241–254. https://doi.org/10. 1007/978-90-481-8863-5_20. McBratney, A., Odeh, I.O.A., Bishop, T.S., Dunbar, M., Shatar, T., 2000. An overview of pedometric techniques for use in soil survey. Geoderma 97, 293–327. https://doi. org/10.1016/S0016-7061(00)00043-4. Mingjun, T., Lixiong, Z., Wenfa, X., Zhiling, H., Zhixiang, Z., Zhaogui, Y., Pengcheng, W., 2017. Spatial variability of soil organic carbon in Three Gorges Reservoir area, China. Sci. Total Environ. 599–600, 1308–1316. National Office of Geology, 1978. Geologic Map of Algeria, Mdaourouch Map. Ministry of Heavy Industry (in French). Nelson, D.W., Sommers, L.E., 1982. Total carbon, organic carbon, and organic matter. In: Page, A.L., Miller, R.H., Keeny, D.R., (Eds.), Methods of soil analysis, second ed. American society of agronomy, Madison, WI, pp. 539–579 (ASA and SSSA). Odeh, I.O.A., McBratney, A., Chittleborough, D.J., 1994. Spatial prediction of soil properties from landform attributes derived from a digital elevated model. Geoderma 63, 197–214. https://doi.org/10.1016/0016-7061(94)90063-9. Odeh, I.O.A., McBratney, A., Chittleborough, D.J., 1995. Further results on prediction of soil properties from terrain attributes: heterotopic cokriging and regression-kriging. Geoderma 67, 215–226. https://doi.org/10.1016/0016-7061(95)00007-B. Olaya, V., 2009. Basic Land-Surface Parameters. Geomorphometry—Concepts, Software, Applications. Developments in Soil Science. Elsevier B.V, Amsterdam, pp. 33. O'Rourke, S., Angers, D., Holden, N., McBratney, A., 2015. Soil organic carbon across scales. Glob. Change Biol. 21. https://doi.org/10.1111/gcb.12959. Pan, G.-X., Li, L.Q., Zhang, X.H., Dai, J.Y., Zhou, Y., Zhang, P.J., 2003. Soil organic carbon storage of China and the sequestration dynamics in agricultural lands. Adv. Earth Sci. 18, 609–618. Piccini, C., Marchetti, A., Francaviglia, R., 2014. Estimation of soil organic matter by geostatistical methods: Use of auxiliary information in agricultural and environmental assessment. Ecol. Ind. 36, 301–314. https://doi.org/10.1016/j.ecolind.2013. 08.009. QGIS Development Team, 2017. QGIS Geographic Information System. Open Source Geospatial Foundation Project. http://qgis.osgeo.org. R Development Core Team, 2016. R: A Language and Environment for Statistical Computing. R Development Core Team, Vienna, Austria (ISBN 3-900051-07-0). http://www.R-project.org (last access 29/07/2019). Roose, E.J., Lal, R., Feller, C., Barthès, B., Stewart, B.A., 2005. Soil Erosion and Carbon Dynamics (Advances in Soil Science). CRC Press, Boca Ratón, Florida, USA, pp. 376. SAGA User Group Association, 2017. SAGA 5.0.0 – System for Automated Geoscientific Analyses. http://www.saga-gis.org (last access 29/07/2019). Sherpa, S., Wolfe, D., van Es, H., 2016. Sampling and data analysis optimization for estimating soil organic carbon stocks in agroecosystems. Soil Sci. Soc. Am. J. 80, 1377–1392. https://doi.org/10.2136/sssaj2016.04.0113. Simbahan, C.G., Dobermann, A., Goovaerts, P., Ping, J., Haddix, M., 2006. Fine-resolution mapping of soil organic carbon based on multivariate secondary data. Geoderma 132, 471–489. https://doi.org/10.1016/j.geoderma.2005.07.001. Soleimani, A., Hosseini, S.M., Massah Bavani, A.R., Jafari, M., Francaviglia, R., 2017. Simulating soil organic carbon stock as affected by land cover change and climate

change, Hyrcanian forests (northern Iran). Sci. Total Environ. 599–600, 1646–1657. https://doi.org/10.1016/j.scitotenv.2017.05.077. Song, X.D., Liu, F., Ju, B., Zhi, J., Li, D., Zhao, Y., Zhang, G., 2017. Mapping soil organic carbon stocks of northeastern China using expert knowledge and GIS-based methods. Chin. Geogr. Sci. 27 (4), 516–528. https://doi.org/10.1007/s11769-017-0869-7. Sun, W., Minasny, B., McBratney, A., 2012. Analysis and prediction of soil properties using local regression-kriging. Geoderma 171–172, 16–23. https://doi.org/10.1016/ j.geoderma.2011.02.010. Teng, M., Lixiong, Z., Xiao, W., Huang, Z., Zhou, Z., Yan, Z., Wang, P., 2017. Spatial variability of soil organic carbon in Three Gorges Reservoir area, China. Sci. Total Environ. 599–600, 1308–1316. https://doi.org/10.1016/j.scitotenv.2017.05.085. Tondoh, J., Ouedraogo, I., Bayala, J., Tamene, L., Sila, A., Vågen, T.-G., Antoine, K., 2016. Soil organic carbon stocks in semi-arid West African drylands: implications for climate change adaptation and mitigation. Soil Discuss. 1–41. https://doi.org/10.5194/ soil-2016-45. United States Department of Agriculture’s Foreign Agricultural Service, 2018. Grain: World Markets and Trade. Wackernagel, H., 1994. Cokriging versus kriging in regionalized multivariate data analysis. Geoderma 62, 83–92. https://doi.org/10.1016/0016-7061(94)90029-9. Wang, Y., Deng, L., Wu, G., Wang, K., Shangguan, Z., 2017. Large-scale soil organic carbon mapping based on multivariate modelling: The case of grasslands on the Loess Plateau. Land Degrad. Dev. https://doi.org/10.1002/ldr.2833. Webster, R., 2001. Statistics to support soil research and their presentation. Eur. J. Soil Sci. 52, 331–340. https://doi.org/10.1046/j.1365-2389.2001.00383.x. Webster, R., Burgess, T.M., 1983. Spatial variation in soil and the role of kriging. Agric. Water Manage. 6, 111–122. https://doi.org/10.1016/0378-3774(83)90003-3. Webster, R., Oliver, M., 2007. Geostatistics for environmental scientists. In: Statistics in Practice, second ed. Wiley, Chichester, pp. 330. Xiao, J., Shen, Y., Tateishi, R., Bayaer, W., 2006. Development of topsoil grain size index for monitoring desertification in arid land using remote sensing. Int. J. Rem. Sens. 12 (27), 2411–2422. Yang, R.M., Zhang, G.L., Liu, F., Lu, Y.Y., Yang, F., Yang, M., Zhao, Y.G., Li, D.C., 2016. Comparison of boosted regression tree and random forest models for mapping topsoil organic carbon concentration in an alpine ecosystem. Ecol. Indicat. 60, 870–878. https://doi.org/10.1016/j.ecolind.2015.08.036. Yang, R.M., Zhang, G.L., Yang, F., Zhi, J.J., Yang, F., Liu, F., Zhao, Y.G., Li, D.C., 2016. Precise estimation of soil organic carbon stocks in the northeast Tibetan Plateau. Sci. Rep. 6, 21842. https://doi.org/10.1038/srep21842. Yorulmaz, A., Aydin, G., Atatanir, L., 2017. Determination of soil organic carbon levels using near infrared spectroscopy (nirs) in saline soils. Adü Ziraat Derg 14 (2), 29–31. https://doi.org/10.25308/aduziraat.310537. Zhang, S., Huang, Y., Shen, C., Ye, H., Du, Y., 2012. Spatial prediction of soil organic matter using terrain indices and categorical variables as auxiliary information. Geoderma 171–172, 35–43. https://doi.org/10.1016/j.geoderma.2011.07.012. Zhou, W., Liu, X.-Y., Yang, Ke, Bai, Z.-K., Cheng, H.-X., 2016. Research on spatial interpolation of soil organic carbon on reclaimed land at Pingshuo opencast coal mine. J. China Coal Soc. 41 (S1), 184–191. Zolfaghari, A.A., Taghizadeh-Mehrjardi, R., Moshki, A.R., Malone, B.P., Weldeyohannes, A.O., Sarmadian, F., Yazdani, M.R., 2016. Using the nonparametric k-nearest neighbor approach for predicting cation exchange capacity. Geoderma 265, 111–119.

11