Simulated Soil Organic Carbon Density Changes from 1980 to 2016 in Shandong Province Dry Farmlands Using the CENTURY Model

The changes in cultivated soil organic carbon (SOC) have significant effects on soil fertility and atmospheric carbon dioxide (CO2) concentration. Shandong Province is an important agricultural and grain production area in China. Dry farmland accounts for 74.15% of the province’s area, so studies on dynamic SOC changes would be helpful to understand its contribution to the Chinese national carbon (C) inventory. Using the spatial overlay analysis of the soil layer (1:10,000,000) and the land use layer (1:10,000,000), 2329 dry farmland soil polygons were obtained to drive the CENTURY model to simulate SOC dynamics in Shandong Province from the period 1980 to 2016. The results showed that the CENTURY model can be used to simulate the dry farmland SOC in Shandong Province. From the period 1980 to 2016, the soil organic carbon storage (SOCS) and soil organic carbon density (SOCD) showed an initial increase and then decreased, especially after reaching a maximum in 2009. In 2016, the SOCS was 290.58 × 106 t, an increase of 26.99 × 106 t compared with 1980. SOCD in the dry farmland increased from 23.69 t C ha−1 in 1980 to 25.94 t C ha−1 in 2016. The dry farmland of Shandong Province was a C sink from 1980 to 2016. Among the four soil orders, inceptisols SOCD dominated, and accounted for 47.81% of the dry farmland, followed by >entisols > vertisols > alfisols. Entisols SOCD growth rate was the highest (0.23 t C ha−1year−1). Compared to 1980, SOCD in 2016 showed an increasing trend in the northeast, northwest and southeast regions, while it followed a downward trend in the southwest.


Introduction
According to the fifth report of the Intergovernmental Panel on Climate Change (IPCC), climate warming is evident [1]. The causes of global warming are diverse, such as the burning of large amounts of fossil fuels and industrial production. With increased global warming, a series of related problems have arisen, such as the accelerated melting of polar glaciers, rising sea levels, reduced global freshwater resources, and reduced crop production [2,3]. Therefore, how to mitigate global warming has become an important issue in related fields. The rational use and management of agricultural soil is an effective way to alleviate climate warming [4,5]. It is believed that a large amount of C is stored in the soil, and its activity may affect atmospheric CO 2 concentration. However, although soil is an important C source for atmospheric CO 2 , historical studies show that land use reduces SOC by 4.1 × 10 10 t-5.5 × 10 10 t [6], and global soil erosion has resulted in an SOC loss of × 10 5 km 2 , running about 420 km from north to south and, at its widest, is about 700 km. In 2017, the total population of the province was 10.06 million. Shandong Province has a warm temperate monsoon climate with an annual average temperature of 11-14 °C. The annual frost-free period increases from the northeast coast to the southwest, sunshine hours are adequate, with an average of 2290-2890 h per year, temperature conditions meet the needs for two crops a year, and average annual precipitation is generally between 550 and 950 mm.

Study Method
The CENTURY 4.5 model was used to simulate SOCD and SOCS changes in dry farmland from 1980 to 2016. This model was originally developed by Parton et al. (1987) [35], and was used to simulate the dynamic processes of C, nitrogen (N), phosphorus (P), sulfur (S) and other elements in grassland ecosystems. Later, after a series of parameter modifications, it was gradually applied to forest, crop, and savanna ecosystems. The primary uses of the model are as a tool for ecosystem analysis, as a test for data consistency, and as an evaluation tool to study the effects of changes in management and climate on ecosystems.
In the study area, we selected dry farmland soil for our research. Cultivated in Shandong Province, it covers 7.07 × 10 6 ha, of which 97% is dry farmland [36]. Initially, we extracted 2329 dry farmland soil polygons using spatial overlay analysis of soil and land use layers (1:10,000,000). Additionally, we used data from 186 soil profiles obtained from the second national soil survey data of Shandong Province. The soil profile data included important soil physical and chemical properties, such as soil clay content, SOM, bulk density (BD), and pH. Furthermore, we obtained data from the Chinese National Meteorological Center for 31 climate stations in Shandong Province, including site monthly precipitation, and highest/lowest average temperature. Finally, we acquired data from the Resources and Environmental Scientific Data Center, Chinese Academy of Sciences (RESDC-CAS) on fertilizer and organic fertilizer use, and crop growth status data of 70 sites in Shandong Province, including crop type and site planting area.
Based on the CENTURY 4.5 model, SOC evolution, and soil profile points were simulated. In ArcGIS, soil profile points were correlated with dry farmland soil polygons, and the evolution of SOCD and SOCS in Shandong Province from 1980 to 2016 were obtained.

