The grasslands of Central Asia are located between 30°-50°N and 40°-105°E (Fig. 1), and belong to typical temperate continental climate. The significant variation in altitude within the area (-173-7347 m a.s.l.) results in distinct vertical zonation of vegetation with grassland types ranging from desert grassland (DG), temperate grassland (TG), and forest meadow (FM) (Han et al., 2014). The annual mean temperature in this area is recorded at 3.80℃, while the mean annual precipitation is 1949.0 mm, with June to September contributing to approximately 60.6% of the annual total precipitation. The mountainous areas in the southeastern Central Asia are much wetter with an average of 500.0 mm/a (Ta et al., 2018). The mean temperature is 15.00℃ in the low latitude areas and below 0.00℃ in the high latitude areas (Mitchell and Jones, 2005). Datasets revealed a significant regional increase in surface air temperature ranging from 0.36℃ to 0.42°C over the past 33 a (Hu et al., 2014).

The DNDC model was originally developed for the estimation of C sequestration and N₂O emissions in agricultural ecosystems. Through long-term application (Li et al., 1992), researchers have used this model to simulate biogeochemical C cycle for almost all terrestrial ecosystems (i.e., farmlands, forest lands, wetlands, and grasslands) (Smakgahn et al., 2009; Li et al., 2012; Katayanagi et al., 2013; Guest et al., 2017; Wu et al., 2018). The model comprises six interconnected sub-models, encompassing soil and climate processes, vegetation growth dynamics, decomposition mechanisms, as well as nitrification, denitrification, and fermentation processes. It is primarily driven by four fundamental ecological factors, i.e., climate conditions, soil properties, vegetation characteristics, and management practices (Li et al., 1992).

