Soils not only provide nutrients such as nitrogen and phosphorus for various terrestrial plants but also contain a large amount of organic carbon (Zhao et al., 2022; Niu and Fu, 2024), which is an important carbon reservoir for regulating global warming. Soil organic carbon (SOC), total nitrogen (TN), and total phosphorus (TP), and their ratios (i.e., ratio of SOC to TN, C:N; ratio of SOC to TP, C:P; and ratio of TN to TP, N:P), are important soil variables (Yang et al., 2007; Zhang et al., 2007; Feng et al., 2017; Liu et al., 2019; Kou et al., 2020; Li et al., 2020; Li et al., 2021). At least one of SOC, TN, TP, C:N, C:P, and N:P is closely correlated with soil respiration (Yu et al., 2019b; Fu and Shen, 2022), N₂O flux (Zhang et al., 2020; Bahram et al., 2022), forage nutrition quality and production (Zha et al., 2022; Han et al., 2023b), soil microbial diversity, and plant diversity (Yu et al., 2019a; Wang et al., 2021b; Zhang and Fu, 2021). Therefore, modeling SOC, TN, TP, C:N, C:P, and N:P on multiple space-time scales is an important challenge for global carbon, nitrogen, and phosphorus cycling. To obtain these six soil variables, predecessors have acquired various methods, including direct observation methods (Fu et al., 2012; Han et al., 2023b), spatial interpolated methods (de Melo et al., 2016; Shit et al., 2016), process models (e.g., CENTURY model) (Mikhailova et al., 2000; Tornquist et al., 2009; Zhuang et al., 2010), and machine learning models (e.g., random forest model and eXtreme gradient boosting model) (Bangelesa et al., 2020; Marcal et al., 2021). These previous studies can provide very important guidance for our current and future research but they have the following two deficiencies. First, of the numerous existing machine learning methods (Reichstein et al., 2019; Hutson, 2022), there is still controversy over which is best at quantifying soil carbon, nitrogen, and phosphorus and their ratios (Bangelesa et al., 2020; Wang et al., 2021a). Second, there are always a lot of model input variables for many of these previous studies (Wang et al., 2020). However, the performance of the model does not always increase with the increase of input variables due to the data accuracy of input variables themselves and other potential reasons (Dai et al., 2023; Wang and Fu, 2023). The number of input variables also affects the speed of the model, increasing costs such as electricity (Tian and Fu, 2022; Wang and Fu, 2023). In addition, any input variable depends on raw observation data, which incurs labor costs and requires resources (Tian and Fu, 2022; Dai et al., 2023). By contrast, models based on a smaller number of input variables but with a higher accuracy may be the potentially optimal model, although this requires further proof. Therefore, further research about finding the best models of the six soil variables is necessary.

Alpine soils on the Xizang Plateau are not only important components of global alpine soils but also a very sensitive reservoir of soil carbon under the background of global change. Accurate quantification of soil carbon is an important basis for the accurate quantification of the carbon sink function of terrestrial ecosystems under the background of global change on the Xizang Plateau, and a highly accurate model is an important basis for the accurate quantification of the soil carbon pool. With the development of machine learning techniques, an increasing number of tools are available for us to choose for the modeling of environmental variables (Tian and Fu, 2022; Dai et al., 2023; Wang and Fu, 2023). Under such a scenario, some studies have tried to find the best models of α-diversity, forage nutrition quality and storage, soil moisture, and pH on the Xizang Plateau from a variety of machine learning methods (Han et al., 2022; Tian and Fu, 2022; Dai et al., 2023; Wang and Fu, 2023). These studies have confirmed that the random forest model (RFM) can be the best method (Han et al., 2022; Tian and Fu, 2022; Dai et al., 2023; Wang and Fu, 2023). However, whether the RFM is better than other methods at quantifying SOC, TN, TP, C:N, C:P, and N:P for grasslands on the Xizang Plateau is unclear. On the other hand, although it is well known that climate change and human activities jointly affect grassland ecological systems (e.g., SOC) (Ganjurjav et al., 2015; Wang et al., 2022), it is still controversial whether climate change or human activities play a dominant role in the change of grassland ecosystems on the Xizang Plateau. To separate the relative impacts of climate change and human activities on SOC, TN, TP, C:N, C:P, and N:P, it is necessary to construct the six variable models driven by pure climate variables and models jointly affected by climate change and human activities, respectively. Therefore, further studies are needed to quantify the six variables by choosing the best machine learning method in grassland areas on the Xizang Plateau.

