Temporal and spatial changes in soil organic carbon and soil inorganic carbon stocks in the semi-arid area of northeast China

Soil organic carbon (SOC) and soil inorganic carbon (SIC) has important effects on soil physical, chemical and biological properties. They play an important role in coordinating the relationship between soil water and air, increasing soil water holding capacity and improving plant productivity. In this study, a boosted regression trees (BRT) model was developed to map the spatial distribution carbon stocks in the semi-arid region of Northeast China in 1990 and 2015. During the two periods, 10-fold cross-validation technology was used to test the performance and uncertainty of BRT model. In order to construct the model, 9 environmental variables (derived from climate, topography and biology) and 173 (1990) and 223 (2015) topsoil (0 – 30 cm) samples were used. The comparison between estimated and observed data shows that the RMSE of SOC and SIC stocks were 0.53 kgm (cid:0) 2 and 0.19 kgm (cid:0) 2 in 1990, and 0.65 kgm (cid:0) 2 and 0.20 kgm (cid:0) 2 in 2015, respectively. Elevation, normalized difference vegetation index, mean annual precipitation and Landsat band 3 were identifies as critical environmental factors for simulating the spatial distribution of SOC, accounting for 76.6 % and 70.3 % of the total relative importance in 1990 and 2015, respectively. Mean annual precipitation, mean annual temperature and topographic wetness index were the critical environmental factors for simulating the spatial variation of SIC during the two periods. Land use change also played an important role in the spatial variability of SOC and SIC stocks. In the past 25 years, soil carbon stocks decreased from 6.2 kg m (cid:0) 2 in 1990 to 5.9 kg m (cid:0) 2 in 2015. The spatial distribution pattern of SOC was high in northeastern area and low in southwestern area during the two periods, while the spatial distribution pattern of SIC was opposite to that of SOC stocks. The mapped soil carbon stock distribution is fundamental to future study of soil carbon cycle and regional carbon balance in semi-arid regions.


