Mapping Soil Organic Carbon Using Local Terrain Attributes: A Comparison of Different Polynomial Models

Accepted Manuscript

Title: Mapping Soil Organic Carbon Using Local Terrain Attributes: A Comparison of Different Polynomial Models

Author: SONG Xiaodong, LIU Feng, ZHANG Ganlin, LI Decheng, ZHAO Yuguo and YANG Jinling

PII: DOI: Reference:

S1002-0160(17)60445-4 10.1016/S1002-0160(17)60445-4 NA

To appear in:

Received date: Revised date: Accepted date:

NA NA NA

Please cite this article as: SONG Xiaodong, LIU Feng, ZHANG Ganlin, LI Decheng, ZHAO Yuguo and YANG Jinling, Mapping Soil Organic Carbon Using Local Terrain Attributes: A Comparison of Different Polynomial Models, Pedosphere (2017), 10.1016/S1002-0160(17)60445-4.

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT PEDOSPHERE Pedosphere ISSN 1002-0160/CN 32-1315/P

doi:10.1016/S1002-0160(17)60445-4

Mapping Soil Organic Carbon Using Local Terrain Attributes: A Comparison of Different Polynomial Models SONG Xiaodong, LIU Feng, ZHANG Ganlin*, LI Decheng, ZHAO Yuguo and YANG Jinling

cr ip

t

State Key Laboratory of Soil and Sustainable Agriculture, Institute of Soil Science, Chinese Academy of Sciences, Nanjing 210008 (China)

ABSTRACT:

Ac

ce pt

ed

M

an us

Derived directly from digital elevation model (DEM), local terrain attributes have been widely applied in digital soil mapping (DSM). This study aimed to evaluate the mapping accuracy of soil organic carbon (SOC) concentration by combining prediction methods with local terrain attributes calculated using different polynomial models. The prediction accuracy was used as a benchmark for those who may be more concerned with how accurately the variability of soil properties is modeled in practice, rather than how morphometric variables and their geomorphologic interpretations are understood and calculated. In this study, two neighborhood shapes (square and circular) and six representative algorithms (Evans-Young, Horn, Zevenbergen-Thorne, Shary, Shi and Florinsky algorithms) were applied. In general, 35 combinations of first- and second-order derivatives were produced as candidate predictors for soil mapping using two mapping methods, respectively (i.e. kriging with an external drift and geographically weighted regression). The results showed that appropriate local terrain attribute algorithms can better capture the spatial variation of SOC concentration in a region where soil properties are strongly influenced by the topography. When different combinations of first- and second-order derivatives were used, there was a best combination that leads to a more accurate estimate. For different prediction methods, the relative improvement in the two zones varied between 0.30% and 9.68%. The SOC maps resulting from the higher-order algorithms (Zevenbergen-Thorne and Florinsky) yielded less interpolation errors. Therefore, it is concluded that the performance of predictive methods which incorporate auxiliary variables could be improved by the attempts of different terrain analysis algorithms. Key Words: cross-validation, digital soil mapping, geographically weighted regression, kriging with an external drift, map accuracy

INTRODUCTION With respect to many chemical and physical processes, soil organic carbon (SOC) plays a major role in terrestrial ecosystems and often varies spatially due to the fast changing environmental conditions. A precise map of SOC concentration is often desirable for soil quality evaluation and soil management. With the rapid development of digital soil mapping in recent decades, various soil *

Corresponding author. E-mail: [email protected]. 2

ACCEPTED MANUSCRIPT

Ac

ce pt

ed

M

an us

cr ip

t

predictors have been recognized in prediction models among which topographic attributes have been most extensively applied (McBratney et al., 2003) as they have imposed a significant impact on soil properties. SOC abundance in topsoil is mainly dominated by net primary production of plants that enters soil upon senescence. Topography has been deemed as one of the most important indicators for the processes governing the lateral distribution of SOC, and hence digital terrain analysis has been used to capture the spatial variation of SOC (Gregorich et al., 1998). The SOC concentration of uncultivated soils mainly depends on the landform position because of soil erosion and leaching, which is predominant on sloping and undulating landscape. This is because relatively dry knolls and humid footslopes are formed by surface runoff and internal downslope flow of water through soils in saturated and unsaturated states. High soil moisture availability in these situations promotes plant growth and then affects the SOC concentration. Generally speaking, fundamental local terrain attributes have made diverse geographical effects on the SOC concentration (Behrens et al., 2010). For example, slope may affect the overland and subsurface flow velocity and runoff rate, and soil water content; aspect may control flow direction; and profile (or vertical) curvature may change the soil erosion and deposition rates (Wilson, 2012). A large number of fundamental topographic attributes have been proposed to quantitatively identify landform classes and features (Wilson, 2012). Many studies have shown that prediction methods incorporating these predictors outperform generic geostatistical models (e.g., ordinary kriging) (McBratney et al., 2003). The mapping performance would vary greatly if the methods incorporated into the same predictor derived from different algorithms or at different scales (Smith et al., 2006; Behrens et al., 2010; Pei et al., 2010; Cavazzi et al., 2013). Terrain attributes become selectively used and important particularly when several algorithms can be applied to generate the same attribute and the prediction performance may vary accordingly. Distinguished from non-local attributes, local terrain attributes are derived directly from digital elevation model (DEM) without additional inputs and are usually calculated by moving a three-by-three window (Shary et al., 2002; Behrens et al., 2010). For terrain attributes, the terms local and non-local are usually used regardless of the study scale or DEM resolution and associated with the mathematical sense of a particular variable (Florinsky, 2011). Some local terrain attributes such as slope, aspect and curvatures have been frequently adopted as predictors for estimating SOC (McBratney et al., 2003; Florinsky, 2011; Minár et al., 2013). Slope and aspect are also referred to as first-order derivatives, and twelve kinds of curvatures are referred to as second-order derivatives (Shary, 1995), as they are defined by the formulae depending on the firstand second-order partial derivatives of altitudes. Meanwhile, multifarious mathematical modeling algorithms have been developed to calculate those derivatives from a gridded DEM focusing on various landscapes (Evans, 1980; Horn, 1981; Shary, 1995; Shary et al., 2002; Shi et al., 2007; Minár et al., 2013). As the accuracy of the terrain attributes is unavoidably influenced by the quality of DEMs and calculation algorithms, numerous studies have been conducted and published focusing on estimating the accuracy of those algorithms (Schmidt et al., 2003; Warren et al., 2004), analyzing the relationships between errors of derived attributes and DEM data characteristics (Zhou and Liu, 2004; Wang et al., 2012a) and comparing calculated slope values with actual field measurements (Giles and Franklin, 1996; Warren et al., 2004; Shi et al., 2012). Studies also have shown the suitability of various polynomials under the specified conditions of complex landforms (Shary et al., 2002; Schmidt et al., 2003). Nevertheless, little research is concerning the context of soil mapping and results are not 3

