Mapping key soil properties to support agricultural production in Eastern China

Accepted Manuscript Modelling and mapping of key soil properties to support agricultural production in Eastern China

Yuxin Ma, Budiman Minasny, Chunfa Wu PII: DOI: Reference:

S2352-0094(17)30061-5 doi: 10.1016/j.geodrs.2017.06.002 GEODRS 128

To appear in:

Geoderma Regional

Received date: Revised date: Accepted date:

24 March 2017 1 June 2017 13 June 2017

Please cite this article as: Yuxin Ma, Budiman Minasny, Chunfa Wu , Modelling and mapping of key soil properties to support agricultural production in Eastern China, Geoderma Regional (2017), doi: 10.1016/j.geodrs.2017.06.002

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 Modelling and mapping of key soil properties to support agricultural production in Eastern China Yuxin Maa*, Budiman Minasnya, Chunfa Wua, b a

Faculty of Agriculture & Environment, The University of Sydney, NSW 2006, Australia

b

Department of Agricultural Resources and Environment, Nanjing University of Information Science

SC

RI

PT

and Technology, Nanjing 210044, China

NU

*Corresponding author: Tel.: +61 0432036507.

AC

CE

PT E

D

MA

E-mail: [email protected]

ACCEPTED MANUSCRIPT Abstract The spatial prediction of soil attributes plays important roles in establishing areas for agricultural development and sustainable land management. This paper aims to model and map soil properties (texture, pH, and SOM content) across a 210,000 km2 region in Eastern China using 2895 topsoil legacy soil observations coupled with freely available covariates (DEM, Landsat, Climate and Radar data) at a

PT

spatial resolution of 90m. A decision tree model (Cubist) and Regression Kriging approach were used

RI

in the modelling. The prediction results showed that local Regression Kriging produces high accuracy in predicting the spatial distribution of soil attributes that are affected by human activities (soil pH

SC

and organic matter content), where the accuracy and precision of local RK appears superior to

NU

the Cubist method. Moreover, DEM at different resolution and climate factors were found to be important predictors in this study area. This study can be readily expanded to the whole

MA

China using existing soil data at a relatively high resolution (90 m) to contribute to the GlobalSoilMap project.

AC

CE

PT E

D

Keywords: Soil properties; Spatial distribution; Environmental factors; Cubist; RK; Uncertainty, Aridisols

ACCEPTED MANUSCRIPT 1. Introduction Soil is the centre of terrestrial ecosystems and interfaces with lithosphere, hydrosphere, atmosphere and biosphere. The knowledge on the spatial prediction of soil attributes plays important roles in establishing areas for agricultural development and sustainable land management. Such information is particularly required in the study area in Yangtze Middle and Lower Plain and Yangtze River Delta

PT

which has a long agricultural history and rich products, such as rice and cotton. It is an important

RI

national grain production base and is known as “the land of fish and rice”. This area covers the

SC

municipality of Shanghai, Zhejiang Province and Jiangsu Province with a total area of about 210,000 km2. From north to south, this area exhibits a diverse climate, topography, soil type and landscape

NU

pattern. Thus, it is necessary to understand the spatial distribution of soil properties in this area.

MA

To understand the continuity and gradual characterization of key soil properties for agricultural development such as soil organic matter (SOM), soil pH and soil texture at a fine spatial resolution, a lot of soil surveys need to be carried out. Traditionally, soil mapping is based on an arduous task of

PT E

D

data collection, field survey, interpretation, field checking, demarcation and mapping (Sculla et al., 2003). In the recent 20 years, rapid development of precision agriculture, environmental management and eco-hydrological modelling require accurate soil maps and soil properties in raster format.

CE

Traditional soil mapping cannot meet this demand any more. To solve this issue, soil scientists all over

AC

the world used Digital Soil Mapping (DSM) methods by developing numerical models from soil observation coupled with environmental variables and statistical models (Lagacherie and McBratney, 2006; Minasny and McBratney, 2016). With the development of DSM since the 1990s, an empirical model was introduced by McBratney et al (2003), which formalised the method for predicting the spatial distribution of soil attributes. The empirical method shows that the attributes of soil can be predicted by soil factors (s), climate (c), organism (o), relief(r), parent materials (p), age (a) and spatial position (n) with spatially correlated errors (e) (McBratney et al., 2003). The relationships between soil attributes and the above 7

ACCEPTED MANUSCRIPT environmental factors are usually analysed by a variety of numerical approaches (Minasny and McBratney, 2016). This study will map soil properties that support agricultural production (texture, pH, and soil organic matter content) for the area surrounding the Yangtze River at a relatively high resolution (90 m x 90 m). This study will use data collected from a traditional soil survey and will demonstrate that such data

PT

can be readily used to produce high resolution maps using covariates that are freely and widely

RI

available. We used a decision tree model (Cubist) and Regression Kriging (RK) method to predict the

SC

