Assessing Soil Organic Carbon Stock Dynamics under Future Climate Change Scenarios in the Middle Qilian Mountains

: Soil organic carbon (SOC) its evaluations of the spatio-temporal dynamics of SOC stock under future climate change are crucial for the adaptive management of regional carbon sequestration. Here, we evaluated the dynamics of SOC stock to a 60 cm depth in the middle Qilian Mountains (1755–5051 m a.s.l.) by combining systematic measurements from 138 sampling sites with a machine learning model. Our results reveal that the combination of systematic measurements with the machine learning model allowed spatially explicit estimates of SOC change to be made. The average SOC stock in the middle Qilian Mountains was expected to decrease under future climate change, while the size and direction of SOC stock changes seemed to be elevation-dependent. Speciﬁcally, in comparison with the 2000s, the mean annual precipitation was projected to increase by 18.37, 19.80 and 30.80 mm, and the mean annual temperature was projected to increase by 1.9, 2.4 and 2.9 ◦ C under the Representative Concentration Pathway (RCP) 2.6 (low-emissions pathway), RCP4.5 (low-to-moderate-emissions pathway), and RCP8.5 (high-emissions pathway) scenarios by the 2050s, respectively. Accordingly, the area-weighted SOC stock and total storage for the whole study area were estimated to decrease by 0.43, 0.63 and 1.01 kg m –2 and 4.55, 6.66 and 10.62 Tg under the RCP2.6, RCP4.5 and RCP8.5 scenarios, respectively. In addition, the mid-elevation zones (3100–3900 m), especially the subalpine shrub-meadow Mollic Leptosols, were projected to experience the most intense carbon loss. However, the higher elevation zones (>3900 m), especially the alpine desert zone, were characterized by signiﬁcant carbon accumulation. As for the low-elevation zones (<2900 m), SOC was projected to be less varied under future climate change scenarios. the mid-elevation CO


Introduction
In alpine regions, elevation-dependent warming probably results in a higher temperature increase over higher elevation regions. Consequently, the patterns and dynamics of soil organic carbon (SOC) in alpine ecosystems are being markedly reshaped [1][2][3]. In addition, a large amount of SOC is stored in alpine ecosystems due to the cold and humid climate conditions, thus even small changes in temperature and precipitation could greatly affect SOC stock dynamics [4][5][6]. However, the size and direction of SOC dynamics under future climate change scenarios in alpine regions are still scarcely understood due to limited site-level measurements under such harsh environments and complex topography [7,8]. Thus, a comprehensive evaluation of the spatio-temporal patterns of SOC stocks in alpine regions may enable us to better understand the carbon-climate feedback [9][10][11]. This knowledge will help us to take dedicated adaptive measures to improve the soil carbon stock in alpine regions under future climate change [12,13].

Site-Level Measurements
We conducted soil sampling campaigns along the topography and vegetation type gradients in the study area from 2012 to 2015 during the growing seasons ( Figure 1). In total, 138 sampling sites were investigated, and 138, 133, 117 and 91 composite soil sam-  Overall, from summit to piedmont plain, the main vegetation types included alpine desert, alpine meadow, subalpine shrub, subalpine meadow, montane forest, montane shrub, montane meadow steppe, montane steppe, and montane desert steppe. Correspondingly, the main soil types were classified as alpine frigid desert soils, alpine meadow soils, subalpine shrubby meadow soils, subalpine meadow soils, gray cinnamon, chernozems, kastanozems, and sierozems according to the Chinese soil genetic classification [31], which were roughly referred to as Glacic Cryosols, Leptosols (Gelic), Eutric Leptosols, Mollic Leptosols, Greyzemic Umbrisols, Haplic Chernozems, Haplic Kastanozems, and Haplic Calcisols according to the World Reference Base for Soil Resources (WRB) [32]. For middleand high-elevation zones, soil depths generally ranged from 20 to 60 cm, and there were clear rock compositions in deeper soil layers. Considering the fundamental functions for regional ecological security, the Qilian Mountains were classified as a national nature reserve and a national park in 1980 and 2017, respectively [33]. Thus, the anthropogenic disturbances were relatively limited, and vegetation and soil dynamics were mainly controlled by climate change.