ACCEPTED MANUSCRIPT

an us

cr ip

t

clearly applicable to knowledge-based digital soil mapping (Shi et al., 2012). After discussing “accuracy” of terrain attributes in detail, Shi et al. (2012) stated that the quality of the calculated values of terrain attributes is entirely applicable and situation-specific, and the evaluation has to be made based on what is used in specific applications. Yimer et al. (2006) addressed the effects of aspect and vegetation community types on soil physical and chemical properties. With ordination methods, Bodaghabadi et al. (2011) attempted to quantify relationships between the distribution of soil classes and terrain attributes. However, none of those studies have ever investigated the effect of local terrain attributes calculated by various mathematical modeling algorithms on predicting soil properties. The purpose of this research was to test and evaluate the role of local terrain attributes based on different polynomial models and their impacts on predicting SOC concentration by empirically using two typical prediction methods, namely, kriging with an external drift (KED) and geographically weighted regression (GWR). The local terrain attributes were derived from six polynomials and two types of neighborhood shapes (square and circular neighborhood). The ultimate objective of this application-based research was to evaluate two hypotheses: i) the overall mapping performance based on different polynomial models will be different; and ii) there is a best polynomial model for calculating first- and second-order local terrain attributes in terms of mapping accuracy. STUDY AREA AND DATA Study area

M

The study was carried out in the two zones (the upper and middle reach zones), which are located

ce pt

ed

along the northeast margin of the Qinghai-Tibetan Plateau in China at the intersection of the Tibetan Plateau, the Inner Mongolia-Xinjiang Plateau, and the Loess Plateau (Fig. 1). With the geographical boundary of about 97°20′ - 101°51′E and 37°41′ - 39°59′N, these two zones stretch for 340 km from the northwest to the southeast and covers an area of approximately 50,810 km2, accounting for about 36% of the total area of the Heihe River Basin. The landform declines gradually from the southwest to the northeast in the region, which ranges from 5,523 to 1,254 m in elevation (Fig. 1). Fig. 1 Location of the study area and distribution of soil sampling sites.

Ac

The main district of runoff formation in the whole Heihe River watershed is defined as the upper reach zone of the Heihe River (Qi and Luo, 2005). The district influenced by the Heihe River in Gansu Province is named as the middle reach zone (Fig. 1). The average slope of the upper reach zone and the middle reach zone is 17.6 and 2.4 degree respectively. The upper reach zone is primarily natural ecosystems and most areas are covered by alpine meadow; the middle reach zone is artificial oases composed primarily of irrigated farmland, and vast areas of Gobi and desert (Wang et al., 2012b). The upper reach zone is characterized by a humid and cold climate. The climate of the middle reach zone is a typical temperate arid environment with low precipitation and high potential evaporation. The distinct landscape and topography make these two zones ideal to study the effects of local terrain attributes in DSM separately. Soil data Part of the soil data (2010–2013), i.e. 251 soil points, were provided by the “Heihe Plan Science 4

ACCEPTED MANUSCRIPT

cr ip

t