spatial distribution of soil attributes. We demonstrated the efficiency and accuracy of the model for

NU

mapping soil properties over a relatively large area in China.

MA

2. Materials and methods

D

2.1 Study area

PT E

The study area is located in the eastern China, covering the municipality of Shanghai, Zhejiang Province and Jiangsu Province. It is delimited by longitude 115°59'27.81" E – 123° 1'49.09" E and latitude 27°7'7.84" N – 35°5'33.73"N with a total area of about 21.08×104 km2. The study area is in the

CE

subtropical and temperate zones of monsoonal climate. The mean annual precipitation is around 704–

AC

2,000 mm, with mean annual maximum and minimum temperatures of 13.6 and 18˚C. The northern and southern parts of the Yangtze River have distinct land use patterns. The northern parts are mostly agricultural areas, while the southern parts are dominated by forests. Forests occupy 26% of the total area particularly in mountainous areas. A small portion of grassland (3%) is distributed in patchy areas of river deltas. In the northern part, dryland agriculture, which accounts for 39% of the total area, is mainly found in the northwest and northeast of the study area. Dryland farms are distributed along the alluvial plains

ACCEPTED MANUSCRIPT (i.e., Subei Plain and Huanghuai Plain). Meanwhile, in the upland, the soils are mainly tidal soils and solonchaks which are naturally alkaline (National Soil Survey Office, 1994). Irrigated fields, which are mostly located in the north and associated with the rivers and lakes of the Yangtze Delta plains, are the second largest land use accounting for 30% of the area. The soils in these areas are paddy soils.

PT

2.2 Soil samples and soil analysis A total of 2,895 topsoil samples (0–20 cm) were available from the first National Soil Survey database.

RI

Samples were collected from October 2006 to February 2008 by environmental monitoring station of

SC

the municipality of Shanghai, Zhejiang Province and Jiangsu Province. The samples were collected using an 8 km × 8 km grid across the northern plains (in agricultural areas) and a 16 km × 16 km grid

NU

in the southern mountainous areas considering land use, topography (Fig 1.) and soil type. Sample

MA

locations were recorded using a hand-held global position system (GPS) device (Garmin Map 76, Garmin Corporation). At each pre-defined location, soils in top layer were collected from 6 to 8 points at each location from an area of about 1–2 km2. The samples were then composited, and divided into

AC

CE

PT E

D

parts of 1–2 kg each. Only one of the parts was brought back to the laboratory for analysis.

Figure1. A total number of 2895 sampling sites, elevation and land use in the Eastern China

ACCEPTED MANUSCRIPT All samples brought back to laboratory, air-dried at room temperature (20–22˚C), removed stones and debris, and then sieved to pass a 2 mm sieve. A portion of the soil sample (about 100 g) were ground in an agate grinder and sieved through a 0.149 mm mesh. The prepared soil samples were then stored in polyethylene bottles for analysis. Soil pH was measured by pH meter (Sartorius Basic pH meter PB10) with a soil/water ratio of 1:2.5. SOM content was determined using the dry ashing method. Soil

PT

texture was determined using the hydrometer method (Agricultural Chemistry Committee of China,

RI

1983), with clay defined as particles < 2 μm, silt (2 – 50 μm), and sand (50-2000 μm).

SC

2.2.1 Soil texture data transformation

Since soil texture (sand, silt, and clay content) is a compositional data, it needs to be transformed

NU

before statistical analysis (Aitchison, 1986; Huang et al., 2014). The Log-ratio transformation method

MA

transforms the compositional data into a logarithmic ratio of components and the results will be normally distributed. In this study, we used the following formula to transform the soil texture data: ln 𝑥𝑖 , 𝑖 = 1,2,3, … , 𝑛 100

(1)

PT E

D

𝑦𝑖 =

CE

The back-transformed formula is calculated as:

𝑥𝑖 =

𝑒𝑥𝑝 𝑦𝑖 𝑛 ∑𝑖=1 exp 𝑦𝑖

(2)

AC

where xi is the i-th compositional data (sand, silt, and clay), yi is the corresponding transformed compositional data.

2.2.2 Soil Organic Carbon Density (SOCD) calculation To calculate soil carbon stock, soil organic carbon density (SOCD) is calculated as follows: SOCD = SOM ∗ 0.58 ∗ BD ∗ h/100

(3)

ACCEPTED MANUSCRIPT where the unit of SOCD is kg m-2, SOM is the soil organic matter content of topsoil (g kg-1), BD is bulk density (g cm-3), 0.58 is Van Bemmelen conversion ratio and h is the thickness of soil horizon (cm). Since we do not have gravel measurement, we assume that the topsoil contains negligible amount of coarse materials, which is often true.

PT

2.3 Environmental covariates The digital soil maps predict soil properties along with their prediction intervals at a spatial resolution

RI

of 90 m, based on the Shuttle Radar Topography Mission (SRTM) data. A number of environmental

SC