Study Method
The CENTURY 4.5 model was used to simulate SOCD and SOCS changes in dry farmland from 1980 to 2016. This model was originally developed by Parton et al. (1987) [35], and was used to simulate the dynamic processes of C, nitrogen (N), phosphorus (P), sulfur (S) and other elements in grassland ecosystems. Later, after a series of parameter modifications, it was gradually applied to forest, crop, and savanna ecosystems. The primary uses of the model are as a tool for ecosystem analysis, as a test for data consistency, and as an evaluation tool to study the effects of changes in management and climate on ecosystems.
In the study area, we selected dry farmland soil for our research. Cultivated in Shandong Province, it covers 7.07 × 10 6 ha, of which 97% is dry farmland [36]. Initially, we extracted 2329 dry farmland soil polygons using spatial overlay analysis of soil and land use layers (1:10,000,000). Additionally, we used data from 186 soil profiles obtained from the second national soil survey data of Shandong Province. The soil profile data included important soil physical and chemical properties, such as soil clay content, SOM, bulk density (BD), and pH. Furthermore, we obtained data from the Chinese National Meteorological Center for 31 climate stations in Shandong Province, including site monthly precipitation, and highest/lowest average temperature. Finally, we acquired data from the Resources and Environmental Scientific Data Center, Chinese Academy of Sciences (RESDC-CAS) on fertilizer and organic fertilizer use, and crop growth status data of 70 sites in Shandong Province, including crop type and site planting area.
Based on the CENTURY 4.5 model, SOC evolution, and soil profile points were simulated. In ArcGIS, soil profile points were correlated with dry farmland soil polygons, and the evolution of SOCD and SOCS in Shandong Province from 1980 to 2016 were obtained.

Soil Data
Shandong Province dry farmland soils include mainly inceptisols, entisols, alfisols, and vertisols, which account for 47. 16, 10.32, 35.02, and 6.95% of the dry farmland soil polygons, respectively. Figure 2 shows the distribution of each soil order. Table 1 includes the physical and chemical properties of the four main soil orders in Shandong Province dry farmland soils [37,38].
Sustainability 2020, 12, x FOR PEER REVIEW 4 of 17 2 shows the distribution of each soil order. Table 1 includes the physical and chemical properties of the four main soil orders in Shandong Province dry farmland soils [37,38].

Model Initialization and Parameterization
Model simulation includes four main processes: model parameterization, model initialization, model simulation and result output. The CENTURY model obtains input values through 12 "100" files, each of which contains a certain subset of variables. These 12 "100" files include crop.100, cult.100, fert.100, fire.100, fix.100, graz.100, harv.100, irri.100, omad.100, tree.100, trem.100, <site>.100. Corresponding to regional crop options, farming options, fertilizer usage options, files with fixed parameters mainly related to organic matter, grazing options, fire options, irrigation conditions, organic matter addition options, tree type options, tree removal options and site-specific parameters. Users can create new options or modify the values in the options through the FILE100 program to meet the parameter settings of different research areas [39]. For example, in the fertilizer addition

Model Initialization and Parameterization
Model simulation includes four main processes: model parameterization, model initialization, model simulation and result output. The CENTURY model obtains input values through 12 "100" files, each of which contains a certain subset of variables. These 12 "100" files include crop.100, cult.100, fert.100, fire.100, fix.100, graz.100, harv.100, irri.100, omad.100, tree.100, trem.100, <site>.100. Corresponding to regional crop options, farming options, fertilizer usage options, files with fixed parameters mainly related to organic matter, grazing options, fire options, irrigation conditions, organic matter addition options, tree type options, tree removal options and site-specific parameters. Users can create new options or modify the values in the options through the FILE100 program to meet Sustainability 2020, 12, 5384 5 of 17 the parameter settings of different research areas [39]. For example, in the fertilizer addition project, the first input is the type of fertilizer (nitrogen fertilizer, phosphate fertilizer, sulfur fertilizer), and the amount of fertilizer (g/m 2 ) input into the soil.
The initialization process of the model is to simulate the state of the soil from the original state to the time when the simulation starts, and prepare for the simulation of the research period. The EVENT 100 program was used to set the event, and it was assumed that the SOC accumulation was only affected by environmental factors, which made the model run for thousands of years. In this study, it was set to 5000 years. The temperate C3 grassland vegetation type and short-term appropriate grazing and combustion management system The C pool reaches equilibrium [40]. From 1000 BC to 1954, the land began to use a small amount of organic fertilizer. From 1955 to 1985, a small amount of chemical fertilizers and organic fertilizers were used in agricultural production, and by adjusting the parameters so that the simulated and measured values until 1980 were similar, it was regarded as completing the initialization process of the CENTURY model. The model simulation flow chart is depicted in Figure 3. project, the first input is the type of fertilizer (nitrogen fertilizer, phosphate fertilizer, sulfur fertilizer), and the amount of fertilizer (g/m 2 ) input into the soil. The initialization process of the model is to simulate the state of the soil from the original state to the time when the simulation starts, and prepare for the simulation of the research period. The EVENT 100 program was used to set the event, and it was assumed that the SOC accumulation was only affected by environmental factors, which made the model run for thousands of years. In this study, it was set to 5000 years. The temperate C3 grassland vegetation type and short-term appropriate grazing and combustion management system The C pool reaches equilibrium [40]. From 1000 BC to 1954, the land began to use a small amount of organic fertilizer. From 1955 to 1985, a small amount of chemical fertilizers and organic fertilizers were used in agricultural production, and by adjusting the parameters so that the simulated and measured values until 1980 were similar, it was regarded as completing the initialization process of the CENTURY model. The model simulation flow chart is depicted in Figure 3.