Data Center, National Natural Science Foundation of China” (http://www.heihedata.org). In addition, we collected data from 404 sampling locations by purposive sampling with the method of Zhu et al. (2008). In this method, representative soil sites are selected from soil-scape units constructed by fuzzy c-means classification of pixels based on the soil forming factors (covariates). The locations of the sampling sites were located by a global positioning system (GPS) receiver. At each site, soil samples were collected at 0-20 cm. A total of 655 topsoil samples were compiled in a digital database. There were 379 and 276 sampling points in the upper and middle reach zones respectively. The samples were subsequently dried, sieved at 2 mm and analyzed using the Walkley–Black procedure for SOC (Walkley and Black, 1934). The Kolmogorov Smirnov (K-S) test (P = 0.000 < 0.05) rejected the null hypothesis of normality for samples. The SOC concentration data were transformed by natural logarithm to create an approximately normal distribution (referred to as LnSOC). Finally, the prediction values of SOC were back-transformed to original units. Environmental variables

Ac

ce pt

ed

M

an us

In this study, environmental variables utilized for the mapping of SOC concentration of 0-20 cm topsoil (g kg-1) were: first-order terrain derivatives (slope and aspect), second-order terrain derivatives (profile curvature, plan curvature, minimal curvature, maximal curvature and tangent curvature), elevation, topographic wetness index (TWI), mean annual air temperature (MAAT), mean annual precipitation (MAP), incoming solar radiation, normalized difference vegetation index (NDVI), land cover and soil type. Beijing-1 multispectral data were obtained from the Watershed Allied Telemetry Experimental Research (Li et al., 2009) at 32 m spatial resolution in 2008. The NDVI was computed in the ENVI software environment described as: NDVI = (NIR - red) / (NIR + red), where NIR is the near-infrared band and red is the red band. MAAT and MAP data were extracted from the database of National Meteorological Information Center, China Meteorological Administration (http://cdc.nmic.cn/home.do). Land cover and soil type data were obtained from the Data Sharing Infrastructure of Earth System Science (http://www.geodata.cn). The category variables used in a regression model should be replaced by dummy variables with the values of either 1 or 0, showing the presence or absence of each category. Based on the results obtained by analyzing variance (ANOVA) and Duncan’s method for multiple comparison (Duncan, 1955), the number of land cover categories was reduced to three: (1) cropland, village and barren, (2) forest, and (3) grassland and wetland. There were 14 soil types in this area. Eight primary soil types covered 89.5% of the study area, while the remaining 6 types of soil covered 10.5% of the study area. Using the aforementioned procedure, three soil type classes were combined: (1) Arenosols, Calcisols, Gypsisols, Phaeozems and Solonchaks, (2) Anthrosols, Fluvisols, Gleysols, Greyzems, Kastanozems, and (3) the others (FAO, IUSS Working Group WRB, 2007). SRTM (shorted for Shuttle Radar Topographical Mission) DEM data (Jarvis et al., 2008) were employed and geo-referenced from three arc second resolution to 30 m × 30 m resolution. The cosine function was selected to transform aspect data to the range from -1 to 1 indicating the north-facing degrees. Potential insolation (incoming solar radiation), expressed in kWatt m-2, mainly depends on elevation, slope and aspect, and thus can be derived directly from DEMs. Average annual potential insolation was calculated by means of the Solar Radiation function in ArcGIS 10.0. One year instead of a long-term average potential insolation was calculated, as the potential insolation only changes based on the scales of obliquity (Loutre et al., 2004). TWI was calculated based on the following 5

ACCEPTED MANUSCRIPT

cr ip

t

equation (Moore et al., 1993): TWI = ln(α / tanβ), where α is the accumulative upslope area per unit contour length computed with the D8 algorithm, and β is the local slope. TWI was extracted in System for Automated Geoscientific Analyses (SAGA) software, which will assign a more realistic, higher potential soil wetness than the TWID8 to grid cells situated in valley floors with a small vertical distance to a channel. The principal differences among most algorithms for the computing of local terrain attributes are the number of grid cells used and the weights given to each of those cell values. In general, most algorithms utilize some elevation values in a 3×3 window centered on the elevation cell in question and thus one can find all the unknown coefficients for a polynomial. However, a third-order polynomial should be fitted over all points in a 5×5 neighborhood for approximation of all the coefficients (Florinsky, 2009; Minár et al., 2013). Two types of the neighborhood (Shi et al., 2007), circular neighborhood and square neighborhood, and six mathematical algorithms were applied in this study (Table I).

Table I

an us

Polynomials and neighborhood types used in this study for the calculation of first- and second-order derivatives. Modified from Schmidt et al. (2003) Algorithm Reference Polynomial Algorithm Neighborhood Abbr. Calculation size and type

Evans-Young Evans (1980) Second (1978) Horn (1981)

Zevenbergen Zevenbergen

Unequal

3×3

weighting

Square

Partial

Closed solution

3×3

fourth

1st, 2nd

sHorn

1st

sZT

1st, 2nd

sShary

2nd

cShi,

1st

Square

ce pt

and Thorne

Second

ed

Horn

cEY,

derivatives

Circular/Square sEY

M

and Young

Least-squares fit 3×3

of

(1987)

Shary

Shary (1995)

Second

Least-squares fit, 3×3 match

Square

of central point

Shi et al.

Partial

(2007)

fourth

Florinsky

Third

Ac

Shi

Florinsky

Average solution 3×3

Circular/Square sShi Least-squares fit 5×5

(2009)

sFlo

1st, 2nd

Square

(1) The basis of the Evans-Young algorithm (Young, 1978; Evans, 1980) is to fit a second-order polynomial by least squares to a 3×3 square-gridded window of each sample point. The first- and second-order partial derivatives of elevation can be achieved by the finite difference formulae. This algorithm has been widely accepted and applied due to high accuracy (Schmidt et al., 2003) and argued to be the most optimal method under certain conditions (Shi et al., 2007). (2) As one of the most popular methods in practice, the Horn algorithm (1981) performs well for both slope and aspect estimation, and is employed in the “Slope” and “Aspect” function of the ArcInfo GRID software package and “r.slope.aspect” module of GRASS GIS. The Horn algorithm has 6

ACCEPTED MANUSCRIPT

Ac

ce pt

ed

M

an us

cr ip

t

already considered the effects of the orientation and assigns different weights to the points in the cardinal and diagonal directions. Therefore, the circular neighborhood was not implemented into this algorithm in this study. (3) The Zevenbergen and Thorne algorithm (1987) is more accurate than other algorithms under certain conditions (Jones, 1998; Shi et al., 2007). Incorporated in the ArcInfo “Curvature” command, this algorithm is also widely used, and was selected as a representative method based on the Lagrange polynomial. The circular neighborhood should not be considered, because this algorithm merely utilizes the elevation from four points in the cardinal directions and the result would be inconsistent with the shape of the neighborhood. (4) Both of Shary (Shary, 1995) and Evans-Young algorithms generate the same values of the first-order derivatives, slope and aspect. In the Shary algorithm, the central point should be matched within the quadratic trend surface to simulate the uniform gravity. Hence, the formulae for second partial derivatives are modified and a system of 12 curvatures was constructed (Shary, 1995). The Shary algorithm can generate different curvature values compared with Evans-Young algorithms, and this algorithm was chosen as a reference in calculating curvatures. (5) The Shi algorithm (2007) is an improved version of the Zevenbergen-Thorne algorithm. The novel feature is that the information from both cardinal and diagonal directions is taken into account and thus it is possible to evaluate the effect of the circular neighborhood with the Lagrange polynomial. The circular neighborhood algorithm has proven to be a good alternative to the traditional square neighborhood (Shi et al., 2007; 2012). (6) The Florinsky algorithm (Florinsky, 2009), as a representative algorithm of third-order polynomial, allows one to derive the first-, second-, and third-order partial derivatives within a 5×5 window. Compared with the Evans-Young algorithm, this algorithm allows one to derive more accurate models of topographic attributes. Local terrain attributes were selectively calculated using the aforementioned algorithms and the circular and square neighborhood (Table I), which finally resulted in a total of 35 combinations of first- and second-order derivatives (7 first-order derivatives, 5 second-order derivatives). Two predictor groups were grouped as follows: (1) Terrain attribute predictors (referred to as TA) include slope, aspect, profile curvature, plan curvature, minimal curvature, maximal curvature, tangent curvature, elevation, and TWI. (2) All predictors (referred to as ALL) contain terrain attribute predictors (TA), and other predictors (MAAT, MAP, solar radiation, NDVI, land cover and soil type). It is noted that each group contains 35 combinations of local terrain attributes. Using KED and GWR methods, two predictor groups resulted in the production of 280 maps for the upper and middle reach zones, and enabled a better understanding of the role of local terrain attributes within DSM. The formula for the aforementioned attributes can be found in literature (Horn, 1981; Zevenbergen and Thorne, 1987; Shary, 1995; Shary et al., 2002; Shi et al., 2007; Florinsky, 2011). Most of the local terrain attributes were generated using the Terrain Analysis function of ArcSIE® (http://www.arcsie.com/) and other algorithms were implemented in C++ using GDAL library. All predictors were resampled at a spatial resolution of 30 m. Step-wise regression was executed aiming to derive the best subset of predictor variables. The center scaling was performed upon all variables using the scale function in R (R Development Core Team, 2010) and thus the mean values and standard deviations of predictors were 0 and 1 respectively, aiming to compare the influence of variables on the spatial variation of SOC. 7

ACCEPTED MANUSCRIPT

METHODS Prediction methods Kriging with an external drift (KED) As one of the most frequently used geostatistical techniques, KED provides a means of integrating predictors into the interpolation process. The predictors are applied to describe the average shape of the dependent variable and the trend-free variogram is also required. Simple kriging is performed on the corresponding residuals after the deriving of the local mean of the dependent variable. With KED, the soil property of interest at a location i is predicted by: (1)

cr ip

t

n

yˆ KED  i   mˆ KED  i    KED  i   y  i   mˆ KED  i   1

an us

where yˆ (i) is the predicted soil property at location i, mˆ KED  i  the local mean estimated as linear combination of predictors and KED  i  the simple kriging weight attached to sampling location iα.

M

Geographically weighted regression (GWR) GWR is an extended version of ordinary least squares regression considering the spatial nonstationarity in predictor relationships and utilizing spatially varying coefficients in local linear models (Fotheringham et al., 2002). In GWR, a local regression model is fitted at each prediction location, and the soil property of interest is predicted by (Fotheringham et al., 2002): K

ed

yˆ  i   ˆ0  i    ˆk  i  xk  i 

(2)

k 1

ce pt

where ˆ0 i  is the estimated intercept at location i, ˆk i  the estimated regression coefficient for predictor k and xk(i) the value for the kth predictor at location i. The location-specific regression coefficients ˆ0 i  and ˆk  i  (k = 1...K) are estimated by weighted least squares, with weights related

Ac

to the distance between the prediction location and the sampling locations. The closer a sampling location is to prediction location i, the larger its weight. The soil property of interest at a location i is predicted as follows: (i) draw a circle of given bandwidth around prediction location i, (ii) weight a measurement of the soil property of interest at a sampling location in accordance with its proximity to prediction point i and (iii) estimate the regression coefficients by: 1 βˆ  i    XT W  i  X  XT W  i  y

(3)

where T denotes the transposition of a matrix, X the design matrix formed by the values of predictors xk(i), W(i) a geographical weight matrix for the prediction location i, and y the vector with measurements of the soil property of interest at the sampling locations. The selection of kernel is especially important, as the observations in the estimation process will be determined by the spatial kernel. Two kinds of spatial kernel, namely, fixed kernel and adaptive kernel, are frequently used. The adaptive shape is defined by the number of the nearest neighbors, whose size will increase rapidly to capture enough observations when they are sparsely distributed, and vice versa. Furthermore, the 8

ACCEPTED MANUSCRIPT adaptive spatial kernel was employed for this study in which sample density varied spatially. All the predictive approaches and statistical analysis in this study were performed using the statistical software R (R Development Core Team, 2010) with its add-in packages: “gstat” (Pebesma and Graeler, 2015) and “GWmodel” (Lu et al., 2014). Evaluation

an us

cr ip

t

Leave-One-Out cross-validation (LOOCV), which is a special case of cross-validation where the number of folds equals the number of instances in the dataset, was conducted to evaluate the accuracy of different models through three statistical measurements of the prediction errors. The accuracy of estimates was assessed by the mean absolute error (MAE) and the root mean squared error (RMSE). The values of MAE and RMSE should be as low as possible for accurate interpolation. The relative improvement (RI) (%) for different combinations of local terrain attributes was employed to compare their performance in special cases, where the same predictive method was applied: RI = 100 × (RMSEw - RMSEb) / RMSEw , where RMSEw and RMSEb are the root mean square errors of the worst and best estimation within one case, respectively. The variance inflation factors (VIF) were calculated for the GWR model and thus VIF values for each sampling site can be summarized to diagnose the collinearity problem (Wheeler and Tiefelsdorf, 2005). The global regression, multiple linear regression (MLR) was also carried out for comparison. The VIF is computed as: 1 1  Rk2

M

VIFk 

2

(4)

ed

where Rk denotes the coefficient of determination of predictor k. Generally speaking, a VIF > 10

ce pt

implies a collinearity problem. RESULTS

Spatial variability

Ac

The semivariograms of LnSOC in the two zones are illustrated in Fig. 2. The best-fit variogram model for the variograms was Exponential. The auto-correlation range in the upper reach zone was 0.79 km demonstrating the spatial auto-correlation on a small scale. The nugget/sill ratio of SOC in this zone was 11.76%, which indicates that 11.76% of the variability of SOC consisted of unexplainable or random variations. This nugget/sill ratio reflected that SOC in the upper reach zone had strong spatial dependence. The SOC in the middle reach zone had a larger auto-correlation range (16.12 km) with a nugget/sill ratio of 70.27%. Fig. 2 Semivariogram of natural log-transformed SOC (LnSOC) in the upper and middle reach zones. Evaluation of variable importance The analysis of the local regression coefficients within the GWR model was carried out, which 9

ACCEPTED MANUSCRIPT

cr ip

t

aims to measure the variable importance in the local area where relationships within the variables are usually spatially varying. All variables were standardized so as to measure their importance in the local regression. The higher the absolute value of the mean regression coefficient is, the more important this predictor will be. Based on the ALL predictors, the spatial distribution of coefficients of the GWR model across the two zones is presented in Fig. 3. The GWR coefficients varied spatially, which illustrated that the weights of each variable and the cross correlations between LnSOC and variables were changing at different sites. In the upper reach, the local regression coefficients of ALL predictors were continuous and varied slightly except Soil2 and MAAT. Local regression coefficients of local terrain attributes were with smaller magnitude than other predictors. The results suggested that the local terrain attributes have exerted similar roles in the different prediction methods for the same study area.

an us

Fig. 3 The regression coefficients of geographical weighted regression based on ALL predictors in the upper reach zone (left) and middle reach zone (right). Soil = soil type; Land = land cover; MinCur = minimal curvature; ProCur = profile curvature; MaxCur = maximal curvature; NDVI = normalized difference vegetation index; TWI = topographic wetness index; MAP= mean annual precipitation; MAAT = mean annual air temperature; Solar = incoming solar radiation. Collinearity in local regression

Table II

ce pt

ed

M

Collinearity of predictors may lead to unrealistic estimates of regression coefficients (Wheeler and Tiefelsdorf, 2005), and as a consequence possibly unrealistic predictions with large errors. The VIF values for each predictor in the local regression (GWR) are summarized in Tables II and III respectively. The global regression, multiple linear regression (MLR) was also carried out for comparison. For the upper reach, the mean VIF values of GWR and MLR were smaller than 10 indicating that the collinearity problem was not severe. On the other hand, mean VIF values of MAP and MAAT were much large than 10 showing a somewhat stronger collinearity problem in GWR in the middle reach. This issue could be explained by the moderate spatial variation of MAP and MAAT in this area.

a)

Ac

Variance inflation factors for the MLR and GWR model in the upper reach Itemb) Soil3 Soil2 Land3 Land2 MinCur Slope NDVI TWI MAP MAAT Solar Elevation 2.42 2.42

1.50

1.50 1.16

3.34

1.98

2.84

2.30

2.25

1.65

3.14

GWR Minimum 1.98 1.46

1.39

1.09 1.12

2.76

1.80

2.55

2.00

1.96

1.45

2.21

25%

2.07 1.50

1.39

1.14 1.14

2.87

2.01

2.74

2.34

2.20

1.45

2.30

Median

2.19 1.58

1.45

1.17 1.15

3.03

2.08

2.91

2.59

2.60

1.48

2.48

Mean

2.23 1.63

2.86

2.63 1.20

3.24

2.20

3.07

3.13

3.36

1.82

3.24

75%

2.23 1.77

1.59

1.19 1.17

3.55

2.35

3.46

3.56

3.86

1.56

2.66

Maximum 3.25 2.32 66.95

67.61 2.17

4.99

2.84

4.17

18.23 22.15

5.24

13.64

MLR

a)