covariates were collected in this study including derivatives of digital elevations models (DEM), Landsat 7 images, mean annual temperature (MAT), total annual precipitation (TAP) and sentinel-1

NU

imagery (Table 1.). All covariates were transformed to a common spatial resolution of 90 m using

MA

ArcGIS 10.0.

Table 1. List of the environmental covariates and their source and resolution Environmental covariates Landsat 7 Normalized difference vegetation index (NDVI) Digital elevation Elevation Slope model terrain ruggedness index(TRI) topographic wetness index (TWI) terrain position index (TPI) Multiresolution Index of Valley Bottom Flatness Multiresolution Ridge Top Flatness (MrRTF) (MrVBF) stream power index (SPI) Melton ruggedness number (MRN) Radar sentinel-1 imagery Climate mean annual temperature (MAT) total annual precipitation (TAP)

AC

CE

PT E

D

Type of data Remote sensing

Source Landsat Landsat SRTM SRTM SRTM SRTM SRTM SRTM SRTM SRTM SRTM SAR WorldClim WorldClim

Original 30m Resolution 30m 30m 30m 30m 30m 30m 30m 30m 30m 30m 10m 1km 1km

2.3.1 The digital elevation model (DEM) The DEM was derived from SRTM data at a spatial resolution of 30m. The DEM was resampled (coarse gridded) to 90 m for this study. Terrain attributes, i.e. aspect, terrain ruggedness index(TRI), terrain

ACCEPTED MANUSCRIPT position index (TPI), topographic wetness index (TWI), Multiresolution Index of Valley Bottom Flatness (MrVBF), Multiresolution Ridge Top Flatness (MrRTF), stream power index (SPI) and Melton ruggedness number (MRN) were derived from the DEM using SAGA GIS. The spatial variations in soil properties are scale-dependent and the pattern of distribution may be different at different spatial scales (Corwin et al., 2006; Vanwalleghem et al., 2013). This scale-

PT

dependent effects can also affect the relationship between soil properties and environmental

RI

covariates (Hu and Si, 2013; Zhou et al., 2016). In order to deal with scale-dependent variation in

SC

environmental data and soil properties, 4 elevations and slopes of different resolution were transformed from the 90 m DEM by using the discrete wavelet method via the wavelet toolbox in

NU

MATLAB 8.1. Because wavelet method decomposes objective according to 2N, so the 4 resolutions are 180m (21*90m) , 360m (22*90m), 720m (23*90m) and 1440 m(24*90m). The wavelet method also filter

MA

the DEM from noise.

D

2.3.2 Landsat 7

PT E

Landsat 7 with the Enhanced Thematic Mapper Plus (ETM+) has eight-band multispectral data. ETM+ bands 1 to 7 have a spatial resolution of 30 meters and the panchromatic band 8 has a resolution of

CE

15 meters. Band 1, 2, 3 refer to the blue, green and red, and band 4 is near infrared. The ratio between (b4-b3) and (b4+b3) is called the Normalized Difference Vegetation Index (NDVI). NDVI is an indicator

AC

of the vegetation. 2.3.3 MAT and TAP

The mean annual temperature (MAT) and the total annual temperature (TAP) were obtained from the WorldClim-Global Climate Data with a spatial resolution about 1 km. 2.3.4 Sentinel-1 imagery

ACCEPTED MANUSCRIPT Sentinel-1 is a two satellite constellation which carries a C-SAR sensor, and offers medium and high resolution imaging in all weather conditions. We used two bands from the Sentinel-1 imagery: the VV (vertical transmit/vertical receive) and VH (vertical transmit/horizontal receive) from the radar. The median value for the images from year 2015-2016 was used.

PT

2.4 Prediction method During last decades, DSM techniques have developed from simple linear models to complex machine

RI

learning techniques (Minasny and McBratney, 2016). They include: 1) geostatistical methods, such as

SC

Ordinary Kriging, Co-Kriging and Regression Kriging; 2) linear models, Generalised Linear Models (GLM) (Poggio et al.,2013) and Generalised Addictive Models (GAM) (Scull et al., 2003); 3) Data mining

NU

models, such as Cubist, Random Forests, artificial neural networks (Malone et al., 2009), Genetic

MA

Algorithm (Nelson and Odeh, 2009) and support vector machine (Ahmad et al., 2010; Ballabio, 2009); 4) Expert knowledge models, such as Bayesian inference model and SoLIM (Zhu et al.,1996; Zhu et al.,

D

2001). In this study, we used a widely-used decision tree model called Cubist (Kidd et al., 2015).

PT E

2.4.1 Cubist model

Cubist is a rule-based model which has been found to be quite effective in digital soil mapping. It is a

CE

tree model algorithm based on the M5 theory (Quinlan, 1992) and partition the predictor variates into different subsets according to “if-then” rules. In this study, we fit the Cubist model in the R statistical

AC