In this study, we wanted to find the best models of SOC, TN, TP, C:N, C:P, and N:P at three depths (0–10 cm, 10–20 cm, and 20–30 cm) using only three or four input variables under free-grazing or non-grazing scenarios in the Xizang grassland ecological systems from nine models. The nine models included the RFM, generalized boosted regression model (GBRM), support vector machine model (SVM), multiple linear regression model (MLRM), recursive regression tree model (RRTM), artificial neural network model (ANNM), generalized linear regression model (GLRM), conditional inference tree model (CITM), and eXtreme gradient boosting model (eXGBM) (Han et al., 2022; Wang and Fu, 2023). Previous studies have confirmed that the RFM more accurately models α-diversity and forage nutrition quality, soil pH, and soil moisture in the Xizang grasslands (Han et al., 2022; Tian and Fu, 2022; Dai et al., 2023; Wang and Fu, 2023). Therefore, we hypothesized that the RFM will more accurately model SOC, TN, and TP, and their ratios, in this study.

The study area was located in the grassland areas of the Xizang (Supplementary Figure S1). From west to east, there is a successive distribution of alpine desert steppes, alpine steppes, and alpine meadows (Supplementary Figure S1). All soils were obtained using a soil auger, immediately put in the car refrigerator, and were then transferred to a lab freezer (−20°C) for storage before soil analyses. Soil samples at a depth of 0–10 cm (elevation, 4,279 m to 5,261 m; longitude, 79.46°E to 92.01°E; and latitude 29.28°N to 33.23°N) under non-grazing scenarios were collected in 2011 and 2013–2020. Soil samples at a depth of 10–20 cm (elevation, 4,279 m to 5,261 m; longitude, 79.46°E to 92.01°E; and latitude 29.28 °N to 33.23°N) under non-grazing scenarios were collected in 2011, 2013, 2015, 2017–2018, and 2020. Soil samples at a depth of 20–30 cm (elevation, 4,279 m to 5,261 m; longitude, 84.82°E to 91.07°E; and latitude 29.28 °N to 32.00°N) under non-grazing scenarios were collected in 2011 and 2020. Soil samples at a depth of 0–10 cm (elevation, 2,785 m to 5,330 m; longitude, 79.46°E to 95.68°E; and latitude 28.37°N to 33.23°N) under free-grazing scenarios were collected in 2013 and 2015–2020. Soil samples at a depth of 10–20 cm (elevation, 2,785 m to 5,330 m; longitude: 79.46°E to 95.68°E; and latitude 28.37°N to 33.23°N) under free-grazing scenarios were collected in 2013, 2015, and 2017–2020. Soil samples at a depth of 20–30 cm (elevation, 4,300 m to 5,330 m; longitude, 84.82°E to 91.07°E; and latitude 29.22°N to 32.15°N) under free-grazing scenarios were collected in 2019–2020.

We measured soil organic carbon (SOC), total nitrogen (TN), and total phosphorus (TP) based on the soil samples mentioned above, and then calculated the ratio of SOC to TN (C:N), ratio of SOC to TP (C:P), and ratio of TN to TP (N:P) (Sun et al., 2021; Zha et al., 2022). The potassium dichromate method, Kjeldahl method, and molybdenum antimony resistance colorimetry method were used to analyze SOC, TN, and TP, respectively (Sun et al., 2021). The SOC under non-grazing scenarios at 0–10, 10–20, and 20–30 cm was labeled by SOC_{p_0–10}, SOC_{p_10–20}, and SOC_{p_20–30}, but the SOC under free-grazing scenarios at 0–10, 10–20, and 20–30 cm was labeled by SOC_{a_0–10}, SOC_{a_10–20}, and SOC_{a_20–30}, respectively. Similarly, we labeled TN_{p_0–10}, TN_{p_10–20}, TN_{p_20–30}, TN_{a_0–10}, TN_{a_10–20}, TN_{a_20–30}, TP_{p_0–10}, TP_{p_10–20}, TP_{p_20–30}, TP_{a_0–10}, TP_{a_10–20}, TP_{a_20–30}, C:N_{p_0–10}, C:N_{p_10–20}, C:N_{p_20–30}, C:N_{a_0–10}, C:N_{a_10–20}, C:N_{a_20–30}, C:P_{p_0–10}, C:P_{p_10–20}, C:P_{p_20–30}, C:P_{a_0–10}, C:P_{a_10–20}, C:P_{a_20–30}, N:P_{p_0–10}, N:P_{p_10–20}, N:P_{p_20–30}, N:P_{a_0–10}, N:P_{a_10–20}, and N:P_{a_20–30}.