See Fig. 3 for descriptions of the abbreviations of predictors.

b)

MLR = multiple linear regression; GWR = geographically weighted regression. 10

ACCEPTED MANUSCRIPT

Table III a)

Variance inflation factors for the MLR and GWR model in the middle reach Itemb) Soil3 Soil2 Land3 Land2 ProCurSlope NDVI MaxCur MAP

MAAT Elevation

MLR 1.79

1.79

2.12

3.51 1.35

3.16

8.81

6.80

5.85

GWR Minimum 1.02 1.50

5.16

5.07

1.53

2.10 1.28

1.68

3.77

4.45

2.48

25%

1.04 1.54

5.49

5.25

1.76

3.26 1.30

2.13

10.43

7.09

5.47

Median

1.06 1.56

5.64

5.28

1.91

3.85 1.31

2.62

18.23

11.03

7.28

Mean

1.10 1.55

8.56

8.18

1.98

3.69 1.38

2.66

17.60

11.74

6.92

75%

1.09 1.56

6.24

5.79

2.00

4.00 1.34

2.81

207.65 207.88

5.09

4.91 2.57

6.22

cr ip

Maximum 1.85 1.78

t

1.58 1.58

19.02

11.66

7.45

60.65

49.32

14.85

See Fig. 3 for descriptions of the abbreviations of predictors.

