Spatial-Temporal Changes of Soil Organic Carbon Content in Wafangdian , China

Soil organic carbon (SOC) plays an important role in soil fertility and the global carbon cycle. A better understanding of spatial-temporal changes of SOC content is essential for soil resource management, emission studies, and carbon accounting. In this study, we used a boosted regression trees (BRT) model to map distributions of SOC content in the topsoil (0–20 cm) and evaluated its temporal dynamics from 1990–2010 in Wafangdian City, northeast of China. A set of 110 (1990) and 127 (2010) soil samples were collected and nine environment variables (including topography and vegetation) were used. A 10-fold cross-validation was used to evaluate model performance as well as predictive uncertainty. Accuracy assessments showed that R2 of 0.53 and RMSE (Root-mean-square error) of 9.7 g·kg−1 for 1990, and 0.55, and 5.2 g·kg−1 for 2010. Elevation and NDVI (Normalized Difference Vegetation Index) were the two important variables affecting SOC distribution. Results showed that mean SOC content decreased from 19 ± 14 to 18 ± 8 g·kg−1 over a 20 year period. The maps of SOC represented a decreasing trend from south to north across the study area in both periods. Rapid urbanization and land-use changes were accountable for declining SOC levels. We believe predicted maps of SOC can help local land managers and government agencies to evaluate soil quality and assess carbon sequestration potential and carbon credits.