SOCD Calculation
SOCD was selected as representative of the spatial and temporal distribution characteristics of SOC [41,42]. The SOCD is the SOCS per unit area where SOCi (g/kg) is the organic carbon content of the ith soil polygons, SOMi (g/kg) is the organic matter content of the ith soil polygons, 0.58 is the coefficient of organic matter content converted to organic carbon content, SOCDi (t C ha −1 ) is the organic carbon density of the ith soil polygons, BDi (g/cm 3 ) is soil bulk density of the ith soil polygons, and H (cm) is soil depth. The units of the area of the polygons are hectare. We focused on the top 20 cm of the dry farmland soil. Calculation of SOCS (t = 10 6 g) is where Si is the area of the ith soil polygons. We then calculated the change in soil organic carbon density (ΔSOCD) (t C ha −1 )

SOCD Calculation
SOCD was selected as representative of the spatial and temporal distribution characteristics of SOC [41,42]. The SOCD is the SOCS per unit area where SOC i (g/kg) is the organic carbon content of the ith soil polygons, SOM i (g/kg) is the organic matter content of the ith soil polygons, 0.58 is the coefficient of organic matter content converted to organic carbon content, SOCD i (t C ha −1 ) is the organic carbon density of the ith soil polygons, BD i (g/cm 3 ) is soil bulk density of the ith soil polygons, and H (cm) is soil depth. The units of the area of the polygons are hectare. We focused on the top 20 cm of the dry farmland soil. Calculation of SOCS (t = 10 6 g) is where S i is the area of the ith soil polygons. We then calculated the change in soil organic carbon density (∆SOCD) (t C ha −1 ) where SOCD 1980 is the soil organic carbon density in Shandong Province in 1980, and SOCD 2016 is the soil organic carbon density in Shandong Province in 2016.

Model Validation
The SOC was simulated at the site with the data of the long-term monitoring point on dry farmland, and the simulated SOC value was obtained from 1986 to 2007. Using the dry farmland simulated SOC value and measured SOC value in the corresponding year to calculated the simulation efficiency (ME) and determination coefficient (R 2 ), to determine the simulation results of the long-term monitoring points in the dry farmland. According to the correlation between the simulated SOC value of the point and the measured SOC value of the point, we determined whether the model was suitable for the study area.
The validation point for dry farmland in Shandong Province was located in Jinan, and Table 2 contains validation point data for Jinan. The fitting between the SOC measured and simulated values uses the parameters R 2 and ME as parameters for evaluating simulation accuracy.
For the ME of the validation point measured value, the parameter is as follows For the R 2 of the validation point measured value, the parameter is as follows where P i is the SOC simulation value of the ith dry farmland validation point, O i is the measured SOC value in the ith dry farmland validation point, and O is the average value of the measured SOC value of the dry farmland validation point from the period 1986 to 2007, P is the average value of the simulated SOC value of the dry farmland validation point from 1986-2007. ME varies between 0 and 1, and the closer it is to 1, the better the simulation efficiency [43]. The closer the R 2 is to 1, the better the fitting effect, and the more realistic the fitting function.
To further verify the applicability of the CENTURY model to SOC simulation in Shandong Province, the study area multi-point measured data in 2000 was used to fit with the model simulation data in 2000, to determine the applicability of the CENTURY model in the region. In the regional multi-point verification, 39 sample data collected in Lingxian County of Shandong Province in 2000 were used to fit the model simulation data. Table 3 includes the basic point information.