b)

MLR = multiple linear regression; GWR = geographically weighted regression.

an us

a)

Prediction accuracy

ce pt

ed

M

Results of the validation indices are shown in Fig. 4, where each box is one case including 35 local terrain attributes. There were only slight outliers and long box lengths in these two zones. The box lengths of these two methods were shorter in the upper reach for KED method. KED showed more stable performance for different groups of local terrain attributes both in the upper and middle reach zones. It was suggested that KED accurately modeled the spatial patterns of SOC, in which KED achieved better performance by fitting the randomness (unexplained variability) of global regression residuals, particularly when adequate predictors were applied. Broadly speaking, prediction results with the GWR model improved when more predictors were incorporated in the middle reach. This is because the adaptive spatial kernel was employed and hence the prediction performance was affected by the varied bandwidth optimized by cross-validation.

Ac

Fig. 4 Box plots of SOC concentration prediction accuracy of different cases. There are 35 results in each box, where local terrain attributes were calculated by different algorithms. ALL = all predictors; TA = terrain attribute predictors. The best and worst cross validation results within 280 times prediction are given in Tables IV and V. Most results within one case were moderately different, except KED in the upper reach zone. The relative improvement (RI) values within different groups of predictors revealed that there were moderate differences among these results. Based on the above prediction accuracy analysis, we can accept the first hypothesis, which assumes that the overall mapping performance based on different polynomial models will be different. The RI values of GWR(ALL) in the two zones showed that the local terrain attributes played more important roles in those two methods than KED. However, the groups of local terrain attribute algorithms in each best case were different. This scenario can be explained by the different predictive assumptions of the KED and GWR, which is whether the model residuals are assumed to be dependent and the regression model is fitted at each prediction location. 11

ACCEPTED MANUSCRIPT Furthermore, even if there is the best algorithm for calculating local terrain attributes, the mapping quality may also be affected by the uncertainties resulting from the prediction procedure, soil sampling and DEM data.

Table IV The best and worst cross validation results in the upper reach. Group means the combination of first- (left abbreviation) and second-order (right abbreviation) derivative algorithms, where abbreviations of algorithms were given in Table I. Best case Worst case Method