Introduction
Soil organic carbon (SOC) is the largest terrestrial carbon pool and plays an important role in the global carbon cycle [1][2][3][4].It determines the fate of a terrestrial ecosystem as a potential sink or a source of greenhouse gases [5].In addition, SOC is closely related to soil quality, nutrient availability and crop yield [5].Its decomposition can increase atmospheric CO 2 levels [6], and its sequestration can significantly mitigate global climate change [7].Changes in land use from native vegetation to agriculture are able to cut down on SOC levels [8,9] and increase atmospheric carbon dioxide concentration [10].Therefore, studies focusing on the impact of land-use change and environmental variables controlling spatial-temporal dynamics of SOC are important.Especially in China, which is experiencing population expansion, assessment of spatial-temporal changes of soil properties (e.g., SOC) is considered urgent to make policies to regulate actual needs and ecological demand [7,[11][12][13][14].Moreover, reliable, robust and cost-efficient methods for estimating SOC content and its temporal changes are quite essential.Spatial-temporal variability of SOC content can be estimated by using statistical models that quantify the relationship between SOC content and environmental variables.The number of SOC measurements used in predicting is an issue because it depends on many factors such as scales and costs.It is unlikely to collect sufficient soil samples to reflect the spatial-temporal changes of SOC content at the country scale [7].Digital soil mapping (DSM) can predict soil properties over large areas based on sparse or discrete samples [14].Most DSM studies are based on Jenny's equation [15], which characterizes soils as a function of environmental variables (i.e., climate, organisms, relief, parent material and time).This concept has been employed in numerous DSM techniques such as multiple linear regression [16], linear mixed models [17], random forest (RF) model [18][19][20][21], support vector machines [15], artificial neural networks [22,23], kriging [7] and the boosted regression tree (BRT) model [1,24] for SOC predictions.Despite various DSM techniques available, selecting an optimal method to predict soil properties is not always promising [1,25].
Referring to DSM techniques listed above, BRT is a relatively new tree-based model [22,25] and has rarely been used to model spatial-temporal variability of SOC [26,27].Compared to the single-tree model, BRT is more powerful and efficient [20,26].These studies show that BRT is a reliable model in predicting soil parameters.The fitted model in BRT is a simple linear combination of many trees that are fitted iteratively and boosted to reweight poorly modeled observations [22,27,28].Iterative processes improve the predictive performance and ensure robustness [20][21][22].BRT can deal with various data types, missing values, outliers, irrelevant predictors and interactions between predictors and provide access to evaluate the fitted model [22,27,28].Owing to these merits, BRT has been widely applied in various scientific fields, including soil science [22,29], remote sensing [30], environmental science [31], ecology [25] and epidemiology [32].However, BRT has not yet been applied in spatial-temporal modeling of SOC.
The aim of this study was to use BRT to evaluate spatial-temporal changes in topsoil (0-20 cm) organic carbon content in the northeast of China.The specific objectives were: (1) to map SOC content in 1990 and 2010; (2) to quantify environmental controls of SOC content; and (3) to investigate temporal dynamics of SOC content at a 20 year period (1990-2010).

Study Area
The study was conducted in Wafangdian City, which is located in the southwest of the peninsula of Liaodong in Liaoning province of China (latitude 39 • 20 -40 • 07 N; longitude 21 • 13 -122 • 16 E) (Figure 1).It covers a total area of 3826.6 km 2 of which 71% is under agriculture, 13% as garden plots and the rest as construction land used for construction purposes during rapid urbanization.Geomorphologically, the area is characterized by a complex and undulated hill systems intersected by river valleys, which is mostly agriculture dominated.Relative elevation ranges from 0 m to 830 m ASL (Above Sea Level) with an increasing trend from the southwest to the northeast of Wafangdian.The climate is a sub-humid warm temperate continental monsoon type consisting of four seasons, while, in the summer, it is hot and rainy.The mean annual temperature is 9.3 • C, and the mean annual precipitation ranges between 580 to 750 mm.Precipitation of 60%-70% in the study area is concentrated in the summer with frequent rainstorms.Dominant soils in the study area are Cambisols (58%) and Fluvisols (13%) (World Reference Base for Soil Resources (WRB) [33]).Landuse types follow the common land-use system of the northeast of China-agriculture with corn, and soybeans as main crops, and peanuts and fruits in a rotation.Much of the natural vegetation is warm-temperate deciduous broad-leaved forest where the main tree species are Pinaceae, Taxodiaceae, Cupressaceae and Ginkgoaceae.According to the vegetation type and topography, the Wafangdian City is divided into three geographic areas (Low Mountain and Hilly Area of the Southwest vs. the Hilly-Plain Transitional Area of Central vs. Mountain and Hilly Area of Eastern) (Figure 2).With the Chinese economic reform starting in the 1980s, Wafangdian City started to accelerate urbanization and economic development at an unprecedented pace.Due to the dramatic sprawl of rapid urbanization and persistent increase of population, a large number of forests, grasslands, wetlands and others were converted to construction and farmland.Thus, the 1990s was considered the background of this study.started to accelerate urbanization and economic development at an unprecedented pace.Due to the dramatic sprawl of rapid urbanization and persistent increase of population, a large number of forests, grasslands, wetlands and others were converted to construction and farmland.Thus, the 1990s was considered the background of this study.started to accelerate urbanization and economic development at an unprecedented pace.Due to the dramatic sprawl of rapid urbanization and persistent increase of population, a large number of forests, grasslands, wetlands and others were converted to construction and farmland.Thus, the 1990s was considered the background of this study.

Soil Survey Data from 1990
Soil profile data surveyed in 1990 was obtained from the Agro-Technology Extension Center of Wafangdian City.A total of 110 soil profiles were investigated from all land use and soil types (Figure 1).Soil profile data included information on parent material, soil physical and chemical properties, land use and cropping systems.However, the present study focused on SOC content of the topsoil (0-20 cm) only.

Soil Sampling in 2010
Due to rapid land-use change during the past few decades, re-sampling the same sites for a given land-use type from 1990 was not possible.Thus, we designed a purposive sampling to collect a total of 127 topsoil samples representing variations in elevation, aspect, land use and parent material (Figure 1) in 2010, following the methods described in [34].A handheld global positioning system (GPS) was used to record the geographical location of each sampling site.Soil samples from each location were a composite of five soil cores taken from four corners, and the center of a 1 m × 1 m square area.From the composite sample, 1 kg subsample was retained for laboratory analysis.SOC content of the samples was determined by a wet oxidation method [35] in the Key Laboratory of Agricultural Resources and Environment of Liaoning Province, Shenyang Agricultural University, Shenyang, China.

Predictor Variables
Environmental variables used to predict SOC consisted of terrain attributes and remote sensing (RS) data.Terrain attributes were derived from a 90 m × 90 m Digital Elevation Model (DEM) downloaded from the United States Geological Survey (USGS, Reston, VA, USA).It included elevation, slope aspect, slope gradient, topographic wetness index (TWI) and catchment area (CA) calculated in the System for Automated Geoscientific Analyses (SAGA, Hamburg, Germany) GIS [36].

Land-Use Data
Land-use maps were a vector polygon map created from Landsat images from the year 1990 and 2010.Land-use types were grouped into cultivation land, grassland, and forest land, and SOC content in each land use was derived for both periods.

Prediciton Model Based on Boosted Regression Trees
To quantify the association between environmental variables and SOC, we applied Boosted Regression Trees developed by Friedman et al. [39].To our knowledge, this study is the first attempt to identify spatial-temporal evolution of SOC using the BRT model.BRT combines two statistical techniques-boosting and regression trees-in a single algorithm.Boosting is based on the stochastic gradient improvement of decision trees by minimizing the loss function (deviation) in the splitting of each tree, and regression trees grow with recursive binary splitting of data until there is some stopping criterion.BRT is flexible to use and yet capable of dealing with complex soil-environmental responses, including interactions and nonlinearities [22,27].Compared to the traditional regression approach, BRT showed better results in terms of prediction capability [40,41].BRT performed well in soil properties modeling [42] and has been applied in various studies [27,43].The BRT script provided by Elith et al. [27] with reference package "dismo" version 0.8-17 [44] and "gbm" version 2.1 [45] were used for model calibration in the R environment [22].
Model fitting in BRT usually required the user to set up four parameters: the learning rate (LR), tree complexity (TC), bag fraction (BF) and number of trees (NT) [22,27,41].LR represented the contribution of each tree to the final fitted model [22].TC defined the predictor variables between depth of the tree and the maximum level of interaction [27].BF represented a proportion of data used in each model-the more data used, the less the randomness [26,28].Although BRT was a great extension to model over-fitting, it was still necessary to set the optimum number of trees [27].NT was determined from the combination of LR and TC.For the best predictive performance, various combinations of these parameters were tested to determine optimal settings using 10-fold cross-validation (CV) [27].We determined the optimum values for these parameters by running many simulations covering a range of LR (0.001-0.05),TC (5-10), BF (0.65-0.85) and NT (500-1000) in the two periods.The final optimal values of LR, TC, BF and NT were set as 0.0025, 9, 600, and 0.70 for the 1990 survey and as 0.0025, 9, 0.65 and 800 for the 2010 survey.
The uncertainty associated with the BRT prediction was evaluated with standard deviation (SD) derived from running the model 100 times and the map was generated as an indicator of SOC prediction uncertainty.To obtain the relative importance (RI) of each predictor, the BRT model was repeated for 100 iterations.The RI of variables were measured based on the number of times a variable was selected for modeling and weighted by the square improvement to each split and averaged across all trees [29].The RI of each variable was then scaled so that the sum added to 100 as percentages.A higher percentage of a variable indicated a stronger relative importance of this variable on the response [46].

Model Validation
The performance of the BRT model was evaluated using a 10-fold cross-validation that involved comparisons between predicted and observed SOC values [47].Cross-validation is a useful approach for model identification and assessment of model performance [22].Mean absolute prediction error (MAE), root mean square error (RMSE), coefficient of determination (R 2 ) and Lin's concordance correlation coefficient (LCCC) [48] were considered to validate the model and were calculated as follows: where a i and b i are the predicted and observed SOC values at site i, respectively; a and b are the average values of the predicted and observed SOC values; n is the number of samples; ∂ a and ∂ b are the variances of predicted and observed values; and r is the Pearson correlation coefficient between the observed and predicted values.

Descriptive Statistics of SOC Content
The statistical results of SOC content related to land-use types in 1990 and 2010 are shown in Table 1.SOC content of cultivated land in 1990 ranged from 1.4 to 46.8 g•kg −1 , with a mean of 12.4 g•kg −1 .Average SOC content of cultivated land in 2010 was 10.0 g•kg −1 (±4.2 g•kg −1 ).Compared to SOC content of cultivated land in 1990, it showed a decreasing trend in 2010, with a narrow range of values.In cultivated land, coefficient of variations (CVs) in 1990 and 2010 were 91.1% and 42.2%, respectively, and SD values were 11.3 and 4.2 g•kg −1 .In grassland and forest, SOC content showed slight changes in both periods.SOC content in all land-use types followed a lognormal distribution and was nearly symmetrically distributed around the mean SOC value in both periods.The relationship between SOC content and environmental variables in 1990 and 2010 is shown in Figure 3.

Descriptive Statistics of SOC Content
The statistical results of SOC content related to land-use types in 1990 and 2010 are shown in Table 1.SOC content of cultivated land in 1990 ranged from 1.4 to 46.8 g•kg −1 , with a mean of 12.4 g•kg −1 .Average SOC content of cultivated land in 2010 was 10.0 g•kg −1 (±4.2 g•kg −1 ).Compared to SOC content of cultivated land in 1990, it showed a decreasing trend in 2010, with a narrow range of values.In cultivated land, coefficient of variations (CVs) in 1990 and 2010 were 91.1 % and 42.2%, respectively, and SD values were 11.3 and 4.2 g•kg −1 .In grassland and forest, SOC content showed slight changes in both periods.SOC content in all land-use types followed a lognormal distribution and was nearly symmetrically distributed around the mean SOC value in both periods.The relationship between SOC content and environmental variables in 1990 and 2010 is shown in Figure 3. SOC content varied significantly among geographic areas (Figure 2, Table 2), with a decreasing trend from north to south.Hilly-Plain Transitional Area of the Central had the highest decrease of about 4 g•kg −1 .Areas in the Low Mountain and Hilly Area of the Southwest showed a slight decrease (−0.48 g•kg −1 ).In the eastern mountain area, forests were converted to plantation crops such as apples, cherries and peaches, in which SOC content decreased by 2.8 g•kg −1 .SOC content varied significantly among geographic areas (Figure 2, Table 2), with a decreasing trend from north to south.Hilly-Plain Transitional Area of the Central had the highest decrease of about 4 g•kg −1 .Areas in the Low Mountain and Hilly Area of the Southwest showed a slight decrease (−0.48 g•kg −1 ).In the eastern mountain area, forests were converted to plantation crops such as apples, cherries and peaches, in which SOC content decreased by 2.8 g•kg −1 .Linear correlations between SOC and predictors are shown in Table 3. SOC was positively correlated with elevation.The coefficient was higher in 2010 (r = 0.73) than in 1990 (r = 0.37).Slope aspect, CA and TWI had higher coefficients with SOC (0.48, −0.53 and 0.24) in 2010 compared to those in 1990.Correlations between SOC and imagery indices were significant in 1990.

Model Performance
Table 4 shows the predictive performance of BRT model evaluated with 10−fold cross-validation by 100 runs in both periods.Summary statistics showed that the BRT model had systematically higher R 2 , lower MAE and RMSE.Our results are comparable with previous studies.Ließ et al. [49] reported a superior prediction performance of BRT over Random forest (RF).Yang et al. [22] found a better performance of BRT while mapping SOC in the Qinghai-Tibet Plateau.Razakamanarivo et al. [50] concluded that the best model in predicting SOC stocks in the Malagasy Highlands was BRT.The model prediction uncertainty is shown in Figure 4. Predictions have low uncertainties with a mean SD of 1.7 and 0.7 g•kg −1 in 1990 and 2010, respectively.We calculated the SDs of validation measurements from 100 runs and found that the BRT in both periods produced low SDs of MAE, RMSE, R 2 and LCCC (Table 4).These results revealed that the BRT model had a robust performance to predict SOC distribution in both periods.It is notable that there are unmentioned uncertainties such as sampling error and laboratory measurement error associated with the prediction [51][52][53][54].However, assessments of these errors are beyond the scope of this study.The model prediction uncertainty is shown in Figure 4. Predictions have low uncertainties with a mean SD of 1.7 and 0.7 g•kg −1 in 1990 and 2010, respectively.We calculated the SDs of validation measurements from 100 runs and found that the BRT in both periods produced low SDs of MAE, RMSE, R 2 and LCCC (Table 4).These results revealed that the BRT model had a robust performance to predict SOC distribution in both periods.It is notable that there are unmentioned uncertainties such as sampling error and laboratory measurement error associated with the prediction [51][52][53][54].However, assessments of these errors are beyond the scope of this study.

Relative Importance of Environmental Variables
Variable importance revealed different dominating environmental factors of SOC in 1990 and 2010 (Figure 5).RS indices were the primary variables influencing SOC in 1990, as it was topography in 2010 (Figure 5).In 1990, NDVI (RI > 20%) in both years showed high indicating ability of the

Relative Importance of Environmental Variables
Variable importance revealed different dominating environmental factors of SOC in 1990 and 2010 (Figure 5).RS indices were the primary variables influencing SOC in 1990, as it was topography in 2010 (Figure 5).In 1990, NDVI (RI > 20%) in both years showed high indicating ability of the topsoil organic carbon in coastal mountain areas.This finding is consistent with previous studies [22,[55][56][57][58].It implies that there is a potential of remote sensing in mapping SOC distribution in coastal mountain regions with humid climates.Minnasny et al. [59] reported that such indices were not only excellent predictors, but also played positive roles in the microbial activities in soils affecting the distribution of SOC.
2016, 8, 1154 9 of 16 topsoil organic carbon in coastal mountain areas.This finding is consistent with previous studies [22,[55][56][57][58].It implies that there is a potential of remote sensing in mapping SOC distribution in coastal mountain regions with humid climates.Minnasny et al. [59] reported that such indices were not only excellent predictors, but also played positive roles in the microbial activities in soils affecting the distribution of SOC.The lower influence of NDVI on SOC in 2010 can be attributed to the heavy land-use change such as a large number of forests, grasslands and wetlands converted to agriculture from 1990 to 2010.Some studies [7,12,13] suggested that land-use change had a significant effect on SOC dynamics, since the start of the economic reforms in 1990 in Wafangdian City.This transformation weakened the intrinsic relationship between land-use and organic carbon.Yang et al. [22] concluded that band B3 had an excellent ability to predict SOC in high altitude areas (3200-3600 m).Surprisingly, our study found a lower relative importance of B3 in both periods.Although B3 was useful in distinguishing vegetation types and coverage, it showed the lowest importance.This is likely to be partly mediated by elevation and NDVI, which provide more detailed information in a landscape.
Topography has been recognized as one of the key pedogenetic factors influencing the distribution of SOC.In a high-variation of relief, topography was good predictors for mapping SOC [38,60].We found topographic variables showed higher RI in 2010 (64% in 2010, and 34% in 1990).Elevation was identified as a key topographic factor in both periods.This is consistent with previous studies [60,61].Altitude affects vegetation, local climate, and microbial activity, which have prolonged influence on SOC dynamics.TWI is related to soil moisture status, and it has been successfully used in predicting SOC in several studies [12,[62][63][64].Because of high variation in relief in the area, slope and aspect may have an effect on microclimate on local scales.We found that the influence of slope gradient was different in two periods (2% in 1990; 20% in 2010).In the study, the dataset of 1990 does not include samples in high slopes.This is because sampling is quite difficult in these areas.The purposive sampling strategy used in 2010 involving soil samples from those areas.The predictive model identified the substantial contribution of slope gradient on SOC distribution.Spatial visualization of results was essential to understand the driving forces on SOC distribution.

SOC Content Maps and Spatial-Temporal Changes
Figure 6 shows the predicted SOC in 1990 and 2010.They are means of 100 runs of BRT models.Average SOC content was higher in 1990 (20.8 ± 9.7 g•kg −1 ) than that in 2010 (18.6 ± 6.0 g•kg −1 ).The percentages of SOC content related to levels were shown in Figure 7.In both periods, SOC was high in the east and low in the southwest of Wafangdian City.High SOC content was observed in the The lower influence of NDVI on SOC in 2010 can be attributed to the heavy land-use change such as a large number of forests, grasslands and wetlands converted to agriculture from 1990 to 2010.Some studies [7,12,13] suggested that land-use change had a significant effect on SOC dynamics, since the start of the economic reforms in 1990 in Wafangdian City.This transformation weakened the intrinsic relationship between land-use and organic carbon.Yang et al. [22] concluded that band B3 had an excellent ability to predict SOC in high altitude areas (3200-3600 m).Surprisingly, our study found a lower relative importance of B3 in both periods.Although B3 was useful in distinguishing vegetation types and coverage, it showed the lowest importance.This is likely to be partly mediated by elevation and NDVI, which provide more detailed information in a landscape.
Topography has been recognized as one of the key pedogenetic factors influencing the distribution of SOC.In a high-variation of relief, topography was good predictors for mapping SOC [38,60].We found topographic variables showed higher RI in 2010 (64% in 2010, and 34% in 1990).Elevation was identified as a key topographic factor in both periods.This is consistent with previous studies [60,61].Altitude affects vegetation, local climate, and microbial activity, which have prolonged influence on SOC dynamics.TWI is related to soil moisture status, and it has been successfully used in predicting SOC in several studies [12,[62][63][64].Because of high variation in relief in the area, slope and aspect may have an effect on microclimate on local scales.We found that the influence of slope gradient was different in two periods (2% in 1990; 20% in 2010).In the study, the dataset of 1990 does not include samples in high slopes.This is because sampling is quite difficult in these areas.The purposive sampling strategy used in 2010 involving soil samples from those areas.The predictive model identified the substantial contribution of slope gradient on SOC distribution.Spatial visualization of results was essential to understand the driving forces on SOC distribution.

SOC Content Maps and Spatial-Temporal Changes
Figure 6 shows the predicted SOC in 1990 and 2010.They are means of 100 runs of BRT models.Average SOC content was higher in 1990 (20.8 ± 9.7 g•kg −1 ) than that in 2010 (18.6 ± 6.0 g•kg −1 ).The percentages of SOC content related to levels were shown in Figure 7.In both periods, SOC was high in the east and low in the southwest of Wafangdian City.High SOC content was observed in the regions covered with densely river network and forest.Spatial patterns of SOC content in 1990 were closely related to vegetation, which had been examined in recent studies [55,56,65,66].This spatial pattern was closely related to elevation in 2010.

2016, 8, 1154 10 of 16
regions covered with densely river network and forest.Spatial patterns of SOC content in 1990 were closely related to vegetation, which had been examined in recent studies [55,56,65,66].This spatial pattern was closely related to elevation in 2010.Based on cluster analysis, we obtained five grades of SOC content in two periods.In 1990, SOC content was mainly distributed at level E (>24 g•kg −1) , accounting for about 29% of the total area (Figure 7).The highest SOC content level (E) was mainly observed in the Mountain and Hilly Area of Eastern, and the lowest level (A) (<12 g•kg −1 ) was mainly associated with soils in western and southwestern parts of the coastal areas (Figure 6a).regions covered with densely river network and forest.Spatial patterns of SOC content in 1990 were closely related to vegetation, which had been examined in recent studies [55,56,65,66].This spatial pattern was closely related to elevation in 2010.Based on cluster analysis, we obtained five grades of SOC content in two periods.In 1990, SOC content was mainly distributed at level E (>24 g•kg −1) , accounting for about 29% of the total area (Figure 7).The highest SOC content level (E) was mainly observed in the Mountain and Hilly Area of Eastern, and the lowest level (A) (<12 g•kg −1 ) was mainly associated with soils in western and southwestern parts of the coastal areas (Figure 6a).Based on cluster analysis, we obtained five grades of SOC content in two periods.In 1990, SOC content was mainly distributed at level E (>24 g•kg −1) , accounting for about 29% of the total area (Figure 7).The highest SOC content level (E) was mainly observed in the Mountain and Hilly Area of Eastern, and the lowest level (A) (<12 g•kg −1 ) was mainly associated with soils in western and southwestern parts of the coastal areas (Figure 6a).
In 2010, SOC content was mainly distributed at levels A (<12 g•kg −1 ) and E (>24 g•kg −1 ), accounting for 29% and 43% of the total study area (Figure 7).High SOC content was mainly distributed in the northeast mountainous area and eastern part of Wafangdian City; the low SOC content was mainly in the western parts of the low-hilly land and coastal plains (Figure 6b).From 1990 to 2010, SOC content was improved with the shift from levels B and C to levels A and E. The area of level A and E increased by 10% and 14%, respectively.The area of level B and C decreased by 12% and 13%.
Many studies revealed that land-use change had a significant effect on SOC content distribution [12,13].Land use in Wafangdian City has undergone an intense change during the past decades, manifesting substantial consequences on SOC.Zhao et al. [10] reported an increased SOC level in all geographic areas of Jiangsu Province, China within a 26 year time period.Although there was a decreased trend in most parts of the study area, but few changes happened in the coastal plain.The area with SOC decreasing accounted for 58% of the total area mainly in the eastern mountainous (<−10 g•kg −1 ) and the central plains −5 to 0 g•kg −1 (Figure 8).The highest increase (>10 g•kg −1 ) occurred in the western and eastern plain, accounting for 9% of the total area (Figure 9).Increasing SOC (0 to 5 g•kg −1 ) content was mainly in the southwest coastal areas (28%).The highest decrease (−5 to 0 g•kg −1 ) appeared in the central part of plain and eastern mountainous regions, accounting for about 6% of the total area.2010, SOC content was mainly distributed at levels A (<12 g•kg −1 ) and E (>24 g•kg −1 ), accounting for 29% and 43% of the total study area (Figure 7).High SOC content was mainly distributed in the northeast mountainous area and eastern part of Wafangdian City; the low SOC content was mainly in the western parts of the low-hilly land and coastal plains (Figure 6b).From 1990 to 2010, SOC content was improved with the shift from levels B and C to levels A and E. The area of level A and E increased by 10% and 14%, respectively.The area of level B and C decreased by 12% and 13%.
Many studies revealed that land-use change had a significant effect on SOC content distribution [12,13].Land use in Wafangdian City has undergone an intense change during the past decades, manifesting substantial consequences on SOC.Zhao et al. [10] reported an increased SOC level in all geographic areas of Jiangsu Province, China within a 26 year time period.Although there was a decreased trend in most parts of the study area, but few changes happened in the coastal plain.The area with SOC decreasing accounted for 58% of the total area mainly in the eastern mountainous (<−10 g•kg −1 ) and the central plains −5 to 0 g•kg −1 (Figure 8).The highest increase (>10 g•kg −1 ) occurred in the western and eastern plain, accounting for 9% of the total area (Figure 9).Increasing SOC (0 to 5 g•kg −1 ) content was mainly in the southwest coastal areas (28%).The highest decrease (−5 to 0 g•kg −1 ) appeared in the central part of plain and eastern mountainous regions, accounting for about 6% of the total area.

Effects of Land-Use Change on SOC Content
Land use has significant impacts on SOC by affecting soil properties and biomass production.Dramatic land-use change, such as shifts from forest to cultivation or grassland, resulted in a completely different soil environment in terms of soil moisture and temperature regimes.It significantly influenced organic matter accumulation and decomposition dynamics [6,64].Zhao et al. [7] reported that soil organic matter (SOM) in about 26 year, when upland was changed to paddy field, had the largest increase of 5.41 ± 3.95 g k −1 in Jiangsu Province of China.
Changes in SOC from 1990 to 2010 varied with changes in land use types.A substantial change in SOC was observed when there was no land use change as well.For example, soils in unchanged cultivated land in 1990 and 2010, we observed an average loss of −2.6 g•kg −1 SOC.Soils under unchanged forest and grassland also lost amounts of SOC from 1990 to 2010.These patterns showed a decreasing SOC trend when there was no change in land use.This could be attributed to a substantial amount of production and living activities with no increase in biomass inputs to soils.
In changed land-use types, averages of SOC change ranged from −1.6 g•kg −1 to −3.6 g•kg −1 .The greatest change occurred when grassland was converted to cultivated land (G-C) and forest was converted to cultivated (F-C) (Table 5).The smallest change occurred in F-G (Forest land-Grassland).Land-use types that had been changed to grassland had considerable reduction in SOC content.For example, by shifting from forest to grassland, soil has changed from oxidation status to substitution oxidation-increase environment.It leads to a significant increase in the decomposition rate of soil organic matter (SOM) and a decline in the accumulation rate [7][8][9].

Effects of Land-Use Change on SOC Content
Land use has significant impacts on SOC by affecting soil properties and biomass production.Dramatic land-use change, such as shifts from forest to cultivation or grassland, resulted in a completely different soil environment in terms of soil moisture and temperature regimes.It significantly influenced organic matter accumulation and decomposition dynamics [6,64].Zhao et al. [7] reported that soil organic matter (SOM) in about 26 year, when upland was changed to paddy field, had the largest increase of 5.41 ± 3.95 g k −1 in Jiangsu Province of China.
Changes in SOC from 1990 to 2010 varied with changes in land use types.A substantial change in SOC was observed when there was no land use change as well.For example, soils in unchanged cultivated land in 1990 and 2010, we observed an average loss of −2.6 g•kg −1 SOC.Soils under unchanged forest and grassland also lost amounts of SOC from 1990 to 2010.These patterns showed a decreasing SOC trend when there was no change in land use.This could be attributed to a substantial amount of production and living activities with no increase in biomass inputs to soils.
In changed land-use types, averages of SOC change ranged from −1.6 g•kg −1 to −3.6 g•kg −1 .The greatest change occurred when grassland was converted to cultivated land (G-C) and forest was converted to cultivated (F-C) (Table 5).The smallest change occurred in F-G (Forest land-Grassland).Land-use types that had been changed to grassland had considerable reduction in SOC content.For example, by shifting from forest to grassland, soil has changed from oxidation status to substitution oxidation-increase environment.It leads to a significant increase in the decomposition rate of soil organic matter (SOM) and a decline in the accumulation rate [7][8][9].These results suggested that land use changes played various roles in SOC.In total, the area had experienced SOC declining from 1990 to 2010.This could be linked to the economic development of Wafangdian City during the 1980s, and the city size and population increased so rapidly that a large area of agricultural land was swallowed by construction areas.Ground vegetation was swapped and the areas were isolated from water recharge and gaseous exchange [11].Such transformation caused a tremendous decline in SOC content (Table 2).The regions with steep SOC declines were mostly situated in the eastern mountain area, where most of the forest had been converted to apple, cherry or peach plantations.These changes have negative effects on SOC accumulation (Table 2).It was unexpected that the SOC content in the southwestern hilly area, where national farms dominate and land-use types had not changed greatly, did not show a substantial change during the past 20 years.In general, there was an overall decrease in SOC (Tables 2 and 5) content but with decreasing rates from north to south of the Wafangdian City.

Conclusions
We identified the main factors driving spatial-temporal changes of topsoil organic carbon content at a regional scale in the northeast China.The spatial distributions of SOC content at two periods (1990 and 2010) were mapped by the BRT model.The results demonstrated that the BRT model is an effective and powerful tool for mapping spatial patterns of SOC content.Land-use changes caused by industrialization were the primary reasons that account for spatial-temporal changes of SOC content in the area.Variations in topography also played a key role in the spatial-temporal distribution of SOC content.Spatially, the distribution of SOC content had shown a decreasing trend from north to south between two periods.Remaining unchanged areas were located in the Low Mountain and Hilly Areas in the southwestern parts.Spatial-temporal changes of SOC content were influenced by rapid urbanization and land-use changes.Since the 1990s, land use has been severely changed such as agriculture lands converted to constructions in Wafangdian City.Because of the importance of SOC in regional carbon cycling and soil fertility, accurate spatial mapping of SOC content were essential for evaluating soil quality and carbon sequestration potential in Wafangdian City.

Figure 1 .
Figure 1.Location of sampling sites in the study area.

Figure 2 .
Figure 2. Geographic areas map of Wafangdian City.

Figure 1 .
Figure 1.Location of sampling sites in the study area.

Figure 1 .
Figure 1.Location of sampling sites in the study area.

Figure 2 .
Figure 2. Geographic areas map of Wafangdian City.Figure 2. Geographic areas map of Wafangdian City.

Figure 2 .
Figure 2. Geographic areas map of Wafangdian City.Figure 2. Geographic areas map of Wafangdian City.
where i a and i b are the predicted and observed SOC values at site i , respectively; a and b are the average values of the predicted and observed SOC values; n is the number of samples; a  and b  are the variances of predicted and observed values; and r is the Pearson correlation coefficient between the observed and predicted values.

Figure 4 .
Figure 4. Standard deviation of SOC predicted by the BRT model in the (a) 1990 and (b) 2010 surveys.BRT, boosted regression trees.

Figure 4 .
Figure 4. Standard deviation of SOC predicted by the BRT model in the (a) 1990 and (b) 2010 surveys.BRT, boosted regression trees.

Figure 5 .
Figure 5. Relative importance of each variable as determined from 100 runs of the boosted regression trees (BRT) in 1990 (a) and 2010 (b) surveys, which are shown in decreasing order and normalized to 100%.TWI, SAGA wetness index; CA, catchment area; B3, Landsat TM band 3; B4, Landsat TM band 4; B5, Landsat TM band 5; NDVI, Normalized Difference Vegetation Index.

Figure 5 .
Figure 5. Relative importance of each variable as determined from 100 runs of the boosted regression trees (BRT) in 1990 (a) and 2010 (b) surveys, which are shown in decreasing order and normalized to 100%.TWI, SAGA wetness index; CA, catchment area; B3, Landsat TM band 3; B4, Landsat TM band 4; B5, Landsat TM band 5; NDVI, Normalized Difference Vegetation Index.

Figure 6 .
Figure 6.Mean spatial distribution maps of SOC content predicted by the BRT model in (a) 1990 and (b) 2010 surveys.

Figure 6 .
Figure 6.Mean spatial distribution maps of SOC content predicted by the BRT model in (a) 1990 and (b) 2010 surveys.

Figure 6 .
Figure 6.Mean spatial distribution maps of SOC content predicted by the BRT model in (a) 1990 and (b) 2010 surveys.

Figure 8 .
Figure 8. Spatial distributions of SOC change between the1990 and 2010 surveys.Figure 8. Spatial distributions of SOC change between the1990 and 2010 surveys.

Figure 8 .
Figure 8. Spatial distributions of SOC change between the1990 and 2010 surveys.Figure 8. Spatial distributions of SOC change between the1990 and 2010 surveys.

Figure 9 .
Figure 9. Area percentages of SOC change at different levels between 1990 and 2010 surveys.

Figure 9 .
Figure 9. Area percentages of SOC change at different levels between 1990 and 2010 surveys.

Table 1 .
Summary statistics of SOC (g•kg −1 ) under different land-use patterns in two periods.
Note: SOC, soil organic carbon; SD, standard deviation; CV, coefficient of variation.

Table 2 .
SOC content in different geographic areas in 1990 and 2010 surveys.

Table 3 .
Relationships between observed SOC (g•kg −1 ) with all predictors in 1990 and 2010 surveys.

Table 4 .
Summary statistics of the predictive quality of boosted regression trees (BRT) in 1990 and 2010 for SOC prediction with 100 runs.
Note: MAE, the mean error; RMSE, root mean squared error; R 2 coefficient of determination; LCCC, Lin's concordance correlation coefficient are used to evaluate accuracy; SD, standard deviation.
[50]er R 2 , lower MAE and RMSE.Our results are comparable with previous studies.Ließ et al.[49]reported a superior prediction performance of BRT over Random forest (RF).Yang et al.[22]found a better performance of BRT while mapping SOC in the Qinghai-Tibet Plateau.Razakamanarivo et al.[50]concluded that the best model in predicting SOC stocks in the Malagasy Highlands was BRT. systematically

Table 4 .
Summary statistics of the predictive quality of boosted regression trees (BRT) in 1990 and 2010 for SOC prediction with 100 runs.
Note: MAE, the mean error; RMSE, root mean squared error; R 2 coefficient of determination; LCCC, Lin's concordance correlation coefficient are used to evaluate accuracy; SD, standard deviation.

Table 5 .
Change of SOC content under different land-use patterns during 1990-2010.

Table 5 .
Change of SOC content under different land-use patterns during 1990-2010.