Result of Point Validation
Taking the long-term dry farmland monitoring point in Jinan as a simulation point, relevant input parameters to simulate SOC dynamic changes at that point from 1986 to 2007 were fitted, and the simulated and measured values compared. The correlation analysis between simulated and measured SOCD is depicted in Figure 4.

Result of Point Validation
Taking the long-term dry farmland monitoring point in Jinan as a simulation point, relevant input parameters to simulate SOC dynamic changes at that point from 1986 to 2007 were fitted, and the simulated and measured values compared. The correlation analysis between simulated and measured SOCD is depicted in Figure 4. In Jinan, there is a strong correlation between measured and simulated SOCD (p < 0.01, R 2 = 0.722, ME = 0.37, Figure 4), indicating that the CENTURY model simulation value and the verification point data fit well, and the CENRUTY model can be used to simulate SOC evolution in Shandong Province.

Result of Multi-Point Validation
Comparing the measured SOCD data in Ling County in 2000 with the simulated values at corresponding points in 2000, we obtained Figure 5, which shows SOCD simulated and measured values at multiple points.
It was verified that the simulation efficiency of the measured data and the simulated data was 0.688, and the numerical value was between 0 and 1, indicating that the simulation results are reliable. The R 2 of multi-point linear fitting is 0.765 in the region, and it shows that the simulated SOC value and the measured SOC value fit well.  In Jinan, there is a strong correlation between measured and simulated SOCD (p < 0.01, R 2 = 0.722, ME = 0.37, Figure 4), indicating that the CENTURY model simulation value and the verification point data fit well, and the CENRUTY model can be used to simulate SOC evolution in Shandong Province.

Result of Multi-Point Validation
Comparing the measured SOCD data in Ling County in 2000 with the simulated values at corresponding points in 2000, we obtained Figure 5, which shows SOCD simulated and measured values at multiple points.

SOC Temporal Changes
From 1980 to 2016 on dry farmland in Shandong Province, SOCD varied greatly among the 2392 soil polygons. Table 4 summarizes the SOCD and related statistics of dry farmland soil polygons in 1980,1985,1990,1995,2000,2005, 2010 and 2016. From 1980 to 2016, the SOCD minimum value initially increased and then declined slightly. Maximum SOCD showed a decreasing trend, from 150.06 t C ha −1 in 1980, to 98.12 t C ha −1 in 2016. Beginning in 1980, average SOCD initially decreased and then increased, before slightly dropping back again in 2016. The differences between maximum and minimum SOCD gradually decline from 145.1 to 92.03 t C ha −1 . The annual standard deviation values of SOCD also gradually decreased from 15.11 to 11.89 t C ha −1 , indicating that the dispersion of SOCD data in Shandong Province also decreased. It was verified that the simulation efficiency of the measured data and the simulated data was 0.688, and the numerical value was between 0 and 1, indicating that the simulation results are reliable. The R 2 of multi-point linear fitting is 0.765 in the region, and it shows that the simulated SOC value and the measured SOC value fit well.

SOC Temporal Changes
From 1980 to 2016 on dry farmland in Shandong Province, SOCD varied greatly among the 2392 soil polygons. Table 4 summarizes the SOCD and related statistics of dry farmland soil polygons in 1980,1985,1990,1995,2000,2005, 2010 and 2016. Table 4. Dry farmland soil organic carbon density (SOCD) and related statistics in Shandong Province in 1980,1985,1990,1995,2000,2005,2010  From 1980 to 2016, the SOCD minimum value initially increased and then declined slightly. Maximum SOCD showed a decreasing trend, from 150.06 t C ha −1 in 1980, to 98.12 t C ha −1 in 2016. Beginning in 1980, average SOCD initially decreased and then increased, before slightly dropping back again in 2016. The differences between maximum and minimum SOCD gradually decline from 145.1 to 92.03 t C ha −1 . The annual standard deviation values of SOCD also gradually decreased from 15.11 to 11.89 t C ha −1 , indicating that the dispersion of SOCD data in Shandong Province also decreased.

Temporal Evolution of SOC in Different Soil Orders
SOCD in dry farmland soils across the province was higher in 2016 compared to 1980, where minimum SOCD increased by 0.55 t C ha −1 . Maximum SOCD decreased by 52.52 t C ha −1 , and variation was relatively small, as indicated by lower standard deviation and variance of SOCD in 2016 than in 1980. Figure 7 shows the ΔSOCD values of four main dry farmland soil orders in Shandong Province in <−10, −10-−5, −5-0, 0-5, 5-10, and >10 t C ha −1 classes.