Group

MAEa) RMSEb)

Group

MAE

RMSE

RIc)

g kg-1 g kg-1

g kg-1

g kg-1

g kg-1

%

sFlo.sZT

0.489 0.662

sEvans.cEvans

0.490

0.664

0.301

KED(TA)

sZT.sFlo

0.596 0.808

sEvans.cEvans

0.599

0.814

0.737

GWR(ALL)

sHorn.sEvans

0.533 0.699

cEvans.sShary

0.552

0.738

5.285

GWR(TA)

sZT.sEvans

0.579 0.769

sEvans.sZT

0.609

0.791

2.781

b)

RMSE = root mean squared error

c)

RI = relative improvement.

cr ip

MAE = mean absolute error.

an us

a)

t

KED(ALL)

Table V

Method

Group

MAEa) RMSEb)

KED(TA)

sZT.sEvans

GWR(ALL)

sFlo.sFlo

GWR(TA) a)

sFlo.sZT

MAE

RMSE RIc)

g kg-1

g kg-1

%

0.279

0.415

sShi.sFlo

0.291

0.435

4.598

0.265

0.401

sFlo.sShary

0.269

0.418

4.067

0.283

0.420

sZT.cEvans

0.313

0.465

9.677

sZT.cEvans

0.324

0.489

4.090

ce pt

sFlo.sEvans

Group

g kg-1

ed

g kg-1 KED(ALL)

M

The best and worst cross validation results in the middle reach Best case Worst case

0.312

0.469

MAE = mean absolute error.

RMSE = root mean squared error

c)

RI = relative improvement.

Ac

b)

The sZT and sFlo algorithms performed best for calculating first-order derivatives, and sEvans, sZT and sFlo performed best for the second-order derivatives in terms of RMSE. The successful application of sZT and sFlo algorithms (high order polynomial) concurred with previous studies made by Zevenbergen and Thorne (1987) and Schmidt et al. (2003) that these algorithms can lead to a local denoising and hence greatly enhance the calculation of partial derivatives, even though SRTM DEM data were resampled from 90m spatial resolution with the ±15 m vertical accuracy. The circular neighborhood method did not outperform the square neighborhood methods, which confirmed that this method is more advantageous on a high-resolution DEM (Shi et al., 2007). It may also imply that it is meaningful to discuss the importance of terrain variables in high-relief areas. Therefore, the second hypothesis, concerning whether there is one best polynomial model for calculating first- and second-order local terrain attributes can be accepted under certain conditions. As the ALL predictors achieved best performance, the best cases of different prediction methods 12

ACCEPTED MANUSCRIPT in the two zones were compared independently (Fig. 5). Using ALL predictors, KED outperformed GWR for the two zones. The predictions were back-transformed on the log-scale by exponentiation. Figs. 6 and 7 show the spatial predictions for SOC based on ALL predictors. The maps produced by the two methods reflected much spatial variation, which represents the characteristics of the environmental variables. Both KED and GWR generated similar results in the northwestern part of the middle reach, which was caused by the low variation and strong spatial autocorrelation of SOC concentrations in this area covered by desert and oasis supplied by the Heihe River.

t

Fig. 5 Comparison of predicted/observed LnSOC concentration in the topsoil calculated from best and worst case based on ALL predictors in the upper and middle reach zones.

cr ip

Fig. 6 Predicted soil organic carbon (SOC) maps based on KED(ALL) (a) and GWR(ALL) (b) in the upper reach zone.

an us

Fig. 7 Predicted SOC maps by KED(ALL) (a) and GWR(ALL) (b) in the middle reach zone. DISCUSSION

Ac

ce pt

ed

M

Generally, different combinations of local terrain attributes showed different prediction accuracies in terms of leave-one-out cross validation. Pei et al. (2010) concluded that the accuracy of soil organic carbon maps can be enhanced when better algorithms were adopted to derive topographic wetness index. They attained the results with RI of 7.85% and 2.36% for the model of simple kriging with varying local means and collocated cokriging respectively. Concurring with those results, the RI values achieved in this study ranged widely from 0.30% to 9.68%, proving that the mapping performance can be improved by attempting different digital terrain analysis algorithms. Based on the analysis of the geometrical properties of the land surface, local terrain attributes derived from different algorithms can be applied to enhance the prediction accuracy, even if the two zones have different topography. An interesting finding in the cross-validation of the upper reach was that significant correlations between local terrain attributes and SOC did not obtain high RI values, even if the average slope of the upper reach was relatively larger than that of the middle reach (17.6 vs. 2.4 degree). The large variation of RI values in the middle reach may result from the high uncertainty in this area. It is also worth to note that the “best” terrain calculation algorithm combination is conditional in practice, which depends on what other predictors are included. This is because local terrain attributes played less important roles than those of other predictors, and hence local terrain attributes cannot fully capture the spatial variation of SOC concentration. Furthermore, different groups of predictors might slightly affect the mapping performance. In addition, the uncertainty also will be involved by the prediction techniques, which can model different spatial heterogeneity feature of SOC concentration. The resampling of SRTM DEM data with the 90m spatial resolution may involve certain errors for calculating local terrain attributes. For this scale, Johnson et al. (2009) concluded that the terrain attributes can capture larger scale influences of soil transport and not those occurring at a specific point (i.e. the quantitative pit). Morphologically varied areas may show good performance in soil mapping at coarser resolutions and large window sizes, as well as fine resolutions (i.e. 20 m) with small windows (Smith et al., 2006; Erskine et al., 2007; Behrens et al., 2010; Cavazzi et al., 2013). 13

ACCEPTED MANUSCRIPT

Ac

ce pt

ed

M

an us

cr ip

t