The observed SOC_{a_0–10}, TN_{a_0–10}, TP_{a_0–10}, C:N_{a_0–10}, C:P_{a_0–10}, and N:P_{a_0–10} was 1.27–206.28 g kg⁻¹, 0.10–10.30 g kg⁻¹, 0.11–5.44 g kg⁻¹, 4.37–21.30, 2.75–262.65, and 0.34–13.46, respectively. The observed SOC_{a_10–20}, TN_{a_10–20}, TP_{a_10–20}, C:N_{a_10–20}, C:P_{a_10–20}, and N:P_{a_10–20} was 1.41–117.09 g kg⁻¹, 0.13–8.87 g kg⁻¹, 0.14–4.27 g kg⁻¹, 0.66–22.08, 3.18–249.60, and 0.38–16.69, respectively. The observed SOC_{a_20–30}, TN_{a_20–30}, TP_{a_20–30}, C:N_{a_20–30}, C:P_{a_20–30}, and N:P_{a_20–30} was 1.51–58.47 g kg⁻¹, 0.27–4.92 g kg⁻¹, 0.18–0.96 g kg⁻¹, 3.61–16.42, 2.98–158.50, and 0.68–10.33, respectively. The observed SOC_{p_0–10}, TN_{p_0–10}, TP_{p_0–10}, C:N_{p_0–10}, C:P_{p_0–10}, and N:P_{p_0–10} was 1.89–80.23 g kg⁻¹, 0.30–5.57 g kg⁻¹, 0.18–0.75 g kg⁻¹, 3.47–21.37, 6.16–147.40, and 0.98–11.39, respectively. The observed SOC_{p_10–20}, TN_{p_10–20}, TP_{p_10–20}, C:N_{p_10–20}, C:P_{p_10–20}, and N:P_{p_10–20} was 1.76–43.59 g kg⁻¹, 0.37–4.20 g kg⁻¹, 0.15–0.85 g kg⁻¹, 2.60–14.83, 4.86–152.20, and 1.02–11.44, respectively. The observed SOC_{p_20–30}, TN_{p_20–30}, TP_{p_20–30}, C:N_{p_20–30}, C:P_{p_20–30}, and N:P_{p_20–30} was 1.16–32.74 g kg⁻¹, 0.34–3.32 g kg⁻¹, 0.11–0.85 g kg⁻¹, 2.04–13.13, 3.18–56.23, and 1.00–7.50, respectively. Referring to previous studies (Dai et al., 2023; Wang and Fu, 2023), these soil carbon, nitrogen, and phosphorus datasets were randomly divided into two parts using the sample function of the R software, one of which contains 30 observations for model accuracy testing and the remaining observations for model construction.

Previous studies found that only three climate variables (air temperature, precipitation, and radiation) can obtain high-precision random forest models of soil moisture and pH under fencing conditions, while combining the three climate variables and normalized difference vegetation index (NDVI) can obtain high-precision random forest models of soil moisture and pH under grazing conditions (Dai et al., 2023; Wang and Fu, 2023). In addition, these four input variables are relatively easy to obtain, and the model accuracy does not necessarily improve with more input variables (Dai et al., 2023; Wang and Fu, 2023). In this study, only these four independent variables were used to model soil carbon, nitrogen, and phosphorus.

