According to the bulletin of the third national land survey, in 2020, Chinese grasslands covered an area of 2.645 × 10⁶ km², of which 80.59% were natural pastures and 0.22% were cultivated grasslands. Grasslands in China are mainly distributed in six provinces (Tibetan Autonomous Region, Inner Mongolia Autonomous Region, Xinjiang Uygur Autonomous Region, Qinghai Province, Gansu Province, and Sichuan Province), accounting for 94% of the total grassland area. China is a vast country with complex and diverse climate types, and SOCD is unevenly distributed (Zhang Y et al., 2022). To better investigate the distribution and variation of SOCD, based on climate conditions, Chinese grasslands were divided into four regions in this study: northwestern, northern, southern, and Qinghai–Tibetan regions ( Figure 1 ). The northern and southern regions are in the north and south of the monsoon climate zone, respectively. The northwestern region is in the non-monsoon climate zone, and the Qinghai–Tibetan region is a unique geographical area with a distinctive climate type (Zhang Y et al., 2022).

The basic flow of SOCD simulation modeling and subsequent analysis are presented in Figure 2 . Based on 552 observed data and the modeling factors, we built a grassland SOCD prediction model using an RF model. During the training run, the model was adjusted for 200 iterations, and we used 10-fold cross-validation to evaluate the accuracy of the model training results. After the model training, we input 15 modeling factor layers and simulated SOCD on each raster to obtain the national SOCD layers of 1982, 1990, 2000, 2010, and 2020. At the same time, we calculated the interquartile (the difference between the 75th and 25th percentiles) of 200 iterations of simulated SOCD for each raster, which was used to represent the uncertainty of the simulated SOCD (Ding et al., 2016) ( Supplementary Figure 1 ). The Chinese grassland SOCD layers were extracted by mask using the Chinese grassland vector layers. Following Liang et al. (2019), we examined the spatial simulation performance and validated the spatial accuracy of the grassland SOCD generated from the RF model using 100 measured samples from the original data and their corresponding values from the 2010 grassland SOCD layer. As there were few measured samples for a single year to be obtained, the temporal accuracy of the grassland SOCD generated from the RF model was validated as follows: we selected the measured samples between 2000 and 2005 (n = 100) and treated these soil samples as the measured samples in 2000. Then, we extracted the values of grassland SOCD layer in 2000 and their corresponding measured samples to test the performance of the model following Li H et al. (2022).

The observed data were used as the initial data to determine the contribution of each factor to grassland SOCD using the RF model. The factors with cumulative contributions larger than 70% were the dominant factors driving grassland SOCD change. For the dominant dynamic factors, using the grassland vector layers, we mask-extracted these factor layers to obtain the factor layers for the grassland section and used the zonal statistics (mean value calculation) to calculate the mean values of factors for the four regions (northwestern, northern, southern, and Qinghai–Tibetan regions). Then, we examined the trends of these factors over the past 39 years in four regions using the MK test. The mask extraction and zonal statistics were carried out on ArcGIS 10.2. The RF model and MK trend test were run on the R platform (R version 4.1.2, R Core Team, 2022), using the “randomForest” (Liaw and Wiener, 2002) and the “trend” (Pohlert, 2020) package, respectively.

The accuracy of the RF model training dataset and the results of 10-fold cross-validation are shown in Figure 3 . The R ² of the training dataset was 0.923 and the RMSE was 3.384 kg C m⁻² ( Figure 3A ). The results of the 10-fold cross-validation showed that the R ² of the model was 0.674 and the RMSE was 7.096 kg C m⁻² ( Figure 3B ). These results indicated a good correlation between observed SOCD and predicted SOCD. The R ² and RMSE of the spatial accuracy of the grassland SOCD prediction are shown in Figure 4A , where the R ² and RMSE were 0.239 and 13.302 kg C m⁻², respectively. The temporal accuracy of the grassland SOCD prediction is presented in Figure 4B , with an R ² of 0.256 and an RMSE of 6.063 kg C m⁻². According to the contribution ranking of the RF model, we considered the top 4 factors (cumulative contribution > 70%) as the dominant factors driving changes in grassland SOCD, and they were MAT, NDVI, elevation, and MAWS ( Table 2 ).

At the national level, grassland SOCD simulated by the RF model was low in the west and north and high in the east and south ( Figure 5 ). The high SOCD values were distributed in the eastern Qinghai–Tibetan region and near the Tianshan Mountains of the Xinjiang Uygur Autonomous Region, while the low SOCD values were mainly located in the western Qinghai–Tibetan region and central northwestern region. In the southern and northern regions, the values of grassland SOCD generally ranged from 4 to 12 kg C m⁻². In the northwestern region, the SOCD values were relatively high in the eastern Inner Mongolia and northern Xinjiang (greater than 4 kg C m⁻²), and relatively low in the rest of this region (less than 4 kg C m⁻²). The spatial heterogeneity of grassland SOCD in the Qinghai–Tibetan region was large, showing a gradual increase from northwest to southeast. In this region, SOCD values were higher than 12 kg C m⁻² in southern Tibet, southern Qinghai, and western Sichuan, and the values were lower than 12 kg C m⁻² in the rest of the region. During the period 1982–2020, the largest net increase of SOCD occurred in eastern Inner Mongolia, with an increase over 4 kg C m⁻². Conversely, the SOCD decreased the most in southern Tibet where the reduction was greater than 4 kg C m⁻².