software. Cubist algorithm has three parameters: rules, committees and neighbours. We used the default number of rules and 10 committees, which mean that 10 boosting iteration should be used to get the optimal prediction and the contribution of the environmental predictors. 2.4.2 Regression Kriging Regression Kriging is a hybridised spatial interpolations modelling approach that combines a regression modelling of the dependent variable on environmental covariates (such as DEM, Landsat, MAT and TAP) and interpolated regression residuals. The Regression Kriging method has been

ACCEPTED MANUSCRIPT reported to improve model performance (Odeh et al., 1995; Hengl et al., 2004). In this study, the residual values obtained from Cubist model prediction were interpolated using simple kriging, assuming the mean of zero. We run punctual kriging using local variograms using the VESPER 1.5 Software (Minasny et al., 2002). The final Regression Kriging results were obtained from the sum of Cubist prediction value and the interpolated residuals. The advantage of RK over linear mixed model

PT

is the ability to use to a broader range of regression techniques.

RI

2.4.3 Model training and validation

SC

The data was partitioned into a training set and a validation set to assess the quality of the prediction model. 70% samples were selected randomly as a training set and the remaining as a validation set.

NU

The model performance was measured by the root mean square error (RMSE), coefficient of

MA

determination (R2), bias and Lin’s concordance correlation coefficient (Lin, 1989). 2.5 Uncertainty Analysis

D

In the modelling process, model structure, model parameter and model input can be the three main

PT E

sources which can lead to the output uncertainty (Huang et al., 2015). In this study, we use the empirical uncertainty quantification through fuzzy clustering and cross validation method to quantify

CE

the prediction interval (Malone et al., 2011). In this approach, prediction interval is expressed in the form of quantiles on the basis of the underlying distribution of regression kriging residuals. It includes

AC

4 main steps: 1) Leave-one-out cross validation (LOCV) of estimate of error from regression kriging; 2) Partition the environmental data into clusters using fuzzy k-means method with extragrades method; 3) Construct cluster prediction interval (PI) for each of the clusters and derive uncertainty; 4) Project the clusters spatially and map spatial predictions and uncertainties (upper and lower prediction intervals). We refer to Malone et al. (2011) for a detail description of this method. 3 Results and discussion 3.1 The selection of optimal model

ACCEPTED MANUSCRIPT Table 2 shows the performance of Cubist model on calibration data (CC), Cubist model on validation data (CV) and Regression Kriging on validation data (RKV) of soil properties across the study area. Clayt, sandt and siltt are the log-transformed results of soil texture data in Table 2. According to the model performances in Table 2, the R2 of Cubist model on validation data (CV) for SOM, SOCD, pH and soil texture (clay, silt and sand) are 0.25, 0.20, 0.64, 0.24, 0.26 and 0.18

PT

respectively. Including kriging of the Cubist models residuals, we obtain the RK predictions, which

RI

show R2 increase to 0.31, 0.25, 0.69, 0.40, 0.11 and 0.33, respectively. The percentage of R2 increases

SC

are 24%, 25%, 8%, 67%, -58%, and 83%. It is evident that Regression Kriging model improves the

NU

prediction capability for most soil properties, except for silt.

MA

Table 2. A comparison of the performance of the prediction of various soil properties using Cubist calibration (CC), Cubist validation (CV), and regression kriging on validation (RKV)

t

D

Concordance CC CV RKV 0.45 0.39 0.49 0.41 0.33 0.45 0.83 0.78 0.82 0.37 0.34 0.59 0.30 0.28 0.28 0.45 0.28 0.50

PT E

RKV 0.31 0.25 0.69 0.40 0.11 0.33

CC 1.05 0.11 0.77 0.48 0.38 0.34

RMSE CV 1.12 0.12 0.89 0.47 0.40 0.43

RKV 1.08 0.12 0.82 0.41 0.42 0.39

CC -0.15 -0.01 -0.03 0.03 0.08 0.04

bias CV -0.09 -0.01 -0.03 -0.00 0.09 0.06

RKV 0.04 -0.00 -0.01 -0.02 0.03 0.01

CE

SOM SOC pH D Clayt Siltt Sand

CC 0.37 0.31 0.71 0.30 0.29 0.37

R2 CV 0.25 0.20 0.64 0.24 0.26 0.18

AC

Figure 2 further visually compares the distribution of SOM, SOCD, pH and soil texture (clay, silt and sand) on the validation dataset using the Cubist model (CV) and regression kriging model (RKV). The xaxis is observed value and the y axis is predicted value. The line is 1:1 line. In Figure 2 the observed values are closer to predicted values using RK compared to only using Cubist for all soil properties except silt. Regression kriging takes into account the spatial autocorrelation of the residuals which cannot be explained by the Cubist model and environmental covariates, and the unbiased spatial estimation with minimized variance (Meng et al., 2013), which make the accuracy of RK outperforms Cubist model.