Introduction
Soil carbon pools include soil organic carbon (SOC) and soil inorganic carbon (SIC), which are the largest carbon reservoir in terrestrial ecosystems (Chang et al., 2012). The size of soil carbon is so large that even a small variation of its accumulation and decomposition can lead to a great fluctuation of atmospheric concentrations, directly affecting the global climate and environmental changes (Lal, 2004). Thus, soil carbon studies have attracted an extensive attention from the scientific community (Wan et al., 2019, Stockmann et al., 2015. With the significant change in global climate and environment, soil carbon pool and its changes have attracted extensive attention in the scientific community (Lal, 2004;Huang et al., 2019). It is one of the foci and hotspots of terrestrial ecosystem carbon cycle research and global change research, and also one of the core issues of a series of global change research programs such as the Global Carbon Project and the World Climate Research Programme (Schiffer and Rossow, 1983;Le Quéré et al., 2018).
Arid ecosystems occupy 47 % of the total land area (Lal, 2001), including arid, semi-arid, and semi-humid areas (Reynolds et al., 2007), accounting for about 1/3 of the global carbon sink loss (Allen et al., 1996). As an important part of global carbon cycle, carbon cycling in arid area plays an important role in regional soil carbon pools and carbon budget, which is of importance in the global biogeochemical cycle (Li et al., 2007;Lagacherie et al., 2008). Currently, the research on soil carbon cycle mainly focused more on the role and status of SOC than SIC (Kunkel et al., 2011). However, the SIC stocks are large in semi-arid areas, contributing to the regional carbon cycling (Yu et al., 2020b). Thus considering the role of both SOC and SIC in carbon cycling is important in semi-arid areas. Most SOC dynamics studies focus on SOC stocks and carbon fixation potentials in humid and semi-humid areas (Lal, 2004;Martin et al., 2011;Filippi et al., 2020) while research on SIC dynamics is relatively less in semi-arid regions (Kunkel et al., 2011). Especially in the semi-arid regions of northwest Liaoning, where the ecological environment was gradually deteriorating due to drought and little rain in spring and cold and little snow in winter. Therefore, the study of soil carbon dynamics in semi-arid areas will contribute to the understanding of carbon sequestration at the global scale.
The spatial-temporal variability of soil attributes is the result of the interaction of natural and human factors, which can be classified into five major elements: climate, terrain, parent material, time and biology (McBratney et al., 2003). Because there are many influencing factors, it is challenging to carry out spatial prediction of soil properties on the regional scale. Digital Soil Mapping (DSM) technology provides a fast and economical method for spatial prediction of SOC and SIC in largescale areas based on a small number of sampling point data and major environmental covariates (Martin et al., 2011;McBratney et al., 2003;Yang et al., 2020b). In DSM methods, the spatial-temporal change of soil carbon stocks was estimated by using multi-period soil carbon sample data. In order to explore the relationship between SOC and SIC values and environmental factors (such as annual average temperature, annual precipitation, elevation, etc.), to find powerful prediction factors or indicators related to the distribution of SOC and SIC values, and to draw the spatial change map of SOC and SIC according to the field observation results, some researches have been carried out. In general, the relationship between soil attributes and environmental factors is nonlinear and complex, so the machine learning algorithm which can effectively solve these problems is introduced into the spatial prediction of soil attributes (Martin et al., 2011;Song et al., 2016;Wang et al., 2018;Yang et al., 2020c).
Among different DSM technologies, tree-based models are widely used to predict soil properties, such as soil salt, pH, SOC, and texture (Martin et al., 2011;Wang et al., 2017;Wang et al., 2018). Compared with traditional methods, tree-based models have better performance and effectiveness (Carslaw and Taylor, 2009;Wang et al., 2019;Yang et al., 2020a). The most reliable prediction model is included in DSM toolbox. The boosted regression trees (BRT) model can improve the performance of the model by training multiple regression tree models and combining them to predict, which can effectively avoid the problem of transition fitting, effectively deal with nonlinear and complex problems, and solve qualitative and quantitative problems (Cheong et al., 2014;Yang et al., 2016a,b). Thus BRT model is widely used in remote sensing science, epidemiology, ecology, and fishery science (Müller et al., 2013;Cheong et al., 2014;Wang et al., 2016). However, there are few studies on the spatiotemporal changes of SOC and SIC with the method.
The purpose of this study is to apply the relevant principles of quantitative pedology to analyze the spatiotemporal changes in SOC and SIC stocks in northwest Liaoning. The specific research objective are to: (1) establish a mathematical model to predict SOC and SIC stocks based on environmental impact factors and sampling data in 1990 and 2015; (2) quantify the key environmental variables of SOC and SIC stock changes; (3) reveal the spatial distribution characteristics of soil carbon stocks under different land use patterns in 1990 and 2015; and (4) analyze the spatiotemporal change rules of soil carbon stocks in the region during the past 25 years .

Study area
The study area is located in the semi-arid area of northwest Liaoning in northeast China (latitude 39.98 • -43.48 • N; longitude 118.83 • -124.43 • E), referring to Chaoyang, Tieling and Fuxin. Northwest Liaoning refers to Chaoyang, which is located in the long and narrow area on the south edge of Horqin sandy land, the largest in China. The terrain is dominated by hills and mountains, and the terrain is stepped down from southwest to northeast. This area belongs to the temperate monsoon climate area of semi humid to semi-arid transition. It is dry and windy in spring, hot and rainy in summer, cool and dry in autumn and cold and dry in winter. The mean annual temperature (MAT) is 6.4-8.5 ℃, the annual sunshine hours are 2850-2960 h, the mean annual precipitation (MAP) is 300-700 mm, the annual average evaporation and gale days are 1000-1500 mm and 22-32 d, respectively. The annual accumulated temperature is higher and the precipitation is less. It is known as the "nine droughts in ten years" and is a heavy drought area in Liaoning Province. There are three major water systems in the area, namely Liaohe River, Daling River and Raoyang River. According to the Chinese Soil Taxonomy (Gong, 2002), the dominant soil types are Aridosols (48 %) and Primosols (23 %), and the soil layer is relatively thin (mostly 10-30 cm). The vegetation has the characteristics of both North China flora and Inner Mongolia flora temperate steppe Quercus forest. The area is primarily dominated with sparsely-distributed shrubs and herbs.

Soil survey data for 1990
In 1990, the soil survey data was obtained from Agricultural Technology Extension Centers in Chaoyang, Fuxin and Tieling. The database comprised 171 soil profile information databases (Fig. 1), including soil physical and chemical properties, land-use types, soil types, parent material information, slope gradient, slope aspect, mean annual precipitation (MAT), mean annual temperature (MAT) and other information. In this database, the SOC and SIC contents were measured using the Walkley black wet combination method (Nelson and Sommers, 1982) and a gas volume method (Raskin, 1983), respectively. To estimate dry bulk density, 100 cm3 of soil cores dried for 48 h at 105 • C for bulk density measurement. The soil sample sits were mainly covered with cultivated land (n = 73), forest (n = 52) and grassland (n = 46). This study focused on the SOC and SIC of topsoil (0-30 cm) in the dataset. We used Pedo-Transfer Functions (PTFs) to fill some missing bulk density values in the database with formula: where BD represents bulk density; SOC represents SOC content in topsoil layer (0-30 cm).

Soil sampling in 2015
Due to the large East-West span of the study area, and the terrain is mainly mountainous and hilly areas, considering the time and economic cost, it is difficult to collect the samples in large quantities. In order to accurately obtain the spatial variation of SOC and SIC in the complex geographical landscape units in the region, a purpose sampling method (Yang et al., 2019;Zhu et al., 2008) was selected to design the sampling scheme. Firstly, we selected the environmental variables of MAP, MAT, elevation, slope gradient, normalized difference vegetation index, and topographic wetness index, and clustered them by using the fuzzy cmeans clustering method. As a result, we obtained 16 landscape units. Secondly, depending on the land use patterns, soil type and road accessibility, for each landscape unit, 15-20 sampling points are collected, and 223 sampling points are finally obtained (Fig. 1). Finally, the number of sampling sites for cultivated land, grassland and forest was 104, 68 and 51, respectively. A hand-held GPS was used to locate each sampling geographic coordinate. The sampling depth of each point is 30 cm, and 100 cm 3 of undisturbed soil cores were collected at the center (15 cm) for the laboratory to determine the soil bulk density. After that, about 1 kg samples were collected from each site at the same location for later SOC and SIC measurements in the Central Laboratory of the college of land and environment, Shenyang Agricultural University, China. The determination methods of SOC, SIC and dry bulk density were the same as those in 1990.

Environmental variables
Nine environmental variables including MAP, MAT, elevation (ELE), slope gradient (SG), topographic wetness index (TWI), Landsat band 3 (B3), Landsat band 4 (B4), Landsat band 5(B5), and normalized difference vegetation index (NDVI), which represent climate, topography and biology among the five soil forming factors, are used to predict SOC and SIC reserves in the surface layer of Northwest Liaoning Province. Because these variables are obtained from different departments and platforms, we used Arcgis10.2 (ESRI Inc., USA) software to resample them to 90 m × 90 m spatial resolution, and unify the projection coordinates of environmental variables into WGS_1984_UTM_Zone_50N coordinate system for later modeling and analysis.
Climatic datasets of 1990 and 2015 were obtained from China Meteorological Data Service Center (https://data.cma.cn/en), composed of MAT and MAP. These climatic data were based on 1 km grid data generated by Kriging interpolation from 673 weather stations across China. A 90 m spatial resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) data was downloaded from the United States Geological Survey (USGS). The elevation gradient of the study area ranged from 1 m to 1235 m. The high elevation was mainly distributed in the southwest low mountains and hills, and the low altitude was mainly distributed in the northwest. Two topographic variables (ELE, SG) were derived from SRTM DEM in ArcGIS 10.2 using a spatial analysis model. The System for Automated Geoscientific Analyses Geographic Information System (SAGA GIS) software was used to calculate TWI by SRTM DEM. B3, B4 and B5 representing vegetation formation, coverage and biomass (Wang et al., 2018) were obtained from Landsat TM and ETM + dataset and downloaded from the USGS with a 30 m × 30 m spatial resolution, covering the growing season from July to September in 1990 and 2015. NDVI was derived from B3 and B4 to detect the vegetation growth status and vegetation coverage:

Boosted regression trees and their uncertainty
The BRT model was proposed by Friedman et al. (2000), which consists of boosting and regression trees (Elith et al., 2008). Boosting technology was an improvement based on random gradient of decision tree (Wang et al., 2016). It used all samples at once and changes the weight of the samples in each round of training. The goal of the next round of training was to find a function to adapt to the residual error of the previous round. When the residual was small enough or reached the maximum number of iterations, it stopped iterating. BRT model could deal with soil environmental problems in complex landscape area flexibly and avoid nonlinear interactions between variables. In comparison with the traditional regression model, BRT model had better prediction performance, especially in the spatial simulation of soil properties, and have been widely used in spatial prediction research. In this study, we use the "demo" software package developed by Elith et al. (2008) and build the model in R language environment (R Development Core Team, 2013). To evaluate the prediction uncertainty of SOC and SIC stocks during the two periods, we iterated 100 times of the BRT model and obtained an average standard deviation map.

Model validation
A 10-fold cross-validation method with four common validation indices including mean absolute error (MAE), root mean square error (RMSE), coefficient of determination (R 2 ), and Lin's concordance correlation coefficient (LCCC) (Lin, 1989) was used to evaluate the performance of the BRT model in the two periods. Among them, MAE was used to evaluate the deviation degree of prediction value. The closer to 0, the better the prediction result. RMSE was used to evaluate the overall accuracy of the prediction results. The smaller the value was, the higher the prediction accuracy of the model was. R 2 was used to evaluate the goodness of fit of the model. The closer the value was to 1, the higher the reference value of the model was. LCCC was used to measure the degree of 1-to-1 line distribution between the predicted value and the measured value. The closer the LCCC was to 1, the higher the degree of agreement between the predicted value and the observed value was, and the stronger the prediction ability of the model was. The specific calculation formula are: and predicted values; n is the number of sampling points; r is the Pearson correlation coefficient.

Calculation of SOC and SIC stocks
This study focused on the temporal and spatial variation of soil carbon stocks in the semi-arid area of northeast Liaoning Province. Therefore, it was necessary to accurately calculate the SOC and SIC stocks of the topsoil (0-30 cm) at the sampling location. The calculation formulas are as follows: where SOC density , SIC density , SOC content , and SIC content are the SOC density, SIC density, SOC content, and SIC content, respectively; BD j , D j, and S j represent the bulk density, the soil layer thickness, and the coarse fraction content greater than 2 mm, and the coarse fraction content was calculated by a weighting method; k is the layer of soil profile; j is the specific level, this study is limited to 30 cm. To further reveal the causes of these changes in this study region, we calculated the changes in soil carbon stocks under the main land use types including cultivated land, forest land, and grasslands. In order to further explore the variable trend of soil carbon stocks in Northwest Liaoning Province in the past 25 years, we calculated the spatial variation maps of soil carbon stocks in the two periods.

Descriptive statistics
The descriptive statistical results of SOC and SIC stocks with selected environmental variables at all sampling sites were shown in Figs. 2 and 3. In 1990, the SOC and SIC stocks were 0.38-14.82 kg m − 2 and 0.18-3.13 kg m − 2 , with the average values were 4.79 kg m − 2 and 1.33 kg m − 2 , respectively. The corresponding average values of SOC and SIC in 2015 were 4.69 kg m − 2 and 1.64 kg m − 2 , respectively. In addition, we also calculated the skew coefficients of SOC and SIC stocks, which were 1.14 vs 1.08 in 1990 and 0.78 vs 1.46 in 2015, respectively. Therefore, the data did not follow the standard normal distribution. So in the process of building model, it is necessary to carry out logarithmic transformation to make it conform to the normal distribution. Relationships between logarithmic conversion of SOC stocks (kg m − 2 ) and SIC stocks (kg m − 2 ) with all predictors were shown in Table 1. During the two periods, MAP, NDVI, ELE and SG were positively correlated with the logarithmic conversion of SOC stocks, while MAT, TWI, B3, and B4 were negatively correlated. Correspondingly, MAT, TWI, and B4 were positively correlated with the logarithmic conversion of SIC stocks, while MAP, ELE, SG, and NDVI were negatively correlated.

Model performance and uncertainty
Twenty percent of all sampling points were randomly selected as independent testing dataset to evaluate the prediction performance of BRT model during the two periods (Table 2). In order to ensure the stability of model performance, BRT model was iterated 100 times and its average value was calculated as the final prediction result. The results of model validation showed that the BRT model was efficient and robust in predicting the spatial distribution of SOC and SIC stocks in both periods, because the BRT model had higher R 2 and LCCC, lower MAE and RMSE. In addition, in order to evaluate the uncertainty of BRT model, we calculated the average standard deviations (SDs) of prediction results using BRT model with 100 iterations. In 1990, the SDs of SOC and SIC stocks were 0.047 kg m − 2 and 0.046 kg m − 2 , respectively. In 2015, the SDs of SOC and SIC stocks were 0.049 kg m − 2 and 0.066 kg m − 2 , respectively. The BRT model presented low uncertainty in both periods (Fig. 4).

Relative importance of environmental variables
By iterating BRT model for 100 times, the average relative importance (RI) of 9 environment variables was calculated and normalized to 100 %. We found that each environmental variable had different RI in SOC and SIC prediction (Fig. 5). ELE, NDVI, MAP B3, and MAT were the key environmental factors affecting SOC stocks, accounting for 76.6 % and 70.3 % of the total RI in 1990 and 2015, respectively. Correspondingly, MAP, MAT, elevation and TWI were the key environmental factors affecting the spatial distribution of SIC stocks, accounting for 68.1 % and 69.7 % of the RI in 1990 and 2015, respectively.

Spatial prediction of soil carbon stocks
In this study, 100 iterations of BRT were selected for spatial prediction of SOC and SIC stocks, and the average value of their results was calculated as the final spatial distribution map of SOC and SIC stocks (Fig. 6). SOC stocks in 2015 (4.2 ± 1.4 kg m − 2 ) were lower than those in 1990 (4.6 ± 1.5 kg m − 2 ). Similarly, SIC stocks in 2015 (1.8 ± 0.6 kg m − 2 ) were lower than those in 1990 (1.5 ± 0.5 kg m − 2 ). In order to better xplore the spatial and temporal changes of soil carbon stocks in this region, we added SOC and SIC maps during the two periods in ArcGIS 10.2 using grid calculation module to obtain the spatial distribution map of soil carbon stocks in 1990 and 2015 (Fig. 7), respectively. In addition, soil carbon stocks decreased from 6.2 ± 1.9 kg m − 2 in 1990 to 5.9 ± 1.9 kg m − 2 in 2015. According to the spatial distribution maps  of soil carbon stocks in the two periods, the soil carbon stocks gradually decreased from northeast to southwest.
The spatial variation maps of soil carbon stocks in the two periods was calculated (Fig. 8), the results showed that soil carbon stocks decreased mainly in the northeast and central regions, accounting for 70 % of the total area of the study area, and the decline range was between − 2-1 kg m − 2 and − 1-0 kg m − 2 , accounting for 66 % of the total area (Fig. 8b). Soil carbon stocks increased mainly in the southwest of Northwest Liaoning Province, accounting for 30 % of the total area of the study area.

Changes of soil carbon stocks under different land use patterns
The changes in soil carbon stocks under the main land use types was calculated (Table 3), we found that the decrease of soil carbon stocks was mainly due to the cultivated lands with no land-use changes, and the decrease was 5.54 Tg C. The increase mainly came from forest and grassland with no land use change, increasing by 1.44 and 1.42 Tg C, respectively.

Effects of environmental factors on SOC and SIC stocks
Climatic factors had important impacts on the spatial distribution of SOC and SIC stocks (Wynn et al., 2006;Wang et al., 2018). Previous studies (Batjes, 1996;Papatheodorou et al.,2004;Na et al., 2010;Filippi et al., 2020) have shown that SOC and SIC were significantly correlated with climatic variables, which is further verified in this study (Fig. 5). In the processes of carbon input and output, vegetation productivity and microbial decomposition and transformation are affected by climate conditions such as rainfall and temperature . SOC stocks were strongly affected by rainfall and temperature. SOC content in natural ecosystem decreased exponentially with the increase of temperature (Willaarts et al., 2016). Compared with the arid area, the semiarid area had more MAP, larger cultivated land area, and the soil erosion intensity increases correspondingly (Na et al., 2010). The organic matter attached to the light clay particles of the topsoil was also lost with the soil, which reduced the accumulation of soil organic matter (Wynn et al., 2006). However, the vegetation coverage in the arid area was low, but a large number of biological crusts grow on the soil surface, which not only makes the surface soil tend to be fixed, but also strengthens the soil (Yang et al., 2016;Wang et al., 2018;Filippi et al., 2020). Drought environment is not conducive to plant growth; high temperature can promote microbial activity and accelerate organic matter mineralization, but is not conducive to the synthesis and accumulation of soil humus, resulting in low SOC stocks (Filippi et al., 2020).
MAP was an important component in the process of calcium carbonate deposition (CaCO 3 ) (Bolinder et al., 2007). In semi-arid area, appropriate soil moisture content was conducive to the formation of CaCO 3 . Lal (2004) conducted that vegetation and microbial activities could significantly change the infiltration of water, thus changing the leaching of SIC and precipitation of secondary carbonate. The profile distribution of SIC affected by leaching is also considered to be more complex. Even under the combined action of heavy rainfall and violent biological activities, carbonate is almost completely leached (Na et al.,

Table 1
Relationships between logarithmic conversion of SOC stocks (kg m − 2 ) and SIC stocks (kg m − 2 ) with all predictors in the 1990 and 2015 surveys.  2008). However, the leaching loss of natural soil base was less in semiarid or arid areas. On the contrary, due to the large evaporation of soil water, with the increase of capillary water, the base materials in the lower layer tend to accumulate to the upper layer, resulting in calcium reaction in the soil (Papatheodorou et al., 2004). Therefore, climatic factors such as MAT and MAP had obvious effects on the distribution of SOC and SIC stocks, which were important covariates for spatial distribution prediction. Topography was an important factor in the process of soil formation (Wang et al., 2016). It not only controls the redistribution of water and heat resources, but also affects the material circulation process and intensity of soil ecosystem, and has a profound impact on soil properties (Martin et al., 2011;Yang et al., 2016;Wang et al., 2018). Topography is an important factor to affect ecosystem material and energy flows (McBratney et al., 2003). As the most important topography variable, ELE does not directly affect SOC and SIC stocks, but changes the distribution of soil organic matter through altering bioclimatic factors such as rainfall and temperature (Yang et al., 2016). In addition, in the southwest of Liaoning Province, the mountainous and hilly areas are the main ones. The terrain factors such as SG and slope position were easy to cause geological subsidence and soil erosion, and a large amount of soil nutrients are lost, which was not conducive to the accumulation of organic matter, resulting in the minimum SOC stocks in this region.
The biological variables were the most important factor affecting the spatial variation of SOC stocks in both periods (Fig. 5). This conclusion is consistent with previous findings (Wynn et al., 2006;Willaarts et al., 2016;Yang et al., 2020). Wang et al. (2017) found that biological related variables were closely related to the spatial distribution of SOC stocks, and they were the main source of topsoil SOC stocks and controlled the amount of organic matter entering the soil. In this study, we also found that NDVI and B3 are the most important environmental factors among all biological variables in predicting SOC stocks, and they represent the biomass and productivity of vegetation to a certain extent. In the Qilian Mountains of China, Yang et al (2016) used BRT and RF models to predict the topsoil SOC in an alpine ecosystem, and concluded that NDVI and B3 are the powerful environmental factors among the four biomass variables (NDVI, B3, B4, and B5). Surprisingly, biological variables showed the lowed relative importance in the spatial simulation of SIC stocks. Through analysis, we believed that the source and distribution of SIC stocks were often influenced by multiple factors such as pedogenic parent material, soil moisture, climate, salinity and soil type. Similar conclusions have been reached in previous studies Wang et al., 2010, Rong et al., 2012.

Effects of land use change on soil carbon stocks
Land use is the most direct result of human activities (Adhikari et al., 2019). Under different land use patterns, the disturbance degree of surface vegetation and soil is different, resulting in significant differences in soil carbon stocks (Wang et al.,2016;Li et al., 2010;Yu et al., 2020a). Our study further clarified the effects of land-use change on soil carbon stocks. We found that forest and grassland without land-use change increased carbon stocks by 1.44 and 1.42 Tg C in the past 25 years, respectively, while cultivated land carbon decreased by 5.54 Tg C (Table 3). Changes in land-use patterns will lead to changes in vegetation types, soil microorganisms, and soil physical and soil chemical properties, which will affect the amount and decomposition rate of soil organic matter, and then affect SOC stocks at different degrees (Song et al., 2016;Yigini and Panagos, 2016;Wang et al., 2017). SIC mainly refers to carbonate existing in arid and semi-arid soils, and its content is related to the lithology of parent material (Chang et al., 2012;Du and Gao, 2020;Yu et al., 2020a). Land use change will affect the leaching and deposition process of SIC, and then affect the SIC stocks Du and Gao, 2020). In the Wisconsin of USA, Huang et al. (2019) applied a modified space-for-time substitution method to estimate the spatial and temporal variations of SOC stocks over a 150-year period. They concluded that the spatial variation of SOC stocks at 0-30 cm was mainly affected by land cover and soil types with the largest SOC stocks found in forest and wetland and spodosols. Wu et al. (2009) used soil profile data (China's second national soil survey) to investigate the spatial distribution of SIC for the entire country under present day conditions as well as changes in SIC under historical land use conditions. They also found that human activity may have had a great impact on SIC stocks.
The regional soil carbon stocks were decreased from 1990 to 2015 (Table 3) due to the main land use type as cultivated land in this area, which was China's commodity grain base (Wang et al., 2017). Soil organic matter formed by the residual organic matter of cultivated crops in soils was not lower than the organic matter consumed due to mineralization (Nyssen et al., 2008;Chang et al., 2012). This leads to the decrease of soil fertility every year, especially in recent years with the gradual decrease of cultivated land area and the continuous increase of population (Na et al., 2010;Wang et al., 2016). In the southeast of the region, the trend of increasing soil carbon stocks was mainly due to the implementation of the policy of returning farmland to forest and protecting natural forest (Fig. 7). Over the past 25 years, 2262 km 2 of cultivated had been converted to forest, and soil carbon stocks had decreased by 0.22 Tg C (Table 3). Since the implementation of this policy, the accumulation of litter and the increase of underground root biomass have led to more soil carbon stocks. In addition, forest increased the number and activity of soil microorganisms and animals, accelerated the turnover of organic matter, deepened the rooting system of trees, and was conducive to the accumulation of soil carbon stocks. However, as this region is the main commodity grain base in China, the implementation time of the conversion of farmland to forests in this region varies. In relatively large regions, this policy was implemented only in recent 5-10 years (Wang et al., 2018). Overall, there was not an increasing trend in soil carbon storage changes. In a lake retreat area of the Ethiopian Rift Valley, Nyssen et al. (2008) concluded that land use and cover changes lead to loss of vegetation cover and subsequent change in SOC and soil quality. In general, soil carbon stocks increased after returning cultivated land to forest and grassland (Table 3), but the change to cultivated land will reduce soil carbon stocks.

Relationships between SOC and SIC stocks
We found that SOC and SIC stocks have a significant negative correlation in the semi-arid (Table 1). They showed opposite correlations with the selected environmental variables in the two periods (Fig. 6). This finding was confirmed in previous studies ;  Zhang et al., 2013;Du and Cao, 2020). By comparing the spatial distribution of SOC and SIC stocks in the two periods, we find that there was an opposite and different trend of change in the topsoil, while the southwest of the study area had higher SOC stocks but lower SIC stocks. At the edge of Gurbantunggut Desert in Northwest China, Rong et al. (2012) found that the SOC increased from the topsoil to the bottom, but there was an opposite change trend of SIC. Chang et al. (2012) believed that afforestation affected SOC and SIC, and afforestation within 20 cm of the profile increased SOC and decreased SIC. However, some studies showed that there was a positive correlation between SOC and SIC stocks (Chadwick et al., 1994;Li et al., 2010;Cao, 2012). Li et al. (2010) found that, compared with fallow soil, long-term continuous cropping or rotation of different crops significantly increased the mass fraction of various forms of carbon in 0-40 cm soil layer, increased SOC by 47-139 %, and increased SIC by 20-26 %. In the Tengger Desert of China, Zhang et al. (2009) concluded that the distribution pattern of organic matter in light brown calcareous soil was similar to that of calcium carbonate in natural environment, and there was a significant positive correlation between SOC and SIC, and the change trend was the same. Therefore, the SOC and SIC dynamics were complex. In the future research on soil carbon stocks in arid and semi-arid areas, both SOC and SIC dynamics and their influence factors should be considered.

Limitations in the current study
BRT model well predicted SOC and SIC stocks in the semi-arid area of northwest Liaoning in northeast China. However, there were some uncertainties in this study. First, some soil profile data in 1990 were missing bulk density data, which was supplemented with Pedo-Transfer Functions (PTFs). This may lead to prediction error. Second, the soil data in 2015 were collected and analyzed by different groups, which may lead to sampling and measurement errors. Third, environmental data were obtained from different departments and platforms, and the data accuracy was different. Therefore, some data may be lost in the process of resampling. Fourth, the content and distribution of SIC are affected by multiple factors such as soil parent material, soil moisture, climate, salinity and soil type. This kind of environmental data should be fully introduced into the spatial prediction of SIC stocks. But due to the difficulty of data acquisition, this study only obtained climatic data, which may cause prediction errors in predicting SIC stocks. Finally, this study was only limited to the study of topsoil carbon stocks (0-30 cm) in Northwest Liaoning Province. This will underestimate the total soil carbon stocks, because usually more soil carbon stocks are in deep soils instead of topsoil.

Conclusion
BRT model was used to predict the spatial distribution of SOC and SIC stocks in the semi-arid region of Northeast China. The main environmental factors affecting the spatial and temporal changes of SOC and SIC stocks were identified and the temporal and spatial changes of soil carbon stocks were estimated. We found that BRT model was able to predict the spatial distribution of SOC and SIC stocks with high R 2 and LCCC and low MAE and RMSE in comparison with the observation data. In addition, SOC and SIC stocks showed opposite spatial distribution characteristics in the two periods. During the 25 years, SOC stocks decreased from 4.9 kgm − 2 in 1990 to 3.9 kgm − 2 in 2015, corresponding SIC stocks increased from 2.9 kg m − 2 in 1990 to 4.2 kg m − 2 in 2015. However, soil carbon stocks decreased by 1 kg m − 2 during the 25 years, the spatial distribution pattern gradually increased from southwest to northeast. In addition, we found ELE, NDVI, MAP and B3 were the key environmental factors for simulating SOC stocks, correspondingly, MAP,  MAT, elevation and TWI were the key environmental factors for simulating the spatial distribution of SIC stocks during the two periods. We also found that land use change played an important role in the spatial variability of SOC and SIC stocks. In general, the prediction accuracy was high and reasonable. The predicted temporal and spatial distribution of soil carbon stocks is valuable information for soil and water conservation, environmental protection and agricultural production planning in Northwest Liaoning Province.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The authors do not have permission to share data.