Site-Level Measurements
We conducted soil sampling campaigns along the topography and vegetation type gradients in the study area from 2012 to 2015 during the growing seasons ( Figure 1). In total, 138 sampling sites were investigated, and 138, 133, 117 and 91 composite soil samples were obtained for soil depths of 0-10, 10-20, 20-40 and 40-60 cm, respectively, and the topographic, geographic and vegetative variables for each sampling site were recorded. SOC concentration and bulk density were determined by the wet Walkley-Black method and the core ring method, respectively. Detailed field investigation and sample analysis processes can also be found in Zhu et al. [13]. The SOC stock was calculated using the following equation according to Ding et al. [5]: where SOCS is the SOC stock (kg m -2 ), BD i is the bulk density (g cm -3 ), D i is the soil thickness (cm), SOCC i is the SOC concentration (g kg -1 ), k is the number of layers divided in a soil profile, and C i is the volume percentage of the gravel > 2 mm in diameter at layer i.

Model Predictors
It has been suggested by previous studies that the topographic, climatic, pedological and vegetative variables are good predictors in modelling SOC spatial distribution across the Qilian Mountains [13,34]. Generally, the environmental covariates used for modelling SOC were highly self-correlated, and redundant inputs probably lowered model computational efficiency, especially at a high spatial resolution. Thus, only ten covariates, i.e., the elevation, aspect, slope gradient, topographic position index (TPI), longitude, latitude, temperature, precipitation, soil type, and normalized difference vegetation index (NDVI), were finally employed to develop the data-driven prediction model for SOC. Topography and locations were assumed to be stable at the time scale, whereas temperature, precipitation and NDVI varied on a decadal basis. Mean annual temperature and precipitation datasets from WordClim (https://worldclim.org/ (accessed on 20 February 2021)) were available at a 1 km resolution for the 2000s  and the 2050s (2041-2060). The future climate data of the 2050s were the averages of five global climate models (GCMs), i.e., the BCC-CSM1-1 [35], NorESM1-M [36,37], HadGEM2-AO [38], IPSL-CM5A-LR [39] and GISS-E2-R [40], and data under the RCP2.6 (low emission pathway), RCP4.5 (low-to-moderate emission pathway) and RCP8.5 (high emission pathway) scenarios were employed for SOC dynamics prediction. The elevation, aspect, slope and TPI were derived from the ASTGTM 30 m digital elevation model (DEM) using ArcGIS 10.2 (ESRI Inc., RedLands, CA, USA) and SAGA GIS (http://www.saga-gis.org/ (accessed on 13 February 2021)). The soil type data at a 1 km resolution were derived from the Harmonized World Soil Database (HWSD 1.2) product (https://www.fao.org/soils-portal/data-hub/en/ (accessed on 27 Novmber 2021)), and was reclassified according to the international standard soil classification WRB [32]. In addition, we obtained the vegetation type map at a 1:100,000 scale for the upper reaches of the Heihe River basin from the National Tibetan Plateau Data Center (http://www.tpdc.ac.cn/ (accessed on 23 February 2021)).

Data-Driven Model
We employed a machine learning technique, i.e., the random forest algorithm [41], to construct a data-driven model which could simulate the nonlinear relationships between site-level SOC stock and time-stable (e.g., elevation, aspect, slope, TPI, longitude, latitude) as well as time-varying covariates (e.g., temperature, precipitation and NDVI). During the training, we set the number of trees (ntree) to 1000 and the minimum size of terminal nodes (nodesize) at 5, according to previous studies focusing on SOC mapping in the middle Qilian Mountains with the random forest algorithm [13,42]. To minimize the uncertainty in extracting NDVI values from raster products for each of the sampling sites, we employed the averages of 30 m monthly NDVI layers (20 images in total) during the growing seasons (from May to September) of 2012-2015 (i.e., the soil sampling period). In addition, the temperature and precipitation values for each sampling site were not extracted from the raster products at a resolution of 1 km, and this was mainly because the climate in the Qilian Mountains was markedly shaped by elevation, and direct extraction from coarse resolution climate products may bring high levels of uncertainty to the simulation of carbon-environment relationships. We finally calculated mean annual precipitation and temperature at each sampling site according to the multiple linear regression equations expressed as functions of elevation, longitude, and latitude, as shown by Zhao et al. [43,44]. It should be noted that although Zhao et al. [43] also demonstrated that ordinary kriging (OK) yielded smaller prediction errors in growing-season temperature simulation than multiple linear regression did, we still employed the regression method to estimate the site-level annual temperature and precipitation. This was mainly because the OK method was more suitable for climatic variables interpolation in plains or regions with large spatial extents, while it was less suitable in mountainous regions where climates could vary sharply along the elevation gradient within a very short distance.
The 'randomForest' package was employed to conduct the data-driven model by the R language (R version 3.5.1, R Development Core Team, 2018). To further validate the accuracy of the random forest algorithm in modelling the non-linear relationships for site-level measurements, the ten-fold cross-validation procedure was conducted, and the coefficient of determination (R 2 ) and root mean square error (RMSE) were used to indicate the performance of the model. For details about the ten-fold cross-validation procedure and statistics calculation, see Zhu et al. [13]. It should be noted that soil depth in Glacis Cryosols was generally less than 20 cm, and thus the SOC stock in this soil was only modeled and projected with measurements from depths pf 0-20 cm. For other soils, we modeled SOC dynamics at a depth of 60 cm. The final SOC maps were actually produced by merging the 0-20 cm SOC stock in Glacis Cryosols and the 0-60 cm SOC stock in other soil types.

SOC Prediction
The data-driven model that simulated the non-linear relationships between SOC and environmental variables at a spatial scale was further used to project the SOC changes over time, and such an approach is called space-for-time substitution. However, the efficiency of this approach largely relies on the extents of environmental gradients of site-level measurements, i.e., the values of SOC and other environmental covariates should cover wide varying ranges. The huge elevation range (>3000 m) in the study area markedly shaped the water and heat conditions, and created diverse vegetation communities and soil types, thus making the Qilian Mountains an ideal region to conduct the space-fortime substitution procedure for SOC dynamics modelling. To create spatially explicit estimates of the SOC stock for both the 2000s (regarded as the baseline period) and the 2050s as affected by climate change, the data-driven model was rerun for the study area using the same topographic and geographic and pedological raster layers except for temperature, precipitation and NDVI data. The climate datasets for the 2000s and the 2050s were derived from the same source (i.e., the WorldCim datasets) so as to make the predictions comparable.
The soil type will most likely not change within 30 to 50 years, while the vegetation community and species composition may change significantly over time, especially for steppes and meadows. Thus, the vegetation type maps for the study area by the 2050s under the three RCP scenarios were also predicted. In this study, we developed predictive relationships between vegetation type map and climatic, topographic, and geographic covariates using a random forest model. The vegetation type was coded as 1 to 10 according to elevation before further modelling, which subsequently represented montane desert steppe, montane steppe, montane meadow steppe, montane shrub, montane forest, subalpine meadow, subalpine shrub, subalpine meadow, alpine desert and no vegetation. The predictive relationships were modeled based on all 10,500 grid cell attributes extracted from the vegetation type and other covariates (i.e., temperature, precipitation, topography and locations) maps at a resolution of 1 km × 1 km. Then, the vegetation type maps for the 2000s and the 2050s under the three RCP scenarios were subsequently modeled by the predictive relationships combined with covariates during corresponding periods. Similar procedures were also employed to estimate NDVI in 2050s. Specifically, we used NDVI averages of 1990-2000 for the 2000s, and calculated NDVI changes until the 2050s, which were then added to the NDVI spatial layer of the 2000s to obtain the NDVI for the 2050s. It should be noted that the vegetation type maps for different periods were also used as environmental covariates in projecting NDVI changes over time.
We also calculated the difference between the SOC stock of the 2000s and 2050s under different climate change scenarios for each 1 km × 1 km pixel, and further obtained the spatial distribution of SOC stock changes. In addition, we overlaid the vegetation type maps of both the 2000s and the 2050s and the soil type map over the SOC maps to summarize the means and standard deviations of SOC changes for different vegetation and soil types.

Future Climate and Vegetation Change Characteristics
The climate in the middle Qilian Mountains was expected to getting warmer and wetter from the 2000s to the 2050s (Figures 2 and 3). Specifically, the area-weighted temperature was −3.3 • C for the 2000s, while increased by 1.9 • C, 2.4 • C and 2.9 • C to the 2050s under the RCP2.6, RCP4.5 and RCP8.5 scenarios, and higher elevation regions seemed to be characterized by a slightly larger increase in comparison with the lower elevation zones, especially under the RCP8.5 scenario ( Figure 4). As for the precipitation, it was 344.70 mm in the 2000s, while it increased by 18.37, 19.80 and 30.80 mm to the 2050s under the three scenarios, and precipitation in higher elevation regions tended to experience greater increases ( Figure 5). The vegetation type distribution tended to be less varied over time ( Figure 6). The area of subalpine meadow, montane shrub and montane forest increased by 199.29, 40.46 and 3.47 km 2 , respectively, from the 2000s to the 2050s, while the area of alpine meadow and alpine desert shrank by 124.37 and 57.69 km 2 under the RCP2.6 scenario, respectively. By contrast, the spatial distribution of NDVI varied significantly from the 2000s to the 2050s (Figure 7), and NDVI was estimated to increase in high-elevation zones while decrease in mid-elevation zones (Table 1). Overall, the mean NDVI in the 2050s was 0.500, 0.507 and 0.492 under the RCP2.6, RCP4.5 and RCP8.5 scenarios, which was higher than that in the 2000s by 0.012, 0.019 and 0.004, respectively. creased by 199.29, 40.46 and 3.47 km 2 , respectively, from the 2000s to the 2050s, while the area of alpine meadow and alpine desert shrank by 124.37 and 57.69 km 2 under the RCP2.6 scenario, respectively. By contrast, the spatial distribution of NDVI varied significantly from the 2000s to the 2050s (Figure 7), and NDVI was estimated to increase in high-elevation zones while decrease in mid-elevation zones (Table 1). Overall, the mean NDVI in the 2050s was 0.500, 0.507 and 0.492 under the RCP2.6, RCP4.5 and RCP8.5 scenarios, which was higher than that in the 2000s by 0.012, 0.019 and 0.004, respectively.       Variation of mean annual temperature in the 2000s, and temperature change (difference for temperature between the 2050s and the 2000s) under the RCP2.6, RCP4.5, and RCP8.5 scenarios, along the elevation gradient for each 1 km × 1 km pixel. Note: The blue lines indicate the area-weighted averages of temperature or temperature difference between the 2050s and the 2000s for the whole study area. Each vertical gray line represents a 1 km × 1 km pixel over the study area, and the height of the gray line indicates the temperature value at that pixel. Most parts of the study area were located in regions with elevations ranging from 3000 to 3800 m, and thus, gray lines are much denser in mid-elevation zones, while sparser in high-and low-elevation zones.

SOC Prediction under Different Climate Change Scenarios
The data-driven model based on the random forest algorithm could accurately simulate the non-linear relationships between SOC stock and temperature, precipitation, vegetation, terrain attributes, soil type and geographic positions, with a slope of 0.77 and an intercept of 5.64 kg m -2 between the measured and predicted values for the test datasets ( Figure 8). The ten-fold cross-validation results further show that the R 2 and RMSE values of the model were 0.78 and 5.10 kg m −2 , respectively ( Table 2). The estimated SOC stock maps demonstrated large heterogeneity of SOC in alpine regions, ranging from 0 to 38.97 kg m −2 in the 2000s, and 0 to 38.77 kg m −2 , 0 to 38.62 kg m −2 and 0 to 36.27 kg m −2 in the 2050s for the RCP2.6, RCP4.8 and RCP8.5 scenarios, respectively ( Figure 9). In addition, the spatial patterns of SOC stock as predicted by the model were closely related to elevation in both the 2000s and the 2050s, i.e., the SOC stock in mid-elevation zones was significantly higher than that in lower and higher elevation zones ( Figure 9). Variation of mean annual temperature in the 2000s, and temperature change (difference for temperature between the 2050s and the 2000s) under the RCP2.6, RCP4.5, and RCP8.5 scenarios, along the elevation gradient for each 1 km × 1 km pixel. Note: The blue lines indicate the area-weighted averages of temperature or temperature difference between the 2050s and the 2000s for the whole study area. Each vertical gray line represents a 1 km × 1 km pixel over the study area, and the height of the gray line indicates the temperature value at that pixel. Most parts of the study area were located in regions with elevations ranging from 3000 to 3800 m, and thus, gray lines are much denser in mid-elevation zones, while sparser in high-and low-elevation zones.

SOC Prediction under Different Climate Change Scenarios
The data-driven model based on the random forest algorithm could accurately simulate the non-linear relationships between SOC stock and temperature, precipitation, vegetation, terrain attributes, soil type and geographic positions, with a slope of 0.77 and an intercept of 5.64 kg m -2 between the measured and predicted values for the test datasets ( Figure 8). The ten-fold cross-validation results further show that the R 2 and RMSE values of the model were 0.78 and 5.10 kg m −2 , respectively ( Table 2). The estimated SOC stock maps demonstrated large heterogeneity of SOC in alpine regions, ranging from 0 to 38.97 kg m −2 in the 2000s, and 0 to 38.77 kg m −2 , 0 to 38.62 kg m −2 and 0 to 36.27 kg m −2 in the 2050s for the RCP2.6, RCP4.8 and RCP8.5 scenarios, respectively ( Figure 9). In addition, the spatial patterns of SOC stock as predicted by the model were closely related to elevation in both the 2000s and the 2050s, i.e., the SOC stock in mid-elevation zones was significantly higher than that in lower and higher elevation zones ( Figure 9).

Elevation-Dependent SOC Dynamics in the Qilian Mountains
The area-weighted SOC stock was predicted to decrease significantly from the 2000s to the 2050s, with a loss of 0.43, 0.63 and 1.01 kg m -2 carbon (i.e., 4.55, 6.66 and 10.62 Tg for SOC storage over the basin) under the RCP2.6, RCP4.5 and RCP8.5 scenarios, respectively (Tables 3 and 4). However, the direction of SOC changes varied a lot at the basin scale, i.e., SOC losses mainly occurred at the mid-west and east parts of the basins ( Figure  10). Further analysis showed that the size and direction of SOC changes seemed to be elevation-dependent ( Figure 11). As for the RCP2.6 scenario, the SOC stock increased in regions with an elevation <3100 m or >3900, while it decreased at 3100-3900 m and varied less below 2900 m. Similar distribution patterns for SOC dynamics along the elevation gradient were also found under the RCP4.5 and RCP8.5 scenarios.
SOC changes also varied significantly with vegetation and soil types. The subalpine meadow was estimated to experience a decrease of 2.85, 3.24 and 3.64 kg m -2 for SOC stock (i.e., 5.03, 5.76 and 6.60 Tg for SOC storage) under the RCP2.6, RCP4.5 and RCP8.5 scenarios, respectively (Tables 3 and 4). However, the alpine desert tended to significantly accumulate organic carbon in soils, reaching 0.76, 0.99 and 1.06 kg m -2 for SOC stock, and 1.32, 1.72 and 1.86 Tg for SOC storage. As for soil types, the SOC stock in the Eutric Leptosols was expected to be characterized by the largest decrease, reaching 1.90, 1.96 and 2.38 kg m -2 for the RCP2.6, RCP4.5 and RCP8.5 scenarios, respectively (Table 5). Nevertheless, the largest total amount of carbon loss was contributed by the Mollic Leptosols, which were projected to lose 4.15, 5.33 and 6.71 Tg of carbon under the three scenarios by the 2050s. By contrast, the Glacic Cryosols were estimated to significantly accumulate organic carbon, reaching 0.83, 0.95 and 0.84 kg m -2 for SOC stock, and 1.60, 1.83 and 1.62 Tg for SOC storage by the 2050s under the three scenarios (Table 5).

Elevation-Dependent SOC Dynamics in the Qilian Mountains
The area-weighted SOC stock was predicted to decrease significantly from the 2000s to the 2050s, with a loss of 0.43, 0.63 and 1.01 kg m -2 carbon (i.e., 4.55, 6.66 and 10.62 Tg for SOC storage over the basin) under the RCP2.6, RCP4.5 and RCP8.5 scenarios, respectively (Tables 3 and 4). However, the direction of SOC changes varied a lot at the basin scale, i.e., SOC losses mainly occurred at the mid-west and east parts of the basins ( Figure 10). Further analysis showed that the size and direction of SOC changes seemed to be elevationdependent ( Figure 11). As for the RCP2.6 scenario, the SOC stock increased in regions with an elevation <3100 m or >3900, while it decreased at 3100-3900 m and varied less below 2900 m. Similar distribution patterns for SOC dynamics along the elevation gradient were also found under the RCP4.5 and RCP8.5 scenarios.       SOC changes also varied significantly with vegetation and soil types. The subalpine meadow was estimated to experience a decrease of 2.85, 3.24 and 3.64 kg m -2 for SOC stock (i.e., 5.03, 5.76 and 6.60 Tg for SOC storage) under the RCP2.6, RCP4.5 and RCP8.5 scenarios, respectively (Tables 3 and 4). However, the alpine desert tended to significantly accumulate organic carbon in soils, reaching 0.76, 0.99 and 1.06 kg m -2 for SOC stock, and 1.32, 1.72 and 1.86 Tg for SOC storage. As for soil types, the SOC stock in the Eutric Leptosols was expected to be characterized by the largest decrease, reaching 1.90, 1.96 and 2.38 kg m -2 for the RCP2.6, RCP4.5 and RCP8.5 scenarios, respectively (Table 5). Nevertheless, the largest total amount of carbon loss was contributed by the Mollic Leptosols, which were projected to lose 4.15, 5.33 and 6.71 Tg of carbon under the three scenarios by the 2050s. By contrast, the Glacic Cryosols were estimated to significantly accumulate organic carbon, reaching 0.83, 0.95 and 0.84 kg m -2 for SOC stock, and 1.60, 1.83 and 1.62 Tg for SOC storage by the 2050s under the three scenarios (Table 5).

SOC Prediction
Considering the ideal performance of the random forest algorithm in capturing the non-linear relationships between the SOC and other environmental covariates, we employed this algorithm to construct the data-driven model for future SOC dynamics projection based on the space-for-time substitution method [45,46]. Similar method was also employed by Gray et al. [14] and Adhikari et al. [24]. However, it should be noted that although land use change was also an important variable in modelling SOC dynamics in previous studies [16,21,47], we did not include it in the model for this study. This was mainly because the anthropogenic disturbances were relatively limited after the National Nature Reserve of the Qilian Mountains was established in 1988, and the human-induced land use change was not intense [48]. Furthermore, the vegetation activity in alpine regions was mainly affected by precipitation and temperature [27], and its effect on SOC dynamics was also considered in the model by using the dynamic NDVI datasets as model inputs in SOC projection. SOC changes over time may also be affected by other environmental factors such as nitrogen deposition, climate extremes, CO 2 enrichment, and ecological restoration practices [22]. However, in regions such as the middle Qilian Mountains characterized by a semiarid alpine climate, water and heat conditions may be the dominant factors affecting ecosystem productivity over time, in comparison with factors such as nitrogen deposition, climate extremes, and CO 2 enrichment. As for ecological restoration practices, the complex topography and limited accessibility made it impossible to conduct large-scale reforestation or oversowing in high-elevation zones. Hence, current ecological restoration practices in the middle Qilian Mountains mainly seek to maintain the authenticity of zonal ecosystems through excluding human disturbances [8,33]. In addition, the objective of this study was to determine the potential responses of the SOC temporal dynamics to future climate change scenarios. Although omitting these factors in modelling may lend uncertainty to these predictions to some degree [23,24,49], we believed that our results about SOC changes in response to future climate change were instructive for future carbon management in the Qilian Mountains. Overall, the sharp changes in precipitation and temperature and limited anthropogenic disturbances in the Qilian Mountains enabled us to use the space-for-time substitution method to model SOC dynamics over time with the carbon-environment relationships derived at the spatial scale. In this study, the ten-fold cross-validation presented a quite high R 2 value of 0.74 for the test dataset, which was slightly higher than that of previous studies focusing on modelling the spatio-temporal patterns of SOC in nearby regions [5,22,34]. Thus, our results indicate that the combination of systematic measurements with the machine learning model allowed spatially explicit estimates of SOC changes to be made.

Elevation-Dependent SOC Dynamics
Our model suggests that the SOC stock in the middle Qilian Mountains is expected to experience overall losses under the three climate change scenarios. Our results coincide with those from previous studies [24,50]. However, as we hypothesized, we found that the size and direction of SOC changes seemed to be elevation-dependent, i.e., the SOC stock decreased significantly in mid-elevation zones, while it increased in high-elevation zones and varied less in low-elevation zones in the study area. Similar results were also demonstrated by Chen et al. [2], who found that soils in the Yunnan province experienced intense carbon losses at 3000-3500 m while significant accumulation above 4500 m during the period of 1986-2015. In addition, Ding et al. [51] found that soils of alpine grasslands of the Tibetan Plateau (>4000 m) were also characterized by significant carbon accumulation from the 2000s to the 2010s. It was widely acknowledged that the increase in temperature could stimulate the activity of soil microorganisms, thus accelerating the loss rate of carbon from the soil into the atmosphere [20,52,53]. However, for high-elevation zones, the increase in temperature could also enhance the accumulated temperature, which prolongs the growing season and further stimulates the vegetation productivity. These effects finally lead to increased carbon inputs into soils from both aboveground biomass and roots. As a result, the warming-induced enhancement in productivity may offset the potential decrease in SOC as affected by increased microbial activities in high-elevation zones. Furthermore, considering the arid climate conditions in the middle Qilian Mountains, the increase in precipitation under the three climate change scenarios could significantly enhance the ecosystem net primary production, which also contributes to carbon accumulation in alpine soils.
The mid-elevation zones, which were mainly occupied by subalpine shrubs, subalpine meadows, and montane forests, were predicted to experience the most intense carbon losses under future climate change scenarios. Presumably, this was associated with the high initial standing SOC stock in the three vegetation types, which were larger than 24 kg m −2 (Table 3). Crowther et al. [54] demonstrated that the effect of warming on SOC losses depend on the initial standing carbon stock and considerable carbon losses occurred in soils with high carbon stock, based on assembling data from field experiments in the Northern Hemisphere. In this study, due to the very high SOC values of alpine ecosystems, increased temperature will strongly enhance the activity of microorganisms, leading to accelerated SOC decomposition and intense carbon losses. In addition, as the water and heat conditions in the mid-elevation zones were relatively favorable, both vegetation productivity and SOC stock were at a very high level, and a large amount of humus was accumulated in topsoils. In this case, the effects of precipitation change on vegetation productivity were less significant in comparison with those in low-elevation zones, where the vegetation growth was strongly affected by water conditions. As a result, the positive effects of wetting were much trivial in comparison with those of warming in the mid-elevation zones. Our results are similar to the findings of Prietzel et al. [50] in the German Alps, who demonstrated that warming significantly accelerated the release of carbon from soils with a high initial carbon content, especially for forest soils. In this study, the low-elevation zones were characterized by a relatively high temperature and limited precipitation, and thus water availability was the dominant environmental factor affecting vegetation productivity. The positive effects of precipitation increase may offset the negative effects induced by warming, and thus the SOC stock changes in low-elevation zones were relatively trivial.
The SOC spatio-temporal dynamics were hypothesized to be closely related to soil types. As an important pedological factor, the soil type was generally considered as a predictor for modelling SOC in previous studies [14,45]. In this study, we found that the inclusion of soil types could improve the R 2 values of random forest models by 4% and reduce the RMSE by 0.40 kg m −2 , suggesting an important role of soil types in shaping the spatio-temporal pattern of SOC. Considering the thin soil layers in alpine regions, we separately calculated the SOC stock at 0-20 cm for Glacic Cryosols, while calculating it at 0-60 cm for other soil types. Such efforts could avoid the overestimation of SOC stocks in Glacic Cryosols, and make the SOC projections in alpine regions more reasonable. In addition, we found that the size and direction of SOC changes under the three climate change scenarios differed by soil type. Specifically, the SOC stock tended to decrease significantly in Leptosols, especially in Eutric Leptosols and Mollic Leptosols, while it increased significantly in Glacic Cryosols. This was mainly because the Leptosols in the study area were located at 3000-3900 m, where the standing carbon stocks were very high (>24 kg m −2 ). In addition, the Leptosols were generally characterized by a thin soil layer, and thus topsoil SOC tended to be more prone to experiencing rapid loss due to climate warming. By contrast, the Glacic Cryosols were generally located in regions close to glaciers, where the heat conditions were the main drivers of vegetation activity. Climate warming could accelerate glaciers retreat, and carbon in Glacic Cryosols may significantly accumulate due to enhanced vegetation activity and litter input. The SOC stocks in Haplic Chernozems, Haplic Kastanozems and Haplic Calcisols were less varied, and this was probably associated with the opposite effects of climate warming and wetting in low-elevation zones, where the precipitation was generally less than 300 mm.

Implications for SOC Management under Climate Change
The elevation-dependent SOC dynamics implied that the carbon management strategy in the Qilian Mountains should vary with elevation zones. Soil carbon stocks in midelevation zones (3100-3900 m), especially in the subalpine shrub-meadow zone, should be given first priority in terms of carbon management under future climate change conditions. The adaptive management practices for carbon pools should mainly aim at effectively reducing carbon losses from soils into the atmosphere. In addition, considering the potential carbon accumulation in high-elevation zones in the future, the adaptive management practices should focus on lowering the intensity of anthropogenic disturbances such as coal mining, overgrazing, and illegal construction, which may reduce the carbon sequestration function of undisturbed alpine soils [33]. SOC in the low-elevation zones was estimated to be less varied under future climate change conditions. Hence, to enhance the carbon stock in the low-elevation zones, additional ecological construction projects instead of simple enclosure or conservation, should be implemented in the future. As the initial standing carbon stocks at low elevations were relatively low due to limited vegetation activity induced by less precipitation, soils here may have large carbon sequestration potential if the vegetation could be well irrigated by additional water [28].
The specific management practices should also differ by elevation zone. For the mid-elevation zones, the rapid losses of SOC were mainly related to the warming-induced increase in soil temperature, and thus measures that could hinder soil temperature increase should be given priority. A previous study showed that the shelter effect of the shrub canopy could significantly lower soil temperature and facilitate SOC accumulation in comparison with that of grasslands under similar topographic conditions in the Qilian Mountains [28]. Therefore, native shrub species such as Caragana jubata, Salix gilashanica and Potentilla fruticosa could be selected for afforestation in grasslands where the climatic and topographic conditions were also suitable for shrub growth. In addition, the number of large livestock should also be further reduced in the Qilian Mountains. The feeding behavior of large livestock such as yaks causes great damage to shrub branches and causes large areas of subalpine shrubs to wither ( Figure 12). As a result, the shelter effects of the shrub canopy were weakened, and the high initial standing carbon stocks of subalpine shrubs became more vulnerable under climate warming. In addition, the restoration of degraded grasslands as well as returning cultivated lands to grasslands could also contribute to the carbon sequestration in mid-elevation zones. For the high-elevation zones, SOC dynamics were mainly shaped by climate change, thus conservation may be a better choice considering the harsh climates and limited accessibility in alpine regions. Carbon sequestration in low-elevation zones was essentially an issue related to the tradeoff of water resources. In other words, an effective balance between the agricultural and ecological water demand is fundamental for future carbon management. In addition, considering the lower amount of total water resources under such arid climates, on the one hand, unconventional water resources such as seasonal floods in gullies should be utilized through a series of water conservancy projects, and on the other hand, high-efficiency water-saving measures such as drip irrigation should also be taken into account during ecological construction projects in low-elevation zones in the Qilian Mountains.

Conclusions
The area-weighted SOC stock in the middle Qilian Mountains was more than 22 kg m -2 , and was predicted to decrease significantly from the 2000s to the 2050s with a loss of 0.43, 0.63 and 1.01 kg m −2 under the RCP2.6, RCP4.5 and RCP8.5 scenarios, respectively. As for the whole basin, the total carbon storage loss was 4.55, 6.66 and 10.62 Tg, respectively. However, SOC dynamics were estimated to be elevation-dependent, i.e., the SOC stock increased in regions with an elevation <3100 m or >3900, while it decreased at 3100-3900 m and varied less below 2900 m. SOC changes also varied significantly with vegetation and soil types. The subalpine meadow was projected to experience the most intense carbon loss, reaching to 2.85, 3.24 and 3.64 kg m −2 , while the alpine desert tended to significantly accumulate organic carbon in soils, reaching to 0.76, 0.99 and 1.06 kg m −2 under the three scenarios. The largest total amount of carbon loss was contributed by the Mollic Figure 12. Effect of different domestic animals on vegetation cover in the subalpine shrub. Note: (a,b) are a full and close shot of the subalpine shrub-meadow zone, respectively. The ridge in (a) is covered by fence with the left side mainly being grazed by yak and the right side by goats. The left side of the fence in (b) was mainly grazed by goats, and the right side was mainly grazed by yaks. The photos were taken at Toutan (3400 m a.s.l.) in the subalpine shrub-meadow zone of the middle Qilian Mountains.

Conclusions
The area-weighted SOC stock in the middle Qilian Mountains was more than 22 kg m -2 , and was predicted to decrease significantly from the 2000s to the 2050s with a loss of 0.43, 0.63 and 1.01 kg m −2 under the RCP2.6, RCP4.5 and RCP8.5 scenarios, respectively. As for the whole basin, the total carbon storage loss was 4.55, 6.66 and 10.62 Tg, respectively. However, SOC dynamics were estimated to be elevation-dependent, i.e., the SOC stock increased in regions with an elevation <3100 m or >3900, while it decreased at 3100-3900 m and varied less below 2900 m. SOC changes also varied significantly with vegetation and soil types. The subalpine meadow was projected to experience the most intense carbon loss, reaching to 2.85, 3.24 and 3.64 kg m −2 , while the alpine desert tended to significantly accumulate organic carbon in soils, reaching to 0.76, 0.99 and 1.06 kg m −2 under the three scenarios. The largest total amount of carbon loss was contributed by the Mollic Leptosols, which were projected to lose 4.15, 5.33 and 6.71 Tg carbon, while the Glacic Cryosols could possibly sequester 1.60, 1.83 and 1.62 Tg carbon under the three scenarios by the 2050s.
Thus, SOC management practices should also differ by elevation zone and vegetation and soil type. The mid-elevation zones, especially the subalpine shrub-meadow and Mollic Leptosols, should be given priority in terms of reducing CO 2 emissions, and the main goal for this region should focus on slowing down the carbon loss rates. Native shrub species such as Caragana jubata, Salix gilashanica and Potentilla fruticosa could be selected for afforestation in grasslands where the climatic and topographic conditions were also suitable for shrub growths. The high-elevation zones, especially the alpine desert and Glacic Cryosols, were expected to experience significant carbon accumulation under future climate change scenarios. Thus, vegetation and soils here should be protected from anthropogenic disturbances to keep the authenticity of the zonal ecosystems. SOC sequestration in low-elevation zones should rely on additional ecological construction projects. Unconventional water resources such as seasonal floods in gullies should be utilized through a series of water conservancy projects, and high-efficiency water-saving measures such as drip irrigation should also be taken into account during ecological construction projects in low-elevation zones in the Qilian Mountains.
Overall, our results highlight the important roles of elevation, vegetation and soil types in soil carbon management under future climate change conditions in the Qilian Mountains.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author for research purposes.

Conflicts of Interest:
The authors declare no conflict of interest.