Temporal Evolution of SOC in Different Soil Orders
SOCD in dry farmland soils across the province was higher in 2016 compared to 1980, where minimum SOCD increased by 0.55 t C ha −1 . Maximum SOCD decreased by 52.52 t C ha −1 , and variation was relatively small, as indicated by lower standard deviation and variance of SOCD in 2016 than in 1980. Figure 7 shows the ∆SOCD values of four main dry farmland soil orders in Shandong Province in <−10, −10-−5, −5-0, 0-5, 5-10, and >10 t C ha −1 classes. From Figure 7: In the ∆SOCD range < −10 t C ha −1 , inceptisols accounted for the largest number of dry farmland soil polygons (106). For ∆SOCD values from −10 to −5, there are 134 polygons of alfisols and 20 polygons of entisols. In the ∆SOCD range −5-0 t C ha −1 , the alfisols accounted for the largest number of dry farmland soil polygons (297), while the inceptisols and entisols followed with 99 and 29, respectively. In the ∆SOCD range from 0 to 5 t C ha −1 , the number of alfisols polygons was the highest, with 416. For ∆SOCD of 5-10 t C ha −1 , the alfisols and inceptisols accounted for the most, with 180 and 115, respectively. In the ∆SOCD range of >10 t C ha −1 , entisols polygons dominated, with 195. Soil polygons of the four main soil orders are mainly distributed in the −5-0, 0-5, 5-10, and >10 t C ha −1 ranges, accounting for 18.55, 30.27, 17.99, and 18.94%, of the total number of dry farmland soil polygons, respectively.

Temporal Evolution of SOC in Different Soil Orders
SOCD in dry farmland soils across the province was higher in 2016 compared to 1980, where minimum SOCD increased by 0.55 t C ha −1 . Maximum SOCD decreased by 52.52 t C ha −1 , and variation was relatively small, as indicated by lower standard deviation and variance of SOCD in 2016 than in 1980. Figure 7 shows the ΔSOCD values of four main dry farmland soil orders in Shandong Province in <−10, −10-−5, −5-0, 0-5, 5-10, and >10 t C ha −1 classes. From Figure 8: of the four studied soil orders, inceptisols SOCD was significantly higher than in the other three soil orders. From 1980 to 2016, inceptisols SOCD generally increased, by about 3.10 t C ha −1 . From 1980 to 2016, entisols SOCD was the second highest, and displayed an increasing trend, rising by 8.15 t C ha −1 , one of the largest increases in SOCD of the four main dry farmland soil orders. Changes in vertisols SOCD were similar to that of entisols, increasing by 6.84 t C ha −1 . From 1980 to 2016, alfisols was the only soil order where SOCD decreased, from 23.99 in 1980 to 23.67 t C ha −1 in 2016. It can be known that the SOCD of alfisols decreases at a rate of 0.01 t C ha −1 year −1 , which is an important factor for the C leakage of dry farmland in Shandong Province. It accounts for 35.02% of the dry farmland area. In the future, dry farmland management in Shandong Province should strengthen the rational use and optimized management of alfisols to promote the transformation of alfisols into C sinks. Using SPSS 22, the SOCD of the four major soil orders was tested by K independent samples non-parametric empirical method, and the p value of the Kruskal-Wallis Test was less than 0.05, so the SOCD of each soil orders exists significant differences. In 2016, the SOCD of the four soil orders in 2016 can be ranked as follows: Inceptisols > entisols > vertisols > alfisols. From Figure 8: of the four studied soil orders, inceptisols SOCD was significantly higher than in the other three soil orders. From 1980 to 2016, inceptisols SOCD generally increased, by about 3.10 t C ha −1 . From 1980 to 2016, entisols SOCD was the second highest, and displayed an increasing trend, rising by 8.15 t C ha −1 , one of the largest increases in SOCD of the four main dry farmland soil orders. Changes in vertisols SOCD were similar to that of entisols, increasing by 6.84 t C ha −1 . From 1980 to 2016, alfisols was the only soil order where SOCD decreased, from 23.99 in 1980 to 23.67 t C ha −1 in 2016. It can be known that the SOCD of alfisols decreases at a rate of 0.01 t C ha −1 year −1 , which is an important factor for the C leakage of dry farmland in Shandong Province. It accounts for 35.02% of the dry farmland area. In the future, dry farmland management in Shandong Province should strengthen the rational use and optimized management of alfisols to promote the transformation of alfisols into C sinks. Using SPSS 22, the SOCD of the four major soil orders was tested by K independent samples non-parametric empirical method, and the p value of the Kruskal-Wallis Test was less than 0.05, so the SOCD of each soil orders exists significant differences. In 2016, the SOCD of the four soil orders in 2016 can be ranked as follows: Inceptisols > entisols > vertisols > alfisols.  Figure 9a shows the spatial distribution of SOCD in dry farmland in Shandong Province in 1980. In 1980, high SOCD occurred in the southwest and a few areas in the northwest, and is due to the concentrated distribution of inceptisols with the highest SOCD. SOCD from 20 to 30 t C ha −1 accounted for 20.06% of the total dry farmland area, and occurred mainly in the northwest and parts of the  Figure 9a shows the spatial distribution of SOCD in dry farmland in Shandong Province in 1980. In 1980, high SOCD occurred in the southwest and a few areas in the northwest, and is due to the concentrated distribution of inceptisols with the highest SOCD. SOCD from 20 to 30 t C ha −1 accounted for 20.06% of the total dry farmland area, and occurred mainly in the northwest and parts of the northeastern coastal areas. SOCD from 10 to 20 t C ha −1 accounted for about half of the dry farmland area, and was mostly concentrated in most parts of the southeast and parts of the northeast. In 1980, there were relatively few dry farmland soils with SOCD < 10 t C ha −1 , accounting for only 5.25% of the total dry farmland soil. They were sporadically distributed in the central and southwest regions.