During the study period, the mean grassland SOCD increased in each period with respect to 1982 across China, with a net increase of 0.734 kg C m⁻² from 1982 to 2020 ( Table 3 ). However, variation in the mean grassland SOCD among regions was different over time. In the northern region, the mean grassland SOCD first decreased and then increased, and it was 9.471 kg C m⁻² in 1982 and 9.299 kg C m⁻² in 2020, decreasing by 0.172 kg C m⁻². The mean grassland SOCD in the northwestern, Qinghai–Tibetan, and southern regions showed a net increase of 1.440 kg C m⁻², 0.915 kg C m⁻², and 0.411 kg C m⁻², respectively. In the three regions, the grassland SOCD of each period showed varying degrees of increase compared to 1982.

Chinese grassland area showed a decreasing trend between 1982 and 2020, with a decrease of 3.985 × 10⁵ km² ( Table 4 ). The most evident decrease in grassland area was in the Qinghai–Tibetan region, with a decrease of 2.494 × 10⁵ km². The least decrease in the area was observed in the southern region, with a reduction of 1.318 × 10⁴ km². Over the whole study period, an area of 8.549 × 10⁵ km² of grassland was mainly converted to unused land (58.14%), woodland (20.22%), and cropland (14.82%) ( Table 5 ). In addition, there was also a small proportion of grassland converted to water body and built-up land.

The total SOCS of Chinese grasslands showed a change featuring first an increase and then a decrease ( Table 6 ). The increase of grassland SOCS was obvious during 1982–1990 with an increase of 1.352 Pg, then gradually decreased by 2.509 Pg during 1990–2020. During the study period, the net decrease of 1.158 Pg in grassland SOCS was mainly caused by the reduction from 1990 to 2020. In four regions, from 1982 to 2020, grassland SOCS of the northwestern region showed an overall upward trend, with a net increase (0.950 Pg), and it showed a net decrease in the southern (0.005 Pg), northern (0.618 Pg), and Qinghai–Tibetan regions (0.923 Pg), respectively. Among these three regions, grassland SOCS in the southern and Qinghai–Tibetan regions increased slightly during 1982–2000 and then decreased until 2020, while the northern region experienced a gradual decrease from 1982 to 2020.

Chinese grassland SOCD was high in the east and south, and low in the west and north ( Figure 5 ), a similar observation reported by Liu et al. (2022). Generally, SOCD is influenced by biological productivity and organic matter mineralization, as well as controlled by hydrothermal conditions (Duchaufour, 1983; NSSO, 1998). In the southern region of the study area, the climate is more humid due to monsoonal circulation (Zhang, 1991; Wu and Peng, 2003). On the contrary, in the northwestern region, the climate condition is arid due to the far distance from the ocean and the blockage of the Qinghai–Tibetan plateau (Wu and Peng, 2003). In the southern and northern regions with a warm climate and abundant precipitation, improved vegetation growth and high primary productivity promote the accumulation of SOC and contribute to higher grassland SOCD (Ontl and Schulte, 2012). In northwestern China, the western part is mainly covered by desert, while the eastern part is mostly covered by grasslands. The dry climatic conditions and scarce precipitation are not conducive to vegetation growth and high NPP, which might contribute to relatively low SOCD in the western part (Li et al., 2004; Ontl and Schulte, 2012; Zhang et al., 2016). The Qinghai–Tibetan region has unique climatic conditions, and the elevation of this region is high in the west and low in the east (Zhang Y et al., 2022; Liu et al., 2019). The cold and relatively humid conditions promoted the accumulation of soil organic matter in the southeastern of the Qinghai–Tibetan region, resulting in higher SOCD (Piao et al., 2011; Zhong et al., 2012).

The mean SOCD of Chinese grasslands was 7.791 kg C m⁻² in 1982 and 8.525 kg C m⁻² in 2020, showing a net increase of 0.734 kg C m⁻² ( Table 3 ). This was consistent with our hypothesis that grassland SOCD showed a net increasing trend. The results from the RF model suggested that MAT, NDVI, elevation, and MAWS were the main factors driving the variation in grassland SOCD across China, with MAT being the most important factor, explaining 44.00% of the total variation, followed by NDVI, elevation, and MAWS, explaining 13.45%, 10.62%, and 5.16%, respectively ( Table 2 ).