Related studies have also suggested that there is a significant topographic control on SOC in the topsoil in humid and semi-arid regions (Schwanghart and Jarmer, 2011). There are mainly two approaches for topography to affect SOC concentration in the topsoil by adding lateral components (Schwanghart and Jarmer, 2011). On the one hand, the water balance is directly influenced by topography and hence the lateral water exchange is facilitated by surface runoff. This process results in various moisture conditions, which influences the growth of vegetation in our study area. On the other hand, soil material separated by soil erosion is transported by surface runoff and laterally relocated in the landscape (Martínez-Mena et al., 2008). At the hillslope scale, the SOC in the topsoil is significantly correlated with local terrain attributes, as surface runoff may play an important role in accumulating SOC and runon may increase local water availability and favor plant growth. How well the land surface representation based on terrain analysis algorithms matches the “real terrain” may greatly affect the quality of the resulting soil maps. Here “real terrain” is controlled by the scientific background, in that the earth’s surface is specific to the application and therefore terrain provides the benchmark for error measurement (Schneider, 2001). Consequently, the aforementioned processes will be indirectly affected by different algorithms applied for calculating local terrain attributes. The prediction accuracy will be enhanced by the local terrain attributes based on polynomial models that better describe the patterns of landforms. One of the overall aims of this study was to evaluate the suitability of terrain analysis algorithms by comparing the accuracies of SOC maps. Different from quantitative surface analysis (Jones, 1998; Zhou and Liu, 2004), we used an application-specific scenario for the mapping of SOC at regional scale. Different combinations of first- and second-order derivatives provided diverse SOC maps due to their describing abilities of the general geomorphometry of the land surface. In this study, the sFlo algorithm achieved a large variation in prediction accuracy, and the sEvans, sZT and sFlo algorithms performed better. This result is similar to that of the comparative studies of polynomial models that modeling results from higher-order algorithms showed higher sensitivity to local variations (Schmidt et al., 2003; Florinsky, 2009). It was also beneficial to arrive at a conclusion that we could achieve an optimal combination of first- and second-order derivatives based on disparate algorithms rather than the same algorithm. In common with one of the objectives of geomorphometry, to a certain extent, digital soil mapping aims to quantitatively describe and model the variation of soil properties in terrestrial ecosystems. Therefore, how to utilize the best algorithm for terrain attribute calculation should be paid much more attention. CONCLUSIONS 1.

Appropriate terrain attribute algorithms can be applied to improve prediction performance. When the first- and second-order derivatives were combined with different predicting methods, there was a best combination that can lead to a more accurate prediction map of SOC. 2. The maps based on higher-order approaches yielded moderately smaller prediction errors, indicating that local denoising enhanced the soil mapping in this study. 3. The results of the current study can be used as a benchmark of the suitability of local terrain attributes based on different algorithms. Nevertheless, it does not mean that the best algorithm for the calculation of local terrain attributes in this study will outperform others with different spatial resolutions and neighborhood sizes, especially when the DEM data are generated differently due to the accuracy of the DEM. 14

ACCEPTED MANUSCRIPT

ACKNOWLEDGMENT The study was supported financially by the National Natural Science Foundation of China (Nos. 41130530, 91325301, 41401237, 41371224, and 41571212) and partly by the Jiangsu Province Science Foundation for Youths (No. BK20141053). The authors are grateful to Prof. Ian S. Evans, University of Durham, for sharing his published materials and to Dr. Florinsky, Institute of Mathematical Problems of Biology, Russian Academy of Sciences, for the interpretation of his method.

t

REFERENCES

Ac

ce pt

ed

M

an us

cr ip

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. Bodaghabadi M B, Salehi M H, Martínez-Casasnovas J A, Mohammadi J, Toomanian N, Borujeni I E. 2011. Using Canonical Correspondence Analysis (CCA) to identify the most important DEM attributes for digital soil mapping applications. Catena. 86: 66--74. Cavazzi S, Corstanje R, Mayr T, Hannam J, Fealy R. 2013. Are fine resolution digital elevation models always the best choice in digital soil mapping? Geoderma. 195-196: 111--121. Duncan D B. 1955. Multiple range and multiple F tests. Biometrics. 11: 1--42. Erskine R H, Green T R, Ramirez J A, MacDonald L H. 2007. Digital elevation accuracy and grid cell size: effects on estimated terrain attributes. SOIL SCI SOC AM J. 71: 1371--1380. Evans I S. 1980. An integrated system of terrain analysis and slope mapping. Zeitschrift für Geomorphologie, N.F., Supplementband. 36: 274--295. FAO, IUSS Working Group WRB, 2007. World Reference Base for Soil Resources. ISRIC, Rome. Florinsky I V. 2009. Computation of the third-order partial derivatives from a digital elevation model. INT J GEOGR INF SCI. 23: 213--231. Florinsky I V. 2011. Digital terrain analysis in soil science and geology. Elsevier Academic Press, Amsterdam, pp. 7--16. Fotheringham A S, Brunsdon C, Charlton M. 2002. Geographically weighted regression: The analysis of spatially varying relationships. Chichester: Wiley. Giles P T, Franklin S E. 1996. Comparison of derivative topographic surfaces of a DEM generated from stereoscopic SPOT images with field measurements. PHOTOGRAMM ENG REM S. 62: 1165--1171. Gregorich E G, Greer K J, Anderson D W, Liang B C. 1998. Carbon distribution and losses: erosion and deposition effects. SOIL TILL RES. 47: 291--302. Horn B K P. 1981. Hill shading and the reflectance map. Proceeding of the IEEE. 69: 14--47. Jarvis A, Reuter H I, Nelson A, Guevara E. 2008. Hole-filled SRTM for globe version 4. Available from CGIAR-CSI SRTM 90 m Database http://srtm.csi.cgiar.org. Johnson K D, Scatena F N, Johnson A H, Pan Y. 2009. Controls on soil organic matter content within a northern hardwood forest. Geoderma. 148: 346--356. Jones K H. 1998. A comparison of algorithms used to compute hill slope as a property of the DEM. COMPUT GEOSCI-UK. 24: 315--323. Li X, Li X, Li Z, Ma M, Wang J, Xiao Q, Liu Q, Che T, Chen E, Yan G, Hu Z, Zhang L, Chu R, Su P, 15

ACCEPTED MANUSCRIPT

Ac

ce pt

ed

M

an us

cr ip

t