SOC Spatial Variation
to 30 t C ha −1 , were scattered mainly along the boundary of the province, and their area had increased significantly compared to 1980 (from 20.06% in 1980 to 36.25% in 2016, a relative increase of 16.19%). SOCDs from 10 to 20 t C ha −1 decreased from 53.25% in 1980 to 35.23% in 2016, and were mainly in most of the northwest, in the northeast and southeast, and also in the southwest. Dry farmland SOCDs < 10 t C ha −1 were scattered in a small part of the central and southwest regions, and their area declined compared with 1980. Figure 9c shows the spatial distribution of the changes in soil organic carbon (ΔSOCD) in dry farmland of Shandong Province from 1980 to 2016. ΔSOCD was divided into <−5, −5-0, 0-5, or >5 t C ha −1 . Firstly, the dry soils with ΔSOCD < −5 t C ha −1 were mainly distributed in the southwest and northwest regions, accounting for 8.04% of the total farmland soil. Those from −5-0 t C ha −1 were mainly distributed in the southwest, and coastal and other parts of the northeast, and accounted for 26.44% of the total. A large proportion of dry farmland was covered by soils from 0 to 5 t C ha −1 (31.95%), primarily distributed in the northeast, in parts of the northwest and southeast, and also a small distribution in the southwest of Shandong Province. The area of ΔSOCD with >5 t C ha −1 was 33.57%, mainly distributed in the northwest and southeast regions, with small distributions in the southwest and northeast.   In 2016, SOCD > 30 t C ha −1 were distributed mainly in the southwest, with significant increases in the northwest and a small part of the northeast (4.12% increase compared to 1980). The SOCDs from 20 to 30 t C ha −1 , were scattered mainly along the boundary of the province, and their area had increased significantly compared to 1980 (from 20.06% in 1980 to 36.25% in 2016, a relative increase of 16.19%). SOCDs from 10 to 20 t C ha −1 decreased from 53.25% in 1980 to 35.23% in 2016, and were mainly in most of the northwest, in the northeast and southeast, and also in the southwest. Dry farmland SOCDs < 10 t C ha −1 were scattered in a small part of the central and southwest regions, and their area declined compared with 1980. Figure 9c shows the spatial distribution of the changes in soil organic carbon (∆SOCD) in dry farmland of Shandong Province from 1980 to 2016. ∆SOCD was divided into <−5, −5-0, 0-5, or >5 t C ha −1 . Firstly, the dry soils with ∆SOCD < −5 t C ha −1 were mainly distributed in the southwest and northwest regions, accounting for 8.04% of the total farmland soil. Those from −5-0 t C ha −1 were mainly distributed in the southwest, and coastal and other parts of the northeast, and accounted for 26.44% of the total. A large proportion of dry farmland was covered by soils from 0 to 5 t C ha −1 (31.95%), primarily distributed in the northeast, in parts of the northwest and southeast, and also a small distribution in the southwest of Shandong Province. The area of ∆SOCD with >5 t C ha −1 was 33.57%, mainly distributed in the northwest and southeast regions, with small distributions in the southwest and northeast.