Referring to previous studies (Tian and Fu, 2022; Dai et al., 2023), we modeled SOC, TN, TP, C:N, C:P, and N:P at 0–10, 10–20, and 20–30 cm under non-grazing scenarios using annual temperature (AT), annual precipitation (AP), and annual radiation (ARad) from the nine models. Although the default number of trees in the random forest model is 500, this may not be the optimal value; therefore, this default value was not adopted in this study. This may improve the accuracy of the model. The AT, AP, and ARad data were based on monthly air temperature, precipitation and radiation, respectively (Fu et al., 2022; Wang and Fu, 2023). Monthly climate data with high accuracies were spatial interpolation data from 145 weather stations, and their spatial resolution was 1 km × 1 km (Wang and Fu, 2023).

All the nine models were performed using R.4.1.2. The RFM, RRTM, MLRM, GBRM, and SVMM were based on the randomForest, rpart, stats, gbm, and e1071 packages, respectively (Han et al., 2022; Han et al., 2023a; Dai et al., 2023). The other four models were based on the rminer package (Han et al., 2022; Dai et al., 2023). For the RFM, except ntree and mtry, parameters of the randomForest function were set to default values. To obtain the optimal combination of the ntree and mtry, first we set the range of mtry to 1–3 or 1–4 under fencing or grazing conditions, and the ntree to 1,000, and ran 300 or 400 random forests to filter out the optimal ntree. Then we ran 300 or 400 random forest operations based on the mtry (1–3 or 1–4) and optimal ntree and selected the combination of ntree and mtry based on the largest R ² value and the smallest mean square errors as the final random forest model. For the GBM, first the gbm and gbm.perf functions together were used to find the optimal n.trees. The cv.folds and n.trees of the gbm function were set as 2 and 10,000, respectively, and all the other parameters were set to be default values. The method of the gbm.perf function was set as cv, and all the other parameters were set to default values. Then, we used the gbm function to obtain the final gbm model based on the cv.folds of 2 and optimal n.trees. For the MLRM, RRTM, and SVM, all the parameters of the lm function, rpart function, and svm function were set to default values, respectively. For the ANNM, GLRM, CITM, and eXGBM, the search and scale parameters of fit function were set as heuristic5 and none, and the model parameter was set as mlp, cv.glmnet, ctree, and xgboost, respectively. In addition, all the other parameters were set to default values for the ANNM, GLRM, CITM, and eXGBM.

Similarly, we modeled SOC, TN, TP, C:N, C:P, and N:P at 0–10, 10–20, and 20–30 cm under free-grazing scenarios using AT, AP, ARad, and growing-season maximum normalized difference vegetation index (NDVI_max) from the nine models mentioned above (Tian and Fu, 2022; Dai et al., 2023). The NDVI_max data were maximum monthly NDVI for each year. Monthly NDVI data were obtained from the MOD13A3 NDVI (1 km × 1 km, 1 month) (Ma et al., 2022; Shen et al., 2022).

We combined the ggscatter function of the ggpubr package, the geom_smooth function of the ggplot2 package, the stat_poly_eq function of the ggpmisc package, and the ggarrange function of the ggpubr package to obtain the figures of the linear regression between observed and modeled SOC, TN, TP, C:N, C:P, and N:P, respectively. All statistical analyses were performed using R.4.1.2. In addition, the linear slope and R ² value, relative bias, and root-mean-square error (RMSE) were used to test the accuracy of the various models adopted in this study (Fu et al., 2011; Tian and Fu, 2022).

Overall, the ARad, AP and AT based on the RFM explained the most variation in SOC, TN, TP, C:N, C:P, and N:P, those based on the RRTM explained the second most variation, and those based on the MLRM explained the least variation under the non-grazing scenarios (Supplementary Tables S1, S3, S5). The ARad, AP, AT, and NDVI_max based on the RFM explained the most variation in SOC, TN, TP, C:N, C:P, and N:P, those based on the RRTM explained the second most variation, and those based on the MLRM explained the least variation under the free-grazing scenarios (Supplementary Tables S1, S3, S5). That is, the RFM had a greater ability to explain the six soil variables than the RRTM and MLRM, whereas the RRTM had a greater ability to explain the six soil variables than the MLRM. This finding supported some previous studies that demonstrated that the RFM had a greater ability to explain soil moisture and pH, plant species α-diversity, forage nutrition quality, and storage (Han et al., 2022; Tian and Fu, 2022; Dai et al., 2023; Wang and Fu, 2023). Moreover, the prediction accuracy of RFM was also higher than that of RRTM and MLRM (Tables 1, 2). This finding also supported some previous studies that demonstrated that the RFM had a higher prediction accuracy of soil moisture and pH, plant species α-diversity, forage nutrition quality, and storage than the RRTM and MLRM in grasslands on the Tibetan Plateau (Han et al., 2022; Tian and Fu, 2022; Dai et al., 2023; Wang and Fu, 2023). Therefore, the performance of the RFM was better than the RRTM and MLRM, at least for the grassland areas on the Xizang Plateau.