ACCEPTED MANUSCRIPT Thus, we conclude that the Regression Kriging model is an optimal model for predicting soil properties in this area. Generally all prediction is quite good with best prediction in soil pH (R2 = 0.71) followed by sand (R2 = 0.37), SOM (R2 = 0.37), SOCD (R2 = 0.31), clay (R2 = 0.30) and silt (R2 = 0.29). pH is well predicted because elevation and soil type are the dominant influencing factors. Soil in the north is dominated

PT

by tidal soils and paddy soils which are neutral to slight alkaline and by red acidic soils in the south

RI

thus pH increases from south with high elevation to north with low elevation. The local variation of

SC

soil pH due to rather than cultivation and land use is further taken into account when we used local kriging of the residuals.

NU

Although prediction of SOM is generally good (R2 = 0.37, RMSE = 1.05), it is worth noticing that the

MA

prediction fail to replicate some sites in the southeast with high concentration (> 5 g kg-1). This can be due to local management such as ploughing, crop rotation and straw returning to field, which improve the concentration of SOM in the study area compared to the natural elements (Li et al., 2008). Such

PT E

D

local influences cannot be captured by covariates of 90 m and also by kriging of the residuals.

AC

CE

Cc

Cv

RKv

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

PT

ACCEPTED MANUSCRIPT

RI

Fig 2. The scatters of Cubist calibration (CC), Cubist validation (CV) and Regression Kriging on validation

SC

(RK)

NU

3.2 The importance of environmental covariates

Table 3 shows the contribution of all the environmental covariates in Cubist model for the prediction

MA

of soil properties. We just list the top 10 important covariates in the model among the 26 covariates. For SOM and SOCD prediction, we can see that elevation at 4 levels of decomposition (1-4), Landsat

D

band 3, MAT, TAP and land use are the contributing factors in the model. Elevation 1 is the most

PT E

important predictor in SOM model where it is used 90% in the model, however Total Annual Precipitation (TAP) is the most important predictor in the SOCD model with 90% use. Band 4 and NDVI

CE

contribute only in SOM model and MrVBF and TWI only in SOCD model. The contribution of elevation at 4 levels show that SOM and SOCD are closely related to the elevation at various spatial resolution.

AC

Elevation at multiple resolutions may reflect local physiography, which can be important for predicting soil attributes.

In terms of soil pH, elevation 1-3, slope 1-2, band 2-3, MAT and TWI play important roles in the model. The largest contributor one is TAP (85%) and slope 2 (48%). The percentage of other factors rise from 52% to 79%. Regarding to soil texture, only elevation 1-3, MAT and TAP play roles in the three texture models. Elevation 4 contributes in clay model (39%) and silt model (31%). MrVBF only appears in silt model

ACCEPTED MANUSCRIPT (27%) and TWI appears in silt (27%) and sand model (53%). Band 2-4 contribute in clay model, band 4-7 in silt model, and band 3 in sand model. From the six soil attributes, it is obvious that DEM and climate factors play significant roles in model prediction. Radar images from Sentinel and other topographic indices are not found to be important

PT

predictors.

pH Conds Model 56% 52% 60% 19% 57% 48%

13%

CE

56% 65%

29% 53%

D

84% 78%

67% 63%

79% 90% 65%

10% 91%

72% 85%

58%

12%

79%

PT E

56% 30% 12%

51%

16% 21% 13% 12%

MA

22% 63% 55%

Clayt Conds Model 26% 6% 54% 8% 27% 39% 13%

SC

15%

Model 76% 61% 64% 47%

NU

SOCD Conds 48%

28% 93% 2%

63% 80%

Siltt Conds Model 6% 68% 7% 67% 30% 31%

Sandt Conds Model 66% 60% 9% 45% 22% 41%

6% 45%

24% 89% 21%

34% 54% 38% 66% 69% 67% 27% 27%

41% 33% 84% 20%

12%

56%

AC

Elevation1 Elevation2 Elevation3 Elevation4 Slope1 Slope2 Slope3 Slope4 Band2 Band3 Band4 Band5 Band7 MAT TAP MrVBF MrRTF TWI NDVI Land use Sentinel band1 Sentinel band2 TPI TRI SPI MRN

SOM Conds Model 18% 90% 85% 54% 50% 60%

RI

Table 3. Contributions of the environmental covariates predictors in Cubist model

3.3 Prediction maps Table 4. The statistical characteristics of observed values, Cubist prediction and RK prediction on the map. Min

Mean

Max

39% 61% 39% 53%

ACCEPTED MANUSCRIPT

MA

NU

SC

RI

PT

Observed Cubist RK Observed Cubist RK Observed Cubist RK SOM 0.07 0.62 -0.12 2.17 2.21 2.34 13.00 8.07 12.03 SOCD 0.01 0.15 0.04 0.31 0.34 0.35 0.74 0.56 0.82 pH 3.30 4.14 3.50 6.77 6.63 6.47 9.30 9.02 9.88 Clay 0.60 4.23 4.00 14.31 12.16 12.8 43.20 26.93 27.11 Silt 2.20 15.80 15.92 39.97 39.46 39.35 85.40 58.69 58.43 Sand 1.00 15.34 14.90 45.69 48.42 48.48 96.90 86.35 86.16