Discussion
Applying the CENTURY model, we simulated SOCD and SOCS evolution in dry farmland from 1980 to 2016 in Shandong Province. Research initiatives, however, often have certain limitations and uncertainties, and our study is no different.
Firstly, the study is limited to dry farmlands and is based on soil polygons superimposed onto the land use and soil layers. It did not fully consider land use change caused by other factors, such as the acceleration of urbanization in recent years. This is mainly because the article calculates the SOCD and SOCS of dry farmland, using dry farmland soil polygons as the basic simulation unit. If the land use and soil layers in different years are superimposed to obtain the corresponding simulated basic unit, multiple simulations are required. This will bring a lot of work, and the corresponding verification data are also lacking. It is difficult to estimate the spatial changes in SOCD and SOCS based on years of land use.
Secondly, the CENTURY model itself has certain deficiencies. For example, it uses monthly and annual steps, so simulated crop growth is relatively long, whereas actual growth and changes in biomass occur at much finer time scales [40]. To overcome deficiencies in the model, it could be improved accordingly, for example, by reducing the simulation step size, increasing the amount of simulation calculation, and obtaining more accurate simulation results. Bista et al. [44] used the DayCENT model to simulate the changes in soil organic carbon under various management conditions, and showed that between 1931 and 2080, the conventional dry wheat winter-summer fallow system lost 8.66 × 10 4 -21.92 × 10 4 t C ha −1 of SOCD. Furthermore, the model also has certain operational problems. For example, in 2011, Maren Oelbermann et al. [45] evaluated the SOCS of 19-year-old tropical and 13-year-old temperate alley and single crops using the CENTURY model, but SOC simulated values were lower than measured values because the model failed to account for change in soil BD. Silver et al. [46] used field data and the CENTURY model to estimate reserves of underground C in highly weathered soils in Amazonian forests. They found that simulated SOCD of sandstone was half that of clay, but the measured values indicated that it was 89% of the clay SOCD. Therefore, incorporating initial SOCD and soil clay content into the model input parameters is important for reducing simulation uncertainty [47].
Next, under the current dry farmland management measures in Shandong Province, the soil appears as a C sink. However, the study did not fully consider the greenhouse gas emissions caused by management measures. Studies show that greenhouse gas leakages can offset −241-660% of the benefits of soil C sequestration [48]. For example, returning straw to the field has been vigorously promoted in Shandong Province. In this study, returning straw to the field in Shandong Province increased SOCS, but some studies have suggested that returning straw to the field may not necessarily increase SOC. Even if C content increased, methane produced by returning straw to the field might be discharged into the atmosphere, greatly offsetting the benefits of soil C sequestration [49][50][51]. Therefore, when considering the impact of management measures on soil C sequestration, the impact of greenhouse gas emissions from management measures should be fully considered. Then, adequate model simulation requires the input of appropriate and accurate data for parameters such as climate, soil texture, soil management measures, and so on. These data are often difficult to obtain, and may not even exist, especially for long time periods. Furthermore, the data itself may have inaccuracies in the measurement method, which also increases the uncertainty of the model simulation.
Finally, the calculation of SOC will also bring uncertainty to the research. Liang et al. [52] show that the calculation of China's terrestrial SOCS is also affected by different estimation methods. The difference stems from different SOCD values, different soil depths, calculation formulas, basic data, and different mapping scales. The resulting estimation uncertainty is approximately 25.5-30.4%.
Based on the above research uncertainties, the following prospects can be made for future research. Using the Daycent model, the regional SOC was simulated, using GIS digital soil maps, a 1 m soil depth, and no gravel larger than 2 mm. For the SOCD content calculation formula, you can try to measure SOC density by remote sensing spectrum for model verification [53][54][55][56].
It is instructive to compare our results with other studies, to highlight reasons for the differences between the studies. We showed that SOCS in dry farmland (0-20 cm) in Shandong Province is about 290.58 × 10 6 t, or about 0.29% of the national SOCS [23], and 2.88% of the SOCS of dry farmland in China [57], which is less than the estimated surface SOCS in Shandong Province (350.65 × 10 6 t) by Dai [58]. These differences maybe, in part, due to differences in the scope of the two studies. The SOCD of dry farmland (0-20 cm) in Shandong Province is 25.94 t C ha −1 , which is less than the average value (29.00 t C ha −1 ) estimated by Liu et al. [47]. The two study research methods were the same, but the research objectives differed, focused on dry farmland soil in our study, and all soils in Shandong Province in Liu's research. Other differences may be caused by other soils or different research methods. For example, our average SOCD was greater than the 22.20 t C ha −1 estimated by Dai et al. [59] for surface SOCD in Shandong Province.
Compared with other regional studies, SOCD in the surface soil of dry farmland in Shandong Province is lower than the national average value of 105.3 t C ha −1 [60], 48.6 t C ha −1 in the arid region in western China [28], 33.3 t C ha −1 in Guangxi Province [61] and 74.6 t C ha −1 in Henan Province [62], but similar to the average of 29.29 t C ha −1 in western Chongqing [63]. Therefore, the relatively low SOCD in dry farmland soils in Shandong Province also means that the study area has huge soil C sequestration potential.
SOCD comparisons of dry farmland with the SOCD of other land uses. The main land use type in Shandong Province is dry farmland, and there are relatively few studies on SOC for other land uses. Therefore, we compared the dry farmland SOC in Shandong Province with the land use types in other regions.
Xie et al. [57] found that China's paddy field SOCD was 97.6 t C ha −1 in 1980, which is much higher than that of dry farmland soil in Shandong Province. Wen et al. [64] found paddy field SOCD of 46.21 t C ha −1 in northern Guangdong. Ma et al. [65] estimated that the paddy field SOCD in Taoyuan County is 43.5 t C ha −1 . However, Xu et al. [66] found that the SOCD of paddy soil in Anhui Province was (19 ± 4.91) t C ha −1 , which is lower than that of dry farmland soil in this study. The SOCD of woodlands and forests is usually higher than that of dry farmland soil. The SOCD of eucalyptus, Taxodiaceae and Masson pine plantations in South China was 32.9, 36.0, and 38.6 t C ha −1 , respectively [67]. The average SOCD of the forests in the northeast region is 217.8 t C ha −1 [68]. Yang et al. [69] estimated that the SOCD of the soil in Ziwuling forest area of Loess Plateau is 105.2 t C ha −1 . Grassland SOCD in different areas varies greatly, and the contemporary grassland SOCD in Inner Mongolia is 19.9 t C ha −1 [70]. There are also 299.7 t C ha −1 of SOCD in the source region of the Yellow River [71].
We found that SOCD and SOCS in Shandong Province showed an initial increasing trend, and then slowly decreased from 1980 to 2016. The SOCD decrease from 26.47 t C ha −1 in 2009 to 25.94 t C ha −1 in 2016 is noteworthy and may be explained in the following way: firstly, the input of external SOM had continued to increase in the previous period, but the external SOM input after 2009 did not increase, resulting in greater SOM decomposition than the accumulation of SOC; and secondly, with time, mechanized farming in Shandong Province continues to rise. The main dry farming crops have a mechanical harvesting yield of more than 90% [72], but the quality of mechanized farming is poor, leading to greater soil destruction and accelerated SOC decomposition.