Grassland SOCD and grassland area are important factors influencing changes in grassland SOCS. A reduction in grassland area would lead to the loss of grassland SOCS and make grasslands a source of greenhouse gas emissions. Between 1982 and 2020, the Chinese grassland area declined severely, with a total reduction of 3.985 × 10⁵ km², resulting in a 1.158 Pg reduction in grassland SOCS ( Tables 4 and 6 ). This was contrary to our hypothesis that SOCS of Chinese grasslands increased between 1982 and 2020. The transfer of grasslands to unused lands was the main reason for the reduction in grassland area, accounting for 58.14% of the total changed area ( Table 5 ; Figure 9 ). Various factors, such as population growth, overgrazing, emphasis on livestock over grass, frequent droughts, insects, and other natural disasters, accelerate the process of grassland degradation and desertification (Bardgett et al., 2021; Li et al., 2021). These led to the conversion of some grasslands to unused lands (gobi, sandy land, swampland, gravel and rock, bare land, and other unused lands) and resulted in the loss of SOCS. Furthermore, the cultivation of grassland, the construction of reservoirs, and the development of urban areas led to the conversion of some grasslands to croplands, water bodies, and built-up lands (Bardgett et al., 2021; Chang et al., 2022; Li et al., 2016a). Different land types have different SOCDs, and it is generally accepted that woodlands have the highest SOCD, followed by grasslands, croplands, and built-up lands (Li Y et al., 2022; Lin et al., 2021; Zhu et al., 2021). Therefore, the conversion of grasslands to woodlands with high SOCD usually promoted the accumulation of SOCS. In contrast, the conversion of grasslands to land types with low SOCD (e.g., croplands, water bodies, and built-up lands), resulted in a loss of SOCS, leading to significant CO₂ emission to the atmosphere and accelerating global warming.

As mentioned above, the SOCD of grasslands in China gradually increased during the study period. This implies that more carbon can be sequestrated into the soil of grasslands. However, due to the reduction in grassland area caused by degradation and conversion, grassland SOCS decreased across China in general. Therefore, taking certain measures to mitigate the decline in grassland areas appears to be crucial in increasing grassland SOCS.

To do this, China has implemented several ecological restoration projects since the 1980s, including the Grain to Green Program introduced in 1999 and the Returning Grazing Land to Grassland Program introduced in 2003 (Li et al., 2016a; Zhou W et al., 2017). These programs were aimed to convert agricultural lands into forests and grasslands (Ferraro and Kiss, 2002) and to relieve grazing pressure on degraded grasslands through grazing bans, fallow grazing for grass, fencing, and rotational grazing (Tong et al., 2004). It was shown that the carbon sequestration amount due to these programs’ implementation since 2001 and 2010 was 0.200 Pg and 0.118 Pg, respectively (Lu et al., 2018). Various grassland protection programs, such as banning grassland reclamation and private conversion and destruction of grasslands, were also enforced (Xing, 2016). These measures not only were effective in alleviating the degradation of grasslands and the reduction of SOCS due to grassland area decrease (Li et al., 2016a; Zhou W et al., 2017), but also enhanced the grassland’s ability against wind erosion through increasing vegetation cover (Shao et al., 2016).

Besides this, as temperature and altitude influenced grassland SOCD change, additional work should be done in the future. For example, we should pay more attention to temperature or altitudes that caused a decrease in SOCD. We should try to restrict change from carbon sources to pools in such areas by implementing proper ecological engineering or adopting measures such as water and soil conservation and regulating grassland use practices. However, rarely has any study looked at the improvement of grassland SOCD in specific temperatures and altitudes.

The RF model has been used extensively for modeling. For example, Li H et al. (2022) and Sreenivas et al. (2014) used it to quantify the dynamics of SOC in China and SOCD in southern India, respectively. Their findings demonstrated the feasibility of the model in predicting SOC dynamics. In this study, the results also indicated that the RF model has good spatiotemporal predictive power ( Figure 4 ). However, there are still some uncertainties in the model.

First, there may be more factors that contribute to SOCD changes but were not included in our RF model. These factors include litter quality, nitrogen deposition, and CO₂ enrichment (Zhang W et al., 2022; Sumiyoshi et al., 2017). Future studies should focus on including more factors in the model and intensive sampling in sparsely sampled areas to reduce the uncertainty of model predictions (Liang et al., 2019). Second, the data sources are not consistent. For example, the data sources of some modeling factors used to predict grassland SOCD in 2020 were not consistent with the previous data sources, which may affect the model prediction results (Zhang et al., 2017). Third, we only have the observed soil samples data in 2010 to build a model, which may cause some uncertainty in the results over a long time series of simulations. Fourth, while the datasets containing 552 samples across the study area are large enough to increase the performance of the model and extrapolability of the results, assessment of the spatial coverage of samples was not carried out and is beyond the scope of the study objective. Therefore, in future studies, assessment of spatial coverage of sample locations to capture the variability of environmental conditions within the four regions (e.g., the southern, northern, northwestern, and Qinghai–Tibetan regions) should be considered to increase the extrapolability of the results.