(b) SOCD

AC

CE

PT E

D

(a) SOM

(c) pH

(d) clay

RI

PT

ACCEPTED MANUSCRIPT

(f) sand

SC

(e) silt

NU

3. Distribution of two models of soil properties in the topsoil across the study area

The spatial distribution of Cubist model and Regression Kriging method for soil properties across the

MA

study area are shown in Fig.3 and the statistical parameters for Cubist and RK maps are shown in Table 4. Generally both maps are similar, and the differences are more reflected in local pattern.

PT E

D

SOM varies largely across the study area and increases from north to south (Fig 3.a). In the dryland and irrigated area with low elevation, SOM is relatively low. However, in the forest area with high elevation, SOM is relatively high. This is mainly due to the land use type and elevation variation. In the

CE

soil sampling area with high elevation, the low temperature and less human disturbance are indicators

AC

of SOM accumulation, and vice versa. The dry lands and irrigated areas have significant human disturbances with little accumulation of plant residues. The forest area has a vegetation cover and a high humus content in the top layer. There is slight differences between the Cubist and RK prediction of SOM in this area. The RK prediction shows a large extent of data (-0.12~12.03) (Table 4.) and large variation with higher SOM in the west part and lower in the northern under the same elevation (Fig 3.a). To show the differences more specifically, we selected a typical agricultural area, which across the Yangtze River and includes the Taihu Lake (Fig 4. - Fig6.).

ACCEPTED MANUSCRIPT Figure 4 shows the spatial distribution of SOM in selected area. It is obvious that SOM concentration is high in Taihu Lake Basin. Taihu Lake Basin is one of the highly intensive use of agricultural areas in China, where a lot of fertilizers have been applied and agricultural land use has changed. From 1980s, the large amounts of food crops planting have shifted into the planting of economic crops such as vegetables which needed abundant organic fertilizer and the agricultural management measures such

PT

as straw turnover have been promoted (Liu et al., 2009). These elements all lead to the high SOM

RI

content around Taihu Lake.

SC

Topsoil SOCD is calculated from SOC and soil bulk density, and shows variation similar to SOM. SOCD is lower in the dryland area and higher in the irrigated field and forest area (Fig.3b). Table 4 shows

NU

that the variation of SOCD from the RK model (0.04~0.82) is higer than that of Cubist model (0.15~0.56). Moreover, the spatial distribution of SOCD is more detailed in the RK model. RK shows

AC

CE

PT E

D

MA

that some parts have lowest SOCD in drylands which is not shown in Cubist model.

Fig 4. Distribution of SOM in the topsoil across the selected area in Taihu Lake Basin using the RK model

Soil pH shows obvious differences from south to north and across acidic soil, neutral soil and alkaline soil in the study area (Fig.3c). Soils are mostly acid in the forest area except some costal part.

ACCEPTED MANUSCRIPT Conversely, in the dryland and irrigated areas, soils are neutral to alkaline. The regional differences are mainly due to different parent materials. The tidal soils and solonchaks in northern part are alkaline while the southern forest is dominated by red acidic soil. The data variation of RK (4.14~9.02)

MA

NU

SC

RI

PT

is larger than that of Cubist method (3.50~9.88).

Fig 5. Distribution of soil pH in the topsoil across the selected area in Taihu Lake Basin using the RK

D

model

PT E

In Figure 5, the land use in the selected area is dominated by irrigated field and paddy soils, but its soil pH distribution is quite variable. Soils in Taihu Lake Basin in the south are acid and soils in costal part

CE

(north) are neutral and alkaline. The reason of low soil pH around Taihu Lake are the widely and increased use of chemical fertilizer especially nitrogen fertilizer. In addition, with the rapid

AC

development of industry, a lot of large chemical companies are distributed around Taihu Lake. The large emission of SO2 and NO2 into atmosphere caused acid rain which further lowered the soil pH (Liu et al., 2006).

Soil texture varies largely across the study area. Basically, in the southern part, soil texture has small variation and the silt and sand contents are high (Fig 3.d - f). The southern forest area is dominated by red soils, with sandy texture and low clay content (<15%). The northern part shows a contrasting texture.

RI

PT

ACCEPTED MANUSCRIPT

NU

RK model

SC

Fig 6. Distribution of clay content in the topsoil across the selected area in Taihu Lake Basin using the

In the areas next to Yangtze River Plain, clay content is high due to the marine sediments (Fig 6.).

MA

Paddy fields are mostly found in this area as its soil has a high silt content which is preferred in rice cultivation. The northeastern and northwestern parts have high sand content due to the tidal soil and

D

solonchaks. Compared Table 4 and Figure 3 d-f, it is clear that there is no obvious differences between

PT E