Conclusions
Although we show that, from 1980 to 2016, the dry farmland in Shandong Province was generally a C sink (26.99 × 10 6 t C), the changing SOCD and SOCS trends suggest that, under existing management measures, climatic conditions, etc., remain unchanged, it will be converted into a C source. Therefore, corresponding changes should be made to increase the input of external organic matter, including the amount of organic fertilizer that was applied, increasing the amount of straw returned to the field, and adopting conservation tillage.
The four main soil orders in Shandong Province, ranked by the size of the area covered, can be ranked as: Inceptisols > alfisols > entisols > vertisols. However, based on SOCD they are ranked as: inceptisols > entisols > vertisols > alfisols. Inceptisols contribute the most to the SOC in dry farmland of Shandong Province, and there are significant differences in SOCD among various soil orders.
SOCD gradually increases from the northeastern coast to the inland southwest. In 1980, dry farmland SOCD of Shandong Province was >30 t C ha −1 , and mainly concentrated in the southwest and small parts of the northwest, except for some small coastal areas. Compared with 1980, the SOCD in Shandong Province increased significantly in 2016, especially in the west and in coastal areas. Dry farmland SOCD > 30 t C ha −1 were mainly concentrated in the southwest and northwest, while those in the range of 20-30 t C ha −1 were mainly distributed in the northwest. Declines in SOCD were concentrated mainly in the southwest, although here SOCD was initially higher. The SOCD area increased significantly.
The area with reduced SOCD accounted for 34.48% of the dry farmland area in Shandong Province, and the range of reduced SOCD was mainly concentrated from −5 to 0 t C ha −1 (26.44%). The soils with ∆SOCD in the range of 0-5 and >5 t C ha −1 accounted for 31.95% and 33.57% of the dry farmland area, respectively.