Detailed information on soil parameters, including clay fraction, textural class, drainage rate, bulk density, organic C content, pH value, and cation exchange capacity are extracted from the corresponding grid cell in the Harmonized World Soil Database (HWSD) (https://gaez.fao.org/pages/hwsd) (FAO/IIASA/ISRIC/ISS-CAS/JRC, 2009). Elevation data were derived from the WorldClim data (http://www.worldclim. org/download) (Hijmans et al., 2005). The spatial data were subsequently smoothed to a resolution of 40 km. N deposition data derived from a mosaic Asian anthropogenic emission inventory under the international collaboration framework of the Model Inter-Comparison Study for Asia (MICS-Asia) (Li et al., 2017). Three different long-term gridded atmospheric reanalysis datasets were used including Medium-Range Weather Forecasts (ECMWF) Re-Analysis-Interim (ERA-Interim), Modern-Era Retrospective analysis for Research and Applications (MERRA), and Climate Forecast System Reanalysis (CFSR). ERA-Interim is a global atmospheric reanalysis produced by the European Centre for ECMWF with a spatial resolution of about 80 km (https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-interim). The time frame covered by the dataset starts from 1 January 1979 and continuously updated in near-real time. Detailed description of ERA-Interim data archive is given in Dee et al. (2011). MERRA was generated utilizing the Goddard Earth Observing System Data Assimilation System v.5.0 (GEOS-5), comprising of the GEOS-5 atmospheric model and the Grid-point Statistical Interpolation (GSI) analysis system (Rienecker et al., 2011). The MERRA data set is available from 1979 and onwards with a resolution of 0.50°×0.65° (https://disc.gsfc.nasa.gov/information/mission-project?title=MERRA-2). The CFSR is designed and executed as a globally comprehensive and high-resolution coupled atmosphere-ocean-land surface-sea ice system, aiming to provide the most accurate estimation of the state of these interconnected domains during the study period (https://rda.ucar.edu/datasets/ds093.0/). The CFSR global atmosphere data exhibit a spatial resolution of approximately 38 km, and have been available since 1979 (Saha et al., 2014). The datasets exhibit variations in terms of their spatial and temporal resolution, available time period, as well as the methodology employed for their derivation. To facilitate a direct comparison among simulations, we standardized the datasets to a uniform spatial (40 km×40 km) and temporal scale (daily resolution). The time period of 1979-2014 was adopted as the common interval for air temperature, precipitation, and shortwave radiation serving as shared climate variables across all datasets.

After parameterizing the DNDC model in arid grasslands, a total of 48 net primary productivity (NPP) and 6 vegetation C were used for model validation. Due to the limited availability of validation data in Central Asia, we collected data from grasslands that exhibit similar conditions to our own grassland ecosystem from other literatures (Anwar et al., 2006; Zhao et al., 2006, 2007; Zhang et al., 2008; Fan et al., 2009; Toderich et al., 2009; Yan, 2009; Yang et al., 2010; Zhang et al., 2012). We used R² to verify the accuracy of DNDC model, which can reflect the consistency between simulated and observed values.

For this study, the effects of N deposition, climate change, and grazing on C and water dynamics were evaluated through three distinct simulations. We designed two groups of management practice scenarios for this purpose with and without N deposition. Prior to conducting various simulations, we utilized the mean values of temperature, precipitation, and N deposition rates spanning from 1979 to 2014 in conjunction with other model initial datasets for model spin-up until C and N pools achieved equilibrium.

The leaf temperature data were utilized to calculate the WSI for both the controlled environment and field components of this study. To determine WSI, we employed the empirical method proposed by Gardner et al. (1992), consisting of the following equation:

where dT is the difference between leaf temperature and ambient air temperature (°C); and the subscripts m, LL, and UL represent measured, lower limit, and upper limit temperature differences, respectively.

We conducted a comparative analysis between site simulations and observations to assess the accuracy of NPP and vegetation C estimates, which serve as indicators of ecosystem C dynamics. The R² values between observed and simulated NPP and vegetation C were 0.83 and 0.99, respectively (Fig. 2). The agreement between observed and simulated values was generally satisfactory, indicating that the model is suitable for simulating the spatial and temporal variations of C dynamics in Central Asia.

Figure 3 shows spatial distribution of simulated annual average vegetation C during the period 1979-2014. Total vegetation C in Central Asia was 0.35 (±0.09) Pg C/a for the whole area. The highest vegetation productivity was distributed in FM where average vegetation C was 165.30 (±47.36) g C/(m²•a), while DG and TG had a lower NPP. WSI with values ranging from 0.00 (maximum stress) to 1.00 (no stress), and the average WSI was 0.20 (±0.02) in Central Asia. WSI was also high in FM (0.34) and was the lowest in DG with mean WSI of 0.09 (±0.01). The temporal variations of vegetation C during the period 1979-2014 are shown in Figure 4. The simulated vegetation C increased by 3.27g C/(m²•a) (R²=0.83; Fig. 4a), while WSI showed no upward or downward trend.

Figures 5 and 6 show N influence on vegetation C difference and WSI variation during the period 1979-2014. During the past 36 a, increasing N deposition led to an increase in vegetation C of 65.56 (±83.03) Tg C at regional scale. As simulated in our experiments, the areas in which N deposition had a large effect on vegetation C were mainly located in DG (Figs. 5a-c and 6a). In these areas, N deposition increased vegetation C by more than 1.41 (±0.98) g C/(m²•a). N deposition also evidently increased vegetation C in FM by 0.68 (±0.61) g C/(m²•a), while decreasing trend was found in TG. Changes in N deposition led to vegetation C decreases of more than 31.83 (±75.83) Tg C in TG during the period 1979-2014. During the past 36 a, N deposition slightly decreased WSI in Central Asia. In addition, N deposition-induced decline in WSI was apparent in TG than in other grasslands. N deposition almost had no effect on WSI in the enriched DG (Figs. 5d-f and 6b).

We compared the effects of changing climate (without N deposition effect) on modeled C dynamics and WSI (Fig. 7) during the period 1979-2014. All results simulated by three reanalysis data showed a significant increasing trend in vegetation C, ranging from 2.10 (±0.41) g C/(m²•a) for CFSR to 5.04 (±0.36) g C/(m²•a) for ERA (Fig. 7a). Thus, our simulations suggested that climate change was the major driver for increased terrestrial productivity compared with N deposition. Meanwhile, no significant trend in WSI for the same period was detected.

Three types of grassland in this study showed the same trends but different extent in C and water dynamics over the last three decades. All the types of grassland showed an increasing trend of vegetation C (Fig. 7c), however, FM showed more significant increasing value of 5.31 (±0.25) g C/(m²•a) (R²=0.92), while the increasing values were 2.94 (±0.49) g C/(m²•a) (R²=0.49) and 1.54 (±0.02) g C/(m²•a) (R²=0.57) for TG and DG, respectively. No significant trend in WSI for the same period was detected.

Atmospheric N deposition can significantly influence the C and N cycles of terrestrial ecosystems, thereby impacting their structure and functioning. In recent years, numerous experimental studies have been conducted globally to simulate N deposition effects in various ecosystems (Stevens et al., 2011; Kinugasa et al., 2012; Zhang et al., 2023). Previous research has highlighted the crucial role of N as a limiting factor for grassland growth (Li et al., 2015b; Stevens et al., 2015). However, it is worth noting that in arid and semi-arid areas where vegetation growth is primarily constrained by precipitation, the contribution of N has often been overlooked. Lu et al. (2016) illustrated in their research that increased N deposition raised NPP by 9.60 g C/(m²•a) on average, accounting for around 92.20 Tg C/a of the national total. In Europe, the average NPP under N deposition was 5.00-75.00 g C/g N (de Vries et al., 2009). In China, the average NPP under N deposition in grasslands was 25.00 g C/g N in early 21^st century (Lu et al., 2012). Our study indicated that N deposition led to an increase of 65.56 (±83.03) Tg C in the grasslands in Central Asia. The substantial increase in vegetation C can be attributed to enhanced photosynthesis resulting from the rise in plant N with N addition (Kinugasa et al., 2012). Zhu et al. (2021) also observed that N input in arid and semi-arid grasslands could stimulate root growth for improved nutrient and water uptake. Considering the prevalence of widespread N limitation in Central Asian grasslands as indicated by experimental and monitoring studies, our modeled C sequestration under N deposition further supported the hypothesis of N limitation in this area, although with moderate magnitude diminishing over time.

In terms of N deposition-induced C sequestration, our simulations and other field investigations (Wang et al., 2021) suggest that grasslands serve as a stronger C sink. While there are limited indications of N saturation in our simulations, determining the threshold for N saturation is crucial across different types of grassland ecosystems.

The results showed different trends in C and water dynamics under N deposition over the past 36 a. Different climate zones and soil characteristics associated with each type of grassland affect their trends. Our study indicated that different types of grassland responded variably to N deposition. In detail, N deposition had positive effects on DG productivity, which was consistent with the findings of Li (2021). DG was modeled to exhibit a higher degree of N limitation on vegetation C compared with other types of grassland, which is consistent with empirical observations and modeling results. Kinugasa et al. (2012) imply that increased N deposition can enhance grassland recovery after a drought even in arid areas like the Mongolian Steppe. This greater productivity might be attributed to plants enhancing photorespiration under drought conditions, which stimulates malate production in chloroplasts and generates reductants for nitrate assimilation, making them particularly adept at utilizing soil nitrate as a source of N (Eisenhut et al., 2019).

Notably, effects of N dynamics were negative in TG, whether it was compensated for by other process settings in the model, or reached N saturation was uncertain. The negative effects of N deposition on species richness in temperate grasslands in China have been widely observed (Fang et al., 2012), as well as decreased belowground biomass with N addition. One probably explanation is that N addition increased plant growth, which will increase the rate of water loss from plants and result in the drier condition, causing drought in the long term.

Water and N are two primary limiting factors across many ecosystems (Kamran et al., 2022). It has been suggested that applying N should only be done in areas with sufficient plant-available water since its uptake depends on soil and plant water status (Tilling et al., 2007). In numerous studies, researchers have demonstrated that adequate N nutrition can enhance plant drought resistance and effectively improve water relations when growing under dry soil conditions (Dang et al., 2006; Adamtey et al., 2010). Lin et al. (2012) found no significant effect of N fertilization on crop WSI under different soil drought treatments, which aligns with our findings. In the field of agriculture and forestry, WSI is extensively utilized for irrigation scheduling and N fertilization (Emekli et al., 2007), yet, our finding suggested that the calculation of N and water should be different, depending on different ecosystems.

In this study, we used three sets of climate data to partially correct the uncertainty of parameters, however, due to the time constraints on the data, our analysis is limited to the situation prior to 2015. Simulations driven by different climate datasets yielded different NPP results, when compared with field observation, and simulation ability of different meteorological datasets varied in different areas (Smakgahn et al., 2009; Gaillard et al., 2018). The model structure represents another significant source of uncertainty. We have not considered grasslands reclaimed for farmlands and other land use and land cover (LULC) changes, which have been estimated to result in substantial C emissions (Chang et al., 2022; Feng et al., 2023). Especially in TG, the anthropogenic influence was more obvious. Interactions among N cycle, deposition, and anthropogenic effects have primarily shown a net loss of C from ecosystems, potentially offsetting the effects of N deposition alone. In addition, LULC changes can modify the properties of N deposition by altering the ground surface features and reactive N emission (Lu et al., 2016). In addition, this model also lacks other essential elements like phosphorus.