Cubist and RK model. As texture is governed mostly by environmental factors rather than local land use, the Cubist prediction extent of clay, silt and sand are quite similar to the data extent of soil texture

CE

in RK model (Table 4.).

AC

4 Uncertainty analysis

After assessing the soil attributes, the uncertainty analysis provides the critically best and worse conditions of soil attributes in a landscape, knowing the upper and lower prediction bounds (Xiong et al., 2015). Model predictions are associated with the prediction interval or a range with a defined level of confidence. Figure 7 shows the variability of SOM in the topsoil across the study area with a lower prediction limit, an upper prediction limit, and the predicted map and its prediction limit range. These prediction limits

ACCEPTED MANUSCRIPT constitute a 90% prediction interval. For the most part, 90% of the time, an observed value fitted within its well-defined predicted prediction interval (Malone et al., 2014). In Figure 7, the predicted average content of SOM was 2.35 g kg-1. We are 90% confident that the true

CE

PT E

D

MA

NU

SC

RI

PT

average of a given SOM is between 1.41 and 3.06 g kg-1 based on the lower and upper prediction limits.

AC

Fig 7. SOM predictions and prediction limits derived using a Cubist model together with LOCV and fuzzy classification 5 Conclusions

In this study, we predict the spatial distribution variation of soil properties with DEM, Landsat, Climate and Radar data as environmental factors using Cubist decision tree and Regression Kriging method. We can conclude that:

ACCEPTED MANUSCRIPT 1) DEM and climate factors play significant roles in model prediction. The different resolutions of elevation characterise local physiography and was found to be important predictors. Radar images from Sentinel and other topographic indices were not found to be important predictors in this study area. 2) Local RK is an efficient tool to predict the spatial distribution of soil attributes and the accuracy and

PT

precision of prediction of local RK in this study appears to be superior to just using Cubist regression

RI

tree method The Cubist model can only model the correlation between environmental covariates and

SC

soil attributes without the spatial autocorrelation of the residuals.

3) This study showed that digital soil mapping technique using widely-available covariates can be used

NU

efficiently to map a range of soil properties at a relatively high resolution (90 m) over a large area

MA

(210,000 km2). This study can be readily expanded to the whole China using existing soil data and contribute to the GlobalSoilMap project (Arrouays et al., 2014). In addition, with increasing availability

PT E

a finer spatial resolution (30m).

D

of fine resolution covariates and computing power, future studies will explore the advantage of RK at

Acknowledgments

CE

This study was funded by the scholarship from the Jiangsu government for overseas studies (JS-2014322) and the Environmental Protection Research Project of Jiangsu Province (Grant no. 201264). We

AC

appreciate all the colleagues who collected and analyzed the soil samples.

Reference Ahmad, S., Kalra, A., Stephen, H., 2010. Estimating soil moisture using remote sensing data: A machine learning approach. Advances in Water Resources 33, 69-80. Aitchison, J., 1986. The statistical analysis of compositional data. Chapman and Hall Inc., London-New York, pp. 58-61.

ACCEPTED MANUSCRIPT Arrouays, D., Mckenzie, N., Hempel, J., McBratney A.B., Richer de Forges, A.C., McBratney, A.B., 2013. GlobalSoilMap: Basis of the global spatial soil information system. CRC Press, the Netherlands. Ballabio, C., 2009. Spatial prediction of soil properties in temperate mountain regions using support vector regression. Geoderma 151 (3–4), 338–350. Corwin, D.L., Hopmans, J., Rooij, G.H.de, 2006. From field- to landscape-scale vadose zone processes:

PT

scale issues, modelling, and monitoring. Vadose Zone 5 (1), 129–139.

RI

Hengl, T., Heuvelinkb, G.B.M., Stein, A., 2004. Generic framework for spatial prediction of soil variables based on regression-kriging. Geoderma 120, 75-93.

SC

Hu, W., Si, B.C., 2013. Soil water prediction based on its scale-specific control using multivariate

NU

empirical mode decomposition. Geoderma 193, 180–188.

Huang, J., Subasinghe, R., Triantafilis, J., 2014. Mapping particle-size fractions as a composition using

MA

additive log-ratio transformation and ancillary data. Soil Science Society of America Journal 78(6), 1967-1976.

D

Huang, J., Zare, E., Malik, R.S., Triantafilis, J., 2015. An error budget for soil salinity mapping using

PT E

different ancillary data. Soil Res. 53, 561–575. Kidd, D., Webb, M., Malone, B., Minasny, B., McBratney, A., 2015. Eighty-metre resolution 3D soilattribute maps for Tasmania, Australia. Soil Research 53(8), 932-955.

CE

Lagacherie, P., McBratney, A.B., 2006. Spatial soil information systems and spatial soil inference

AC