Liu Q, Liu S, Wang J, Niu Z, Chen Y, Jin R, Wang W, Ran Y, Xin X, Ren H. 2009. Watershed Allied Telemetry Experimental Research. J GEOPHYS RES. 114: 1--19. Loutre M F, Paillard D, Vimeux F, Cortijo E. 2004. Does mean annual insolation have the potential to change the climate? EARTH PLANET SC LETT. 221: 1--14. Lu B, Harris P, Charlton M, Brunsdon C, Nakaya T, Gollini I. 2014. GWmodel: Geographically weighted models. R package version 1.2-3. Martínez-Mena M, Lopez J, Almagro M, Boix-Fayos C, Albaladejo J. 2008. Effect of water erosion and cultivation on the soil carbon stock in a semiarid area of southeast Spain. SOIL TILL RES. 99: 119--129. McBratney A B, Mendonça Santos M L, Minasny B. 2003. On digital soil mapping. Geoderma. 117: 3--52. Minár J, Jenčo M, Evans I S, Minár J, Kadlec M, Krcho J, Pacina J, Burian L, Benová A. 2013. Third-order geomorphometric variables (derivatives): definition, computation and utilization of changes of curvatures. INT J GEOGR INF SCI. 27: 1381--1402. Moore I D, Gessler P E, Nielsen G A, Peterson G A. 1993. Soil attribute prediction using terrain analysis. SOIL SCI SOC AM J. 57: 443--452. Pebesma E, Graeler B. 2015. gstat: Spatial and Spatio-Temporal Geostatistical Modelling, Prediction and Simulation. R package version 1.1-0. Pei T, Qin C Z, Zhu A X, Yang L, Luo M, Li B L, Zhou C H. 2010. Mapping soil organic matter using the topographic wetness index: A comparative study based on different flow-direction algorithms and kriging methods. ECOL INDIC. 10: 610--619. Qi S Z, Luo F. 2005. Water environmental degradation of the Heihe River basin in arid northwestern China. ENVIRON MONIT ASSESS. 108: 205--215. R Development Core Team. 2010. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, available at: http://www.R-project.org (accessed 10/15/2014). Schmidt J, Evans I S, Brinkmann J. 2003. Comparison of polynomial models for land surface curvature calculation. INT J GEOGR INF SCI. 17: 797--814 Schneider B. 2001. Phenomenon-based specification of the digital representation of terrain surfaces. Transactions in GIS. 5: 39--52. Schwanghart W, Jarmer T. 2011. Linking spatial patterns of soil organic carbon to topography - A case study from south-eastern Spain. Geomorphology. 126: 252--263. Shary P A. 1995. Land surface in gravity points classification by complete system of curvatures. MATH GEOL. 27: 373--390. Shary P A, Sharaya L S, Mitusov A V. 2002. Fundamental quantitative methods of land surface analysis. Geoderma. 107: 1--32. Shi X, Zhu A X, Burt J, Choi W, Wang R X, Pei T, Li B L, Qin C Z. 2007. An experiment using a circular neighborhood to calculate slope gradient from a DEM. PHOTOGRAMM ENG REM S. 73: 143--157. Shi X, Girod L, Long R, DeKett R, Philippe J, Burke T. 2012. A comparison of LiDAR-based DEMs and USGS-sourced DEMs in terrain analysis for knowledge-based digital soil mapping. Geoderma. 170: 217--226. Smith M P, Zhu A X, Burt J E, Stiles C. 2006. The effects of DEM resolution and neighborhood size on digital soil survey. Geoderma. 137: 58--69. 16

ACCEPTED MANUSCRIPT

Ac

ce pt

ed

M

an us

cr ip

t

Walkley A, Black I A. 1934. An examination of the Degtjareff method for determining soil organic matter and a proposed modification of the chromic acid titration method. SOIL SCI. 37: 29--38. Wang C, Yang Q, Guo W, Liu H, Jupp D, Li R, Zhang H. 2012a. Influence of resolution on slope in areas with different topographic characteristics. COMPUT GEOSCI-UK. 41: 156--168. Wang X, Ma M, Huang G, Veroustraete F, Zhang Z, Song Y, Tan J. 2012b. Vegetation primary production estimation at maize and alpine meadow over the Heihe River Basin, China. INT J APPL EARTH OBS. 17: 94--101. Warren S D, Hohmann M G, Auerswald K, Mitasova H. 2004. An evaluation of methods to determine slope using digital elevation data. Catena. 58: 215--233. Wheeler D C, Tiefelsdorf M. 2005. Multicollinearity and correlation among local regression coefficients in geographically weighted regression. Journal of Geographical Systems, 7: 161--187. Wilson J P. 2012. Digital terrain modeling. Geomorphology. 137: 107--121. Yimer F, Ledin S, Abdelkadir A. 2006. Soil property variations in relation to topographic aspect and vegetation community in the south-eastern highlands of Ethiopia. FOREST ECOL MANAG. 232: 90--99. Young M. 1978. Terrain analysis: program documentation. Report 5 on Grant DA-ERO-591-73-G0040, ‘Statistical characterization of altitude matrices by computer’. Department of Geography, University of Durham, England. pp. 27. Zevenbergen L W, Thorne C R. 1987. Quantitative analysis of land surface topography. EARTH SURF PROC LAND. 12: 47--56. Zhou Q, Liu X. 2004. Analysis of errors of derived slope and aspect related to DEM data properties. COMPUT GEOSCI-UK. 30: 369--378. Zhu A X, Yang L, Li B, Qin C, English E, Burt J E, Zhou C H. 2008. Purposive sampling for digital soil mapping for areas with limited data. In: Hartemink, A E, McBratney, A B, Mendonca Santos, ML (Eds.), Digital Soil Mapping with Limited Data. Springer-Verlag, New York, pp. 233--245.

17

ACCEPTED MANUSCRIPT

an us

cr ip

t

FIGURES

Fig. 2

Ac

ce pt

ed

M

Fig. 1

18

Ac

ce pt

ed

M

an us

Fig. 3

cr ip

t

ACCEPTED MANUSCRIPT

Fig. 4

19

an us

cr ip

t

ACCEPTED MANUSCRIPT

Fig. 6

Ac

ce pt

ed

M

Fig. 5

Fig. 7

20