The R ² values under the non-grazing scenarios were not always lower than those under the free-grazing scenarios when we constructed the models of the six soil variables (Supplementary Tables S1, S3, S5). Moreover, the model accuracies under the non-grazing scenarios were not always lower than those under the free-grazing scenarios (Tables 1, 2). This finding implied that the increase in the number of independent variables cannot always improve the accuracy of the models. This finding was similar to several previous studies conducted on (Han et al., 2022; Tian and Fu, 2022; Dai et al., 2023; Wang and Fu, 2023) and outside the Qinghai-Xizang Plateau (Veloso et al., 2022). These previous studies found that the numbers of the input variables were not always positively related to the model accuracies of soil moisture and pH, plant species α-diversity, forage nutrition quality, and storage (Han et al., 2022; Tian and Fu, 2022; Dai et al., 2023; Wang and Fu, 2023). Moreover, a previous study found that the R ² value (<0.86) between the observed and modeled SOC pool in a previous study (Liu et al., 2023) was lower than that between the observed and modeled SOC in this study, but the input variables numbers (>=12) of the previous study (Liu et al., 2023) were much greater than those in this study. Therefore, an increase in the input variables of the model does not always lead to a better model of plant and soil variables, at least for the grasslands of the Xizang Plateau. Moreover, instead of focusing on increasing the number of input variables and how to obtain more accurate input variable data, it is better to focus on how to find the optimal model to better help solve related ecological and environmental problems.

Similar to some previous studies (Dai et al., 2023; Wang and Fu, 2023), the predicted accuracies of the six soil variables based on the RFM were dependent on soil depth. This finding may be due to the following mechanisms. First, soil moisture and pH can regulate their responses of SOC, TN, TP, C:N, C:P and N:P to external disturbance (e.g., warming) (Yu et al., 2019a). The predicted accuracies of the six soil variables should be related to those of soil moisture and soil pH. In addition, soil moisture and pH can generally change with soil depth. Second, solar radiation and precipitation first reach the surface of soils and then gradually reach the deep layer of soils through energy conduction or infiltration. The normalized difference vegetation index is calculated by the spectral information reflected by land surface, and the energy source of air temperature is mainly the long-wave radiation of the land surface. All of these may lead to closer relationships between the six topsoil variables and the four input variables than deeper soil variables.

The RFM and GBRM achieved a greater accuracy for the six soil variables than the other seven models. However, the tree numbers of the GBRM were much greater than those of the RFM for most cases (Supplementary Tables S1, S2). The model complexity can be generally positively correlated with the tree numbers, whereas the calculation speed can be generally negatively correlated with the tree numbers (Han et al., 2022; Tian and Fu, 2022). Therefore, the RFM should be the best model among the nine models when all things are considered in this study, which was in line with the hypothesis. This phenomenon was similar to some previous studies that demonstrated that the RFM performed better than the other models at modeling soil moisture and pH, plant species α-diversity, forage nutrition quality, and storage in grassland ecological systems on (Han et al., 2022; Tian and Fu, 2022; Dai et al., 2023; Wang and Fu, 2023) and outside the Qinghai-Xizang Plateau (Fernandez-Delgado et al., 2014; Nussbaum et al., 2018). All these findings further support the advantage of the RFM in modeling soil and plant variables in grasslands on the Qinghai-Xizang Plateau.