systems: perspectives for Digital Soil Mapping, in Lagacherie, P., McBratney, A.B., Voltz M. (Eds.), Digital Soil Mapping: An Introductory Perspective. Elsevier Inc., Amsterdam, pp. 3-22. Li, X.H., Hao, M.D., Wang, Z.H., Li, L.L., 2008. Factors affecting soil organic carbon in cropland and their regulation. Agricultural Research in the Arid Areas 26(3), 176-181. Lin, L.I., 1989. A concordance correlation coefficient to evaluate reproducibility. Biometrics 45(1), 255268. Liu, F.C., Shi, X.Z., Yu D.S., 2006. Spatial and temporal variability of sol acidity in typical areas of Taihu Lake region in the last 20years. Resources and Environment in the Yangtze Basin 15(6), 740-744.

ACCEPTED MANUSCRIPT Liu, Y.M., Jiang N., 2009. Change of agriculture soil organic carbon in the Tai Lake Region, China. Soil 41(5), 715-718. Malone, B.P., McBratney, A.B., Minasny, B., Laslett, G.M., 2009. Mapping continuous depth functions of soil carbon storage and available water capacity. Geoderma 154 (1–2), 138–152. Malone, B.P., McBratney, A.B., Minasny, B., 2011. Empirical estimates of uncertainty for mapping

PT

continuous depth functions of soil attributes. Geoderma 160, 614-626.

RI

Malone, B., Minasny, B., Odgers, N., McBratney, A.B., 2014. Using model averaging to combine soil property rasters from legacy soil maps and from point data. Geoderm 232-234, 34-44.

SC

McBratney, A.B., Mendonca Santos, M.L., Minasny, B., 2003. On digital soil mapping. Geoderma 117,

NU

3-52.

Meng, Q.M., Liu, Z.J., Borders, B.E., 2013. Assessment of regression kriging for spatial interpolation –

MA

comparisons of seven GIS interpolation methods. Cartography and Geographic Information Science 40, 28-39.

D

Minasny, B., McBratney, A.B., Whelan, B.M., 2002. Vesper (Variogram estimation and spatial

PT E

prediction plus error) version 1.6. Aust. Ctr. for Precision Agric., Sydney. Minasny, B., McBratney, A.B., 2016. Digital soil mapping: A brief history and some lessons. Geoderma 264, 301-311.

CE

National Soil Survey Office, 1994. Soil Taxonomy (in Chinese). China Agriculture Press, Beijing.

AC

Nelson, M.A., Odeh, I.O.A., 2009. Digital soil class mapping using legacy soil profile data: a comparison of a genetic algorithm and classification tree approach. Australian Journal of Soil Research 47, 632649.

Odeh, I.O.A., McBratney, A.B., Chittleborough, D.J., 1995. Further results on prediction of soil properties from terrain attributes: heterotopic cokriging and regression-kriging. Geoderma 67, 215–226. Poggio, L., Gimona, A., Brewer, M.J., 2013. Regional scale mapping of soil properties and their uncertainty with a large number of satellite-derived covariates. Geoderma 209-210, 1-14.

ACCEPTED MANUSCRIPT Quinlan, J.R., 1992. Learning with continuous classes. In Proc. of the Fifth Australian Joint Conference on Artificial Intelligence World Scientific, Singapore, pp. 343–348. Scull, P., Franklin, J., Chadwick, O.A., McArthur, D., 2003. Predictive soil mapping: a review. Progress in Physical Geography 27(2), 171-197. Vanwalleghem, T., Stockmann, U., Minasny, B., McBratney, A.B., 2013. A quantitative model for

PT

integrating landscape evolution and soil formation. J. Geophys. Res. Earth Surf. 118 (2), 331–347.

RI

Xiong, X., Grunwald, S., Myers, D.B., Kim, J., Harris W.G., Bliznyuk N., 2015. Assessing uncertainty in soil organic carbon modelling across a highly heterogeneous landscape. Geoderma 251-252, 105-

SC

116.

NU

Zhou, Y., Biswas, A., Ma, Z., Lu, Y., Chen, Q., Shi, Z., 2016. Revealing the scale-specific controls of soil organic matter at large scale in Northeast and North China Plain. Geoderma 271, 71–79.

Ecological Modelling, 90, 123-145.

MA

Zhu, A.X., Band, L.E., Dutton, B., Nimlos, T.J., 1996. Automated soil inference under fuzzy logic.

D

Zhu, A.X., Hudson, B., Burt, J., Lubich, K., Simonson, D., 2001. Soil mapping using GIS, expert knowledge,

AC

CE

PT E

and fuzzy logic. Soil Science Society of America Journal 65:1463-1472.

ACCEPTED MANUSCRIPT Highlights

 DEM and climate factors play significant roles in model prediction.  Local Regression Kriging is an efficient tool to predict the spatial distribution of soil attributes. 

The accuracy and precision of prediction of local RK in this study appears to be superior to

PT

just using Cubist regression.  Digital soil mapping technique using widely-available covariates can be used efficiently to map

AC

CE

PT E

D

MA

NU

SC

RI

a range of soil properties.