The prediction accuracy of the six variables based on the RFM was greater than those in some previous studies (Yang et al., 2008; Yang et al., 2009; Yang et al., 2010; Shangguan et al., 2013; Kou et al., 2019; Wang et al., 2020; Wang et al., 2021a; Liu et al., 2023). For example, the modeled TN could only explain 34%–41% of the variation in the observed TN, and the RMSEs and relative bias were 1.15–1.20 g kg⁻¹ and 49.17% in China, respectively (Zhou et al., 2020). Moreover, compared with some previous studies (Kou et al., 2019; Wang et al., 2021a; Liu et al., 2023), although the number of independent variables used in the models constructed in this study was small, the accuracy was not reduced. This phenomenon may be due to the introduction of more independent variables that may introduce new uncertainties, considering that any independent variable can have its own data quality. Therefore, the RFM of SOC, TN, TP, C:N, C:P and N:P at the three depths in this study can have relative accuracies and can be used for related studies (e.g., the spatial distribution of these six variables) in grassland areas on the Xizang Plateau.

It is well known that climate change and human activities together affect these soil variables (Yu et al., 2019a; Zhang and Fu, 2021; Zha et al., 2022), but their relative contributions are unclear. Many previous studies only modeled actual soil organic carbon, total nitrogen, total phosphorus, and their ratios (Yang et al., 2009; Kou et al., 2019), so it was difficult to disentangle the relative effects of climate change and human activity on these soil variables. However, in this study, the potential and actual random forest models of these soil variables were constructed respectively, so it was relatively easy to separate the relative contributions of climate change and human activities to these soil variables.

There are some uncertainties or limitations in this study. First, as the NDVI comes from optical remote sensing, it is often affected by the conditions of satellite itself (such as instrument failure), the angle between the satellite remote sensing and the surface, weather conditions, and land surface conditions. The Xizang Plateau became dimming, which coincided with the increase in precipitation during 2000–2020 (Huang and Fu, 2023). During the rainy season on the Xizang Plateau, it is often rainy or cloudy for several consecutive days, which may prevent the collection of near infrared and red band data of plants and the calculation of the NDVI. On the Xizang Plateau, grassland productivity tends to decline from southeast to northwest. The NDVI is easy to saturate in areas with high vegetation productivity, while it is affected by soil interference in areas with low vegetation productivity (Ma et al., 2022; Shen et al., 2022). Second, many grassland areas on the Xizang Plateau have steep slopes, and the soil samples were mainly collected in areas with relatively shallow slopes. This may affect the relationships between these soil variables and climate data or NDVI data. Third, although climate and NDVI data are based on the corresponding year, soil variables are not entirely derived from observations in the same year, and climate change and human activities may alter these soil variables (Zhang et al., 2015; Fu and Shen, 2017). This phenomenon may also interfere with the model accuracy. Fourth, although the climate data used in this study can be highly accurate, there is some uncertainty associated with them. Fifth, only random forest models with soil variables at depths of 0–10, 10–20, and 20–30 cm were built in this study; therefore, future studies may consider building optimal models of soil variables at greater depths. Sixth, the data in this study were all obtained from the Xizang Plateau; therefore, it may be necessary to further test the accuracy of the models constructed in this study when extrapolating them to other regions of the Qinghai-Xizang Plateau. We can consider additional soil sampling in other areas of the plateau to build the optimal model of soil variables in the entire Qinghai-Xizang Plateau region in the future. Finally, only random forest models of the contents of these soil variables were constructed in this study, and the densities of soil organic carbon, total nitrogen, and total phosphorus could better reflect the changes in these soil variables. Therefore, the construction of optimal models of soil organic carbon density, total nitrogen density, and total phosphorus density should be considered in the future.

This study is the first to model SOC, TN, TP, C:N, C:P, and N:P at depths of 0–10, 10–20, and 20–30 cm under non-grazing and free-grazing scenarios based on nine methods. We can build a database of soil carbon, nitrogen, and phosphorus and their ratios over the past few decades or for the future across the whole grasslands of the Xizang Plateau under fencing and free-grazing conditions using the RFM established in this study. This database can be used to solve some of ecological, environmental, and soil management problems in grassland ecosystems in the Xizang Plateau. For example, the relative effects of climate change and human activities on soil organic carbon during the past decades can be distinguished and quantified, thus providing data support for improving soil carbon sequestration and mitigating climate warming.