Large-Scale Analysis of the Spatiotemporal Changes of Net Ecosystem Production in Hindu Kush Himalayan Region

: The Hindu Kush Himalayan (HKH) region is one of the most ecologically vulnerable regions in the world. Several studies have been conducted on the dynamic changes of grassland in the HKH region, but few have considered grassland net ecosystem productivity (NEP). In this study, we quantitatively analyzed the temporal and spatial changes of NEP magnitude and the inﬂuence of climate factors on the HKH region from 2001 to 2018. The NEP magnitude was obtained by calculating the difference between the net primary production (NPP) estimated by the Carnegie–Ames Stanford Approach (CASA) model and the heterotrophic respiration (Rh) estimated by the geostatistical model. The results showed that the grassland ecosystem in the HKH region exhibited weak net carbon uptake with NEP values of 42.03 gC · m − 2 · yr − 1 , and the total net carbon sequestration was 0.077 Pg C. The distribution of NEP gradually increased from west to east, and in the Qinghai–Tibet Plateau, it gradually increased from northwest to southeast. The grassland carbon sources and sinks differed at different altitudes. The grassland was a carbon sink at 3000–5000 m, while grasslands below 3000 m and above 5000 m were carbon sources. Grassland NEP exhibited the strongest correlation with precipitation, and it had a lagging effect on precipitation. The correlation between NEP and the precipitation of the previous year was stronger than that of the current year. NEP was negatively correlated with temperature but not with solar radiation. The study of the temporal and spatial dynamics of NEP in the HKH region can provide a theoretical basis to help herders balance grazing and forage.


Introduction
Grasslands occupy approximately 40% of the world's continental area and contribute approximately 30% to global carbon fixation [1][2][3]. Grasslands are indispensable members of the global terrestrial carbon cycle. Grasping the carbon budget of grassland ecosystems can promote understanding of the service functions provided by grasslands and potential response to the climate system [4], and this information can ultimately help achieve the goal of sustainable and rational use of grassland resources.
Net ecosystem production (NEP) is the most direct indicator of the capacity of carbon sources or sinks and is defined as the difference between net primary production (NPP) and heterotrophic respiration (Rh) [5]. Research on grassland ecosystem carbon exchange has made much progress, suggesting that grassland ecosystems sometimes act as potential carbon sequestrators or function in near equilibrium [6][7][8][9] and could balance the global carbon stock and budget of temperate grasslands [1] and alpine grasslands [10,11]. There is still great uncertainty regarding the carbon sources and sinks of grassland ecosystems because grassland may sometimes behave as carbon sources. Most of the fluctuations in these carbon sources and sinks are caused by dramatic changes in climate [2]. For example, in the event of drought or extreme wetness, grasslands may release carbon dioxide into the atmosphere [12,13]. Climate fluctuations usually limit the growth of grasslands and are the main cause of interannual changes in carbon sinks because they greatly reduce NPP and NEP [14]. However, most studies have focused on small regional scales of eddy covariance flux observations and fixed-point experiments. Although there have been studies using machine learning to upscale based on point observation data in order to achieve large-scale inversion [15][16][17], the upscaling process is a function of the flux data available and can induce uncertainty, particularly in regions that have few flux sites [18]. Therefore, the NEP inversion model based on remote sensing data is the guarantee for realizing large-area research, and it is more conducive to evaluating macroscopic patterns. With the development of remote sensing technology, many vegetation indexes have been found to be closely related to the carbon cycle, and large-scale climate data have been generated, which provide new methods for estimating NEP over large areas.
The Hindu Kush Himalayan (HKH) region is a massive alpine landform feature. Approximately half of the HKH region is composed of grassland, which serves as a resource for supporting 240 million mountain people who are dependent on nomadic livestock grazing [19,20]. In the context of global climate change, grassland ecosystems in the HKH region are relatively fragile and are also highly sensitive to various climate disturbances. In recent years, due to overgrazing and extreme weather events, grassland ecosystems in the HKH region have undergone extensive degradation, causing continuous fluctuations in carbon sources and sinks. Recent studies have shown that changes in grassland dynamics and productivity have occurred in watershed or component areas in the HKH region [20][21][22]. For example, Panday and Chimire studied changes in the normalized difference vegetation index (NDVI) from 1982 to 2016 in the HKH region and found that grassland showed an overall greening trend in NDVI magnitude [21]. Abbas et al. analyzed the response of grassland phenology and productivity to climate change in the Upper Indus Basin and explained that grassland phenology and productivity differed in different bioclimate regions [20]. Qamer et al. found that grassland productivity varied with altitude in the Hindu Kush Karakoram region and that the growth conditions of grassland improved from 2001 to 2012 [22]. However, there have been few studies on the carbon cycle in the HKH region. Only part of the research has focused on the Qinghai-Tibet Plateau and mainly focused on site observations of eddy covariance [23][24][25]. There is lack of large-scale studies on the carbon cycle throughout the HKH region; thus, the role of grassland ecosystems in the HKH region in the global carbon cycle remains unclear. Hence, a systematic assessment is needed to understand the grassland carbon dynamics in the HKH region as a basis for devising appropriate management policies.
In this context, we aimed to estimate the carbon cycle of grassland in the HKH region on a large scale based on remote sensing to help us understand the carbon source and sink situation in this area. Based on this, our study mainly focused on the following three aspects of research: (1) estimate the spatial distribution characteristics of grassland NEP based on remote sensing and climate data and analyze differences in spatial distribution between NEP and NPP; (2) study the interannual fluctuations of grassland NEP in the HKH region through statistical trend analysis; and (3) analyze the responses of grassland NEP and NPP to climate change.

Materials and Methods
In this study, we first collected remote sensing data, meteorological data, and land cover data as input data for the Carnegie-Ames Stanford Approach (CASA) model to estimate NPP. Then, we obtained the soil attribute data and combined it with the geostatistical Rh model to estimate Rh. Finally, NEP was obtained from the difference between NPP and Rh, and partial correlation analysis was used to explore the response of NPP and NEP to climate factors.

Study Area
The HKH region is the highest mountain and plateau region in the world (Figure 1). It includes the world's highest peak, Mount Everest, and has many peaks above 8000 m. The region includes the whole of Bhutan and Nepal and parts of mountainous and hilly areas in six other countries (Afghanistan, Bangladesh, China, India, Myanmar, and Pakistan) [26]. This region covers an area of more than 4.3 million km 2 , and it is the third-largest cryosphere in the world and the source of many large rivers, such as Brahmaputra, Indus, and Yellow Rivers [27]. Approximately 240 million people live in the high mountains and depend on nomadism, and 1.65 billion people live in the lower basin [28]. Precipitation in this area is mainly driven by the Indian monsoon and mid-latitude westerly winds in summer and winter, respectively [29]. In terms of precipitation gradient, the eastern and southern regions of the HKH region are relatively humid [30]. Grasslands are mainly distributed in the Qinghai-Tibet Plateau and the Hindu Kush Mountain. In the context of climate change, dynamic changes of grasslands in this region and the state of the carbon cycle will affect the local nomadic populations. Therefore, studying the carbon sources and sinks of grasslands in this region can provide important information for policymakers. In this study, we first collected remote sensing data, meteorological data, and land cover data as input data for the Carnegie-Ames Stanford Approach (CASA) model to estimate NPP. Then, we obtained the soil attribute data and combined it with the geostatistical Rh model to estimate Rh. Finally, NEP was obtained from the difference between NPP and Rh, and partial correlation analysis was used to explore the response of NPP and NEP to climate factors.

Study Area
The HKH region is the highest mountain and plateau region in the world (Figure 1). It includes the world's highest peak, Mount Everest, and has many peaks above 8000 m. The region includes the whole of Bhutan and Nepal and parts of mountainous and hilly areas in six other countries (Afghanistan, Bangladesh, China, India, Myanmar, and Pakistan) [26]. This region covers an area of more than 4.3 million km 2 , and it is the thirdlargest cryosphere in the world and the source of many large rivers, such as Brahmaputra, Indus, and Yellow Rivers [27]. Approximately 240 million people live in the high mountains and depend on nomadism, and 1.65 billion people live in the lower basin [28]. Precipitation in this area is mainly driven by the Indian monsoon and mid-latitude westerly winds in summer and winter, respectively [29]. In terms of precipitation gradient, the eastern and southern regions of the HKH region are relatively humid [30]. Grasslands are mainly distributed in the Qinghai-Tibet Plateau and the Hindu Kush Mountain. In the context of climate change, dynamic changes of grasslands in this region and the state of the carbon cycle will affect the local nomadic populations. Therefore, studying the carbon sources and sinks of grasslands in this region can provide important information for policymakers.

Normalized Difference Vegetation Index (NDVI)
The NDVI data used in the research were mainly derived from the MOD13A3 product of MODIS (Table 1), which was acquired from NASA EARTHDATA (https://earthdata.nasa.gov/) and used to estimate the fraction of photosynthetically active radiation (FPAR) absorbed by vegetation. The temporal resolution of NDVI data is monthly, and the spatial resolution is 1 km. We obtained the NDVI time series from 2001 to 2018 and then performed preprocessing, including splicing, cropping, projecting coordinates to the UTM WGS-1984 coordinate system, and resampling to 0.25° spatial resolution.

Dataset
Spatial resolution Temporal resolution Source NDVI 1 km monthly MOD13A3

Normalized Difference Vegetation Index (NDVI)
The NDVI data used in the research were mainly derived from the MOD13A3 product of MODIS (Table 1), which was acquired from NASA EARTHDATA (https: //earthdata.nasa.gov/ accessed on 1 March 2021) and used to estimate the fraction of photosynthetically active radiation (FPAR) absorbed by vegetation. The temporal resolution of NDVI data is monthly, and the spatial resolution is 1 km. We obtained the NDVI time series from 2001 to 2018 and then performed preprocessing, including splicing, cropping, projecting coordinates to the UTM WGS-1984 coordinate system, and resampling to 0.25 • spatial resolution.

Meteorological Elements
Meteorological elements, including precipitation, air temperature, and solar radiation (Table 1), were derived from fifth-generation European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data (ERA5) (https://cds.climate.copernicus.eu/ accessed on 1 March 2021), which are an alternative to ERA-Interim reanalysis. We obtained the monthly averaged reanalysis product type from the ERA5 dataset, including 2 m temperature, total precipitation, and surface solar radiation. These data were all stored in NC format, so we performed rasterization processing. We implemented the format conversion Remote Sens. 2021, 13, 1180 4 of 16 from NC to TIFF through Python, cropped the extents, reprojected the data to the UTM WGS-1984 coordinate system, and resampled the data to 0.25 • spatial resolution with the nearest-neighbor method.

Soil Attribute Data
When inverting soil respiration, the variable soil organic carbon (SOC) density is essential, and its inversion needs to be calculated through soil attribute data. Soil attributes ( Table 1) were derived from the Harmonized World Soil Database (HWSD) v.1.2 (http: //www.fao.org/ accessed on 1 March 2021), which was developed by the International Institute of Applied Systems Analysis and Food and Agriculture Organization of the United Nations. The HWSD data are mainly composed of hwsd.bil, HWSD.mdb, and HWSD_META.mdb. Four soil parameters (bulk density, organic carbon, gravel content, and soil depth) are needed to calculate the SOC density. The soil attribute data were rasterized, cropped, resampled, and reprojected to make them consistent with NDVI and meteorological data.

Land Cover Type Data
The land cover type data were obtained from MCD12C1, MODIS, and NASA; the spatial resolution of these data is 0.05 • , and the temporal resolution is yearly (Table 1). These data contain five land cover classification schemes, and the International Geosphere and Biosphere Program (IGBP) classification scheme was selected in our study. After the data were acquired, they needed to be cropped, reprojected, and resampled in Arc GIS 10.6 to obtain the input data required by the model with a 0.25 • spatial resolution and UTM WGS-1984 geographic coordinate system. Finally, we chose grasslands and savannas in the 2010 IGBP classification as the grassland type input in the CASA model.

Net Primary Production (NPP) Estimation Model
The CASA model is a remote sensing-based process model. It performs well in longterm series and large-scale estimation of NPP [31,32]. The two indispensable components in the CASA model are absorbed photosynthetically active radiation (APAR) (MJ·m −2 ) and light use efficiency (LUE) factor ε (gC·MJ −1 ) [33][34][35]. For a particular pixel x at month t, NPP is calculated by the following equation: where SOL is the surface solar radiation downwards (MJ·m −2 ). The constant 0.5 indicates that the solar energy used by vegetation photosynthesis is half of the total incident solar radiation [35]. The FPAR is a linear expression of NDVI and simple ratio (SR): Remote Sens. 2021, 13, 1180 5 of 16 where SR max and SR min correspond to the 95% and 5% quantiles of the normal distribution of all pixels from monthly SR, respectively.
where T(x,t) is the temperature stress coefficients, and W(x,t) is the water stress coefficient. ε max is the maximal LUE under ideal conditions with a value of 0.541 [36].

Geostatistical Model of Heterotrophic Respiration
The annual Rh is separated by annual total soil respiration using an empirical relationship through multiple analyses [37]. From multiple analysis, a global expression between Rh and soil respiration (R s ) was proposed: The geostatistical R s model was developed on the basis of the global empirical soil respiration model proposed by Raich [38]. Then, Yu et al. [39] found that soil organic carbon density and soil respiration are strongly correlated, which is an additional control factor for the spatial change of soil respiration. Therefore, the soil organic carbon density variable was added to the global experience model. The parameters in the model were improved through measured data, and the following model was proposed [39]: where R s is the mean soil respiration rates (gC·m −2 ·day −1 ); SOC density is the soil organic carbon density (kg·m −2 ); T is the mean monthly air temperature (°C); P is the mean monthly precipitation (cm); and α, β, P 0 , and K are all constant with values of 1.83, −0.006, 2.972, and 5.657, respectively [39]. Soil organic carbon density is directly related to soil properties [40,41], and its expression is as follows: where BD is the soil bulk density (g·m −3 ); TOC represents the topsoil organic carbon (%); Depth is the depth of soil layer (cm), which is 20 cm in our study; and Gravel denotes the gravel content (%).

Trend Analysis of Grassland Net Ecosystem Production (NEP)
The ordinary least-squares method was selected to calculate a linear regression of grassland NEP over HKH in order to reveal the temporal trend of NEP. The formula is expressed as follows: where slope is the interannual rate of NEP change; n is 18 for years from 2001 to 2018; i is 1 for the year 2001, 2 for the year 2002, and so on; and NEP i is the value of annual NEP in time of i year.

Partial Correlation Analysis
In order to evaluate the effects of a single climatic factor variable on NEP and to exclude interference from other climatic factors, we used a statistical partial correlation analysis method, which has been widely applied by many previous studies to explore Remote Sens. 2021, 13, 1180 6 of 16 the interaction between vegetation and a single climate element [42][43][44][45][46]. When two other climatic factors are used as dependent variables, the partial correlation coefficients of NEP and individual climatic factor can be calculated by the following equation: where R xy,z is the partial correlation coefficient between x and y excluding the impact of variable z; x and y are dependent variables, and z is the control variable. R xy , R xz , and R yz are the simple correlation coefficients between x, y, and z, respectively. We took R xy as an example to calculate the correlation coefficients between NEP and climatic elements: where x i and y i are the annual grassland NEP and climate factor, respectively; x and y are the average NEP and mean climate factors from 2001 to 2018, respectively. Finally, a bilateral t-test was implemented to assess the significance of the partial correlation coefficients with a significance level of 0.05.

Validation of the NPP Values
We compared the monthly NPP estimated by the CASA model with the monthly NPP observed at the flux sites in the FLUXNET2015 database (https://fluxnet.org/data/ fluxnet2015-dataset accessed on 1 March 2021). The flux site observed gross primary production (GPP), so this value needed to be multiplied by 0.55 (the ratio of NPP/GPP) to convert it into NPP [47]. For all sites (Table 2), the NPP measured at the site and estimated by the model through cross-validation had a regression coefficient of 0.74 and root mean square error (RMSE) of 25.65 gC·m −2 (Figure 2), which indicated a high precision of the CASA model estimation. There are two main ways to verify the accuracy of regional NEP model simulations: one is to compare the simulated and measured data, and the other is to compare the simulation among different models [48]. Because the measured data in the HKH region are limited, we compared the average NEP values of the grassland ecosystems with prior studies. We mainly compared NEP values from the Qinghai-Tibet Plateau because prior studies are mostly concentrated there (Table 3). NEP and net ecosystem carbon exchange (NEE) are numerically equal [49][50][51]. There were differences between NEP values in this study and prior studies, but the NEP values of our study were basically within the range of values given in the published literature. This difference might have been mainly caused by the uncertainty of the Rh model, which was established by Bond-Lamberty et al. [37] through statistical analysis of measured data around the world. The following problems occurred: (1) The database used to establish the model had no measured data throughout the HKH region, which led to deviations in Rh when the model was applied in the HKH region. The annual soil respiration (Rh) of the model was mainly derived from the expansion of the soil respiration (Rh) during the growing season, which would cause an overestimation of Rh. The soil microbial activity in the growing season and the decomposition ability were both strong, which would produce a large amount of carbon. In the nongrowing season, due to the decrease of temperature, the activity of microorganisms would weaken, and the carbon production would decrease. (3) A model based on a global scale might not be suitable for research on a regional scale. When the model was constructed, the strong spatial heterogeneity of the ground surface was not considered, so it might be biased when applied in the HKH region with strong heterogeneity.

Reliability Analysis of the NEP Values
There are two main ways to verify the accuracy of regional NEP model simulations: one is to compare the simulated and measured data, and the other is to compare the simulation among different models [48]. Because the measured data in the HKH region are limited, we compared the average NEP values of the grassland ecosystems with prior studies. We mainly compared NEP values from the Qinghai-Tibet Plateau because prior studies are mostly concentrated there (Table 3). NEP and net ecosystem carbon exchange (NEE) are numerically equal [49][50][51]. There were differences between NEP values in this study and prior studies, but the NEP values of our study were basically within the range of values given in the published literature. This difference might have been mainly caused by the uncertainty of the Rh model, which was established by Bond-Lamberty et al. [37] through statistical analysis of measured data around the world. The following problems occurred: (1) The database used to establish the model had no measured data throughout the HKH region, which led to deviations in Rh when the model was applied in the HKH region. (2) The annual soil respiration (Rh) of the model was mainly derived from the expansion of the soil respiration (Rh) during the growing season, which would cause an overestimation of Rh. The soil microbial activity in the growing season and the decomposition ability were both strong, which would produce a large amount of carbon. In the nongrowing season, due to the decrease of temperature, the activity of microorganisms would weaken, and the carbon production would decrease. (3) A model based on a global scale might not be suitable for research on a regional scale. When the model was constructed, the strong spatial heterogeneity of the ground surface was not considered, so it might be biased when applied in the HKH region with strong heterogeneity.   Using the basic input data sets and the optimized CASA model, the grassland NPP in the HKH region from 2001 to 2018 was estimated and its distribution characteristics were analyzed (Figure 3a). Figure 3a shows that the NPP showed an increasing trend along the direction of Hindu Kush Himalayas from west to east, and the gradient distribution of NPP gradually but very obviously increased from northwest to southeast in the Qinghai-Tibet Plateau of China. The average grassland NPP throughout the HKH region was  (Figure 3b), the NPP of grassland in 2018 was higher than that in 2001, and most regions had a growing trend. Among them, regions with a growth magnitude of 0-100 gC·m −2 accounted for the highest proportion (59.8%), and the proportion of regions with NPP growth magnitude greater than 100 gC·m −2 was 9.1%. The growth trend showed that the hydrothermal conditions in 2018 were better than those in 2001, the photosynthesis of grassland was strong, and productivity was high.

Spatial Characteristics of NPP
Using the basic input data sets and the optimized CASA model, the grassland NPP in the HKH region from 2001 to 2018 was estimated and its distribution characteristics were analyzed (Figure 3a). Figure 3a shows that the NPP showed an increasing trend along the direction of Hindu Kush Himalayas from west to east, and the gradient distribution of NPP gradually but very obviously increased from northwest to southeast in the Qinghai-Tibet Plateau of China. The average grassland NPP throughout the HKH region was 294 ± 216 gC·m −2 during 2001 to 2018. In the western part of the HKH region, namely Afghanistan's Hindu Kush Mountain and the western part of the Tibetan Plateau, NPP was less than 200 gC·m −2 , accounting for 46.3% of the total grassland area. Areas with NPP greater than 600 gC·m −2 were mainly distributed in the eastern part of the Qinghai-Tibet Plateau, accounting for 12.4%. Judging from the difference in NPP between 2018 and 2001 (Figure 3b), the NPP of grassland in 2018 was higher than that in 2001, and most regions had a growing trend. Among them, regions with a growth magnitude of 0-100 gC·m −2 accounted for the highest proportion (59.8%), and the proportion of regions with NPP growth magnitude greater than 100 gC·m −2 was 9.1%. The growth trend showed that the hydrothermal conditions in 2018 were better than those in 2001, the photosynthesis of grassland was strong, and productivity was high.

Spatial Characteristics of NEP
The temporal and spatial distribution of NEP in the HKH region was determined by estimating NPP and soil respiration. A positive value of NEP indicates the area is a carbon sink, and a negative value indicates a carbon source [59]. From the perspective of spatial distribution characteristics (Figure 4a), the interannual NEP accumulation in the HKH region from 2001 to 2018 showed a significant upward trend from west to east, which was consistent with the trend of NPP. The grassland ecosystems in the HKH region showed a weak carbon sink, with an average NEP of 42 ± 197 gC·m −2 , and the total net C sequestration was approximately 0.077 Pg C. In the past 18 years, the average amount of C sequestered by grassland ecosystems in the C sink area was 214.4 gC·m −2 , and the total amount of carbon sequestered was approximately 0.187 Pg C. Carbon storage areas were mainly distributed in northern India, Nepal, and the southeastern Tibetan Plateau, accounting for

Spatial Characteristics of NEP
The temporal and spatial distribution of NEP in the HKH region was determined by estimating NPP and soil respiration. A positive value of NEP indicates the area is a carbon sink, and a negative value indicates a carbon source [59]. From the perspective of spatial distribution characteristics (Figure 4a), the interannual NEP accumulation in the HKH region from 2001 to 2018 showed a significant upward trend from west to east, which was consistent with the trend of NPP. The grassland ecosystems in the HKH region showed a weak carbon sink, with an average NEP of 42 ± 197 gC·m −2 , and the total net C sequestration was approximately 0.077 Pg C. In the past 18 years, the average amount of C sequestered by grassland ecosystems in the C sink area was 214.4 gC·m −2 , and the total amount of carbon sequestered was approximately 0.187 Pg C. Carbon storage areas were mainly distributed in northern India, Nepal, and the southeastern Tibetan Plateau, accounting for approximately 52.2% of the total grassland area. At the same time, the average C emissions of grassland ecosystems in the carbon source sector was 116 gC·m −2 , and the total C emissions were approximately 0.

Interannual Variations
The mean annual NPP and NEP of grasslands in the HKH region showed a growth trend from 2001 to 2018 (

Time Series Change Trend Distribution
The Mann-Kendall significance test was chosen to analyze the slope change trend of NEP and NPP, and significance level α = 0.05 (|z|≥1.96) was selected to test the significance. There were obvious differences between NPP (Figure 6a) and NEP (Figure 6b), and the area of increase in NPP was more than that of NEP. The area that showed an increased trend in NPP (63.7%) was greater than that in NEP (53.4%). The area where NPP exhibited a significant increase accounted for approximately 13.1% of the total grassland area, and the area where NEP showed a significant increase accounted for about 8.6% of the total grassland area. The areas with significant increases of NPP and NEP were mainly concentrated in the south of the Afghanistan's Hindu Kush Mountain and the northeast of the Qinghai-Tibet Plateau. The region where NEP showed a downward trend accounted for approximately 46.7% of the total grassland area and was mainly mapped to the southwest of the Qinghai-Tibet Plateau.

Time Series Change Trend Distribution
The Mann-Kendall significance test was chosen to analyze the slope change trend of NEP and NPP, and significance level α = 0.05 (|z| ≥ 1.96) was selected to test the significance. There were obvious differences between NPP (Figure 6a) and NEP (Figure 6b), and the area of increase in NPP was more than that of NEP. The area that showed an increased trend in NPP (63.7%) was greater than that in NEP (53.4%). The area where NPP exhibited a significant increase accounted for approximately 13.1% of the total grassland area, and the area where NEP showed a significant increase accounted for about 8.6% of the total grassland area. The areas with significant increases of NPP and NEP were mainly concentrated in the south of the Afghanistan's Hindu Kush Mountain and the northeast of the Qinghai-Tibet Plateau. The region where NEP showed a downward trend accounted for approximately 46.7% of the total grassland area and was mainly mapped to the southwest of the Qinghai-Tibet Plateau.
cance. There were obvious differences between NPP (Figure 6a) and NEP (Figure 6b), and the area of increase in NPP was more than that of NEP. The area that showed an increased trend in NPP (63.7%) was greater than that in NEP (53.4%). The area where NPP exhibited a significant increase accounted for approximately 13.1% of the total grassland area, and the area where NEP showed a significant increase accounted for about 8.6% of the total grassland area. The areas with significant increases of NPP and NEP were mainly concentrated in the south of the Afghanistan's Hindu Kush Mountain and the northeast of the Qinghai-Tibet Plateau. The region where NEP showed a downward trend accounted for approximately 46.7% of the total grassland area and was mainly mapped to the southwest of the Qinghai-Tibet Plateau.

The Impact of Climate on NPP
To explore the main influencing factors that cause interannual changes of NPP in different regions, we calculated the partial correlation coefficients between NPP and precipitation, air temperature, and solar radiation (Figure 7a-c). Precipitation was the main

The Impact of Climate on NPP
To explore the main influencing factors that cause interannual changes of NPP in different regions, we calculated the partial correlation coefficients between NPP and precipitation, air temperature, and solar radiation (Figure 7a-c). Precipitation was the main limiting factor for NPP (Figure 7a). The area of steppe with NPP that was positively correlated with the precipitation was approximately 1.22 × 10 6 km 2 , which accounted for about 67.3% of the total grassland area. The areas with a significant positive correlation between NPP and precipitation were mainly distributed in the Hindu Kush Mountains and the western areas of the Qinghai-Tibet Plateau. The partial correlation between NPP and temperature ( Figure 7b) and solar radiation (Figure 7c) was insignificant over most of the area. The total areas of nonsignificant correlation were 1.68 × 10 6 and 1.7 × 10 6 km 2 , accounting for 92.1% and 93.5% of the total grassland area, respectively. The distribution of partial correlation showed that there was no strong correlation between NPP and temperature and solar radiation in most areas of the study area, which indicate that temperature and solar radiation are not the limiting climatic factors that affect the NPP of the HKH grassland ecosystem. limiting factor for NPP (Figure 7a). The area of steppe with NPP that was positively correlated with the precipitation was approximately 1.22 × 10 6 km 2 , which accounted for about 67.3% of the total grassland area. The areas with a significant positive correlation between NPP and precipitation were mainly distributed in the Hindu Kush Mountains and the western areas of the Qinghai-Tibet Plateau. The partial correlation between NPP and temperature ( Figure 7b) and solar radiation (Figure 7c) was insignificant over most of the area. The total areas of nonsignificant correlation were 1.68 × 10 6 and 1.7 × 10 6 km 2 , accounting for 92.1% and 93.5% of the total grassland area, respectively. The distribution of partial correlation showed that there was no strong correlation between NPP and temperature and solar radiation in most areas of the study area, which indicate that temperature and solar radiation are not the limiting climatic factors that affect the NPP of the HKH grassland ecosystem. In order to explore the response of NEP to different climatic factors, we analyzed the partial correlation coefficients between NEP and annual accumulated precipitation, annual average air temperature, and annual total solar radiation from 2001 to 2018 and im-

The Impact of Climate on NEP
In order to explore the response of NEP to different climatic factors, we analyzed the partial correlation coefficients between NEP and annual accumulated precipitation, annual average air temperature, and annual total solar radiation from 2001 to 2018 and implemented a significance test (Figure 7d-f). Same as NPP, precipitation was also the dominant factor in NEP changes (Figure 7d). The area of steppe with NEP that was positively correlated with precipitation was approximately 1.15 × 10 6 km 2 , which accounted for about 63.1% of the total HKH grassland area. In contrast to NPP, the partial correlation between NEP and temperature was mainly negative correlation with a partial correlation coefficient of −0.18 across the HKH grassland ecosystem (Figure 7e). The total area of negative correlation was 1.29 × 10 6 km 2 , accounting for approximately 70.8% of the total grassland area, and the significant negative correlation accounted for 17.7% of the total grassland area, which was mainly distributed in the Hindu Kush Mountains and the western region of the Tibetan Plateau. The correlation between NEP and solar radiation (Figure 7f) was also significantly different from NPP and was mainly in the Hindu Kush Mountains. NEP and solar radiation had a positive correlation, while NPP was negatively correlated with solar radiation.

Sources of Uncertainty
There are two main reasons for the differences in the tower-based flux and CASA estimates of NPP. One reason is the mismatch of the footprint; the flux measurement device is typically a tower-based eddy covariance system with a flux footprint of less than 1 km 2 , and they are subsamples of the CASA model pixel (0.25 • ) [60,61]. The other reason is that GPP measurements have systematic and random errors. Tower GPP is calculated as the difference between NEE and ecosystem respiration (Reco). Estimates of Reco at flux tower sites are typically created using nighttime fluxes of NEE [62]. However, the above-canopy fluxes are decoupled from the surface, leading to underestimates of nighttime respiration; thus, GPP exhibits deviations [63].

The Potential Effect of the NEP Spatial Distribution
The distribution of carbon sources and sinks throughout the HKH region were mainly affected by altitude and hydrothermal conditions, which is consistent with previous studies on the alpine grassland ecosystem of the Qinghai-Tibet Plateau [64,65]. As the altitude increased, the temperature decreased (Figure 8c), the solar radiation increased (Figure 8d), and the precipitation varied with altitude ( Figure 8b). The precipitation at 0-2000 m was the highest, followed by that at 3000-4000 m, and the precipitation above 5000 m was the lowest. With the increase in altitude, the areas of grassland carbon sources and sinks both showed an increasing trend (Figure 8a), reaching a peak in the elevation range of 4000-5000 m, with a total net carbon sequestration of 77.25 kg C, mainly due to the widespread distribution of grassland (46.8%) in this altitude range. Above an altitude of 5000 m, the average NEP of grassland was the lowest because the area has limited precipitation and low temperature, which restricts the growth of grassland and reduces the carbon sequestration capacity of grassland. In the height range of 3000-4000 m, grasslands had the strongest carbon fixation capacity, with an average carbon fixation capacity of 283.4 gC·m −2 . The average precipitation in this altitude range was relatively high at 880 mm, which provided sufficient water to the grasslands and resulted in strong photosynthesis, and these grasslands fixed a large amount of carbon from the atmosphere. In the height range of 0-2000 m, grassland had the strongest ability to release carbon into the atmosphere, with an average carbon emission of 151.3 gC·m −2 . The high precipitation in this area reduces solar radiation, reduces photosynthesis, enhances respiration, and causes grassland to release carbon into the atmosphere. Among all the height classifications, the areas of the carbon sinks were larger than the areas of the carbon sources in the 3000-4000 m and 4000-5000 m height ranges, which were manifested as carbon sinks. In contrast, grasslands below 3000 m and above 5000 m were manifested as carbon sources.
283.4 gC·m . The average precipitation in this altitude range was relatively high at 880 mm, which provided sufficient water to the grasslands and resulted in strong photosynthesis, and these grasslands fixed a large amount of carbon from the atmosphere. In the height range of 0-2000 m, grassland had the strongest ability to release carbon into the atmosphere, with an average carbon emission of 151.3 gC·m −2 . The high precipitation in this area reduces solar radiation, reduces photosynthesis, enhances respiration, and causes grassland to release carbon into the atmosphere. Among all the height classifications, the areas of the carbon sinks were larger than the areas of the carbon sources in the 3000-4000 m and 4000-5000 m height ranges, which were manifested as carbon sinks. In contrast, grasslands below 3000 m and above 5000 m were manifested as carbon sources.

The Potential Effect of NEP Temporal Change
Precipitation is the main factor affecting NEP changes in alpine grassland ecosystems [62,64]. By analyzing the relationship between precipitation and NEP, it was found that changes in the two were not synchronized, and there was a lag effect. The peaks and valleys of NEP were always one year behind the time when the peaks and valleys of precipitation appeared. For example, the peak of precipitation appeared in 2010, while the peak

The Potential Effect of NEP Temporal Change
Precipitation is the main factor affecting NEP changes in alpine grassland ecosystems [62,64]. By analyzing the relationship between precipitation and NEP, it was found that changes in the two were not synchronized, and there was a lag effect. The peaks and valleys of NEP were always one year behind the time when the peaks and valleys of precipitation appeared. For example, the peak of precipitation appeared in 2010, while the peak of NEP appeared in 2011; the valley of precipitation appeared in 2006, and the valley of NEP appeared in 2007. By analyzing the partial correlation between NEP and the precipitation of the current year (Figure 7d-f) and the precipitation of the previous year (Figure 9a-c), it was found that the correlation between NEP and the precipitation of the previous year was strong and positive. The positive correlation covered 74% of the total area, which was 10.8% more than the positive correlation between NEP and the precipitation of the current year. This result also showed that if the precipitation in the previous year was sufficient, grassland growth would increase in the current year, and the grassland would have a strong carbon sequestration ability. The correlation between NEP and the temperature of the previous year was also stronger than that of the current year, with 44.1% of the area showing positive correlation, which was 14.9% more than the correlation between NEP and the temperature of the current year. Therefore, good hydrothermal conditions in the previous year will promote growth of NEP in the current year. The correlation between NEP and the solar radiation of the previous year showed an insignificant correlation (92.3%) in most areas, which was basically consistent with the correlation between NEP and the solar radiation of the current year (94.3%).  Figure  9a-c), it was found that the correlation between NEP and the precipitation of the previous year was strong and positive. The positive correlation covered 74% of the total area, which was 10.8% more than the positive correlation between NEP and the precipitation of the current year. This result also showed that if the precipitation in the previous year was sufficient, grassland growth would increase in the current year, and the grassland would have a strong carbon sequestration ability. The correlation between NEP and the temperature of the previous year was also stronger than that of the current year, with 44.1% of the area showing positive correlation, which was 14.9% more than the correlation between NEP and the temperature of the current year. Therefore, good hydrothermal conditions in the previous year will promote growth of NEP in the current year. The correlation between NEP and the solar radiation of the previous year showed an insignificant correlation (92.3%) in most areas, which was basically consistent with the correlation between NEP and the solar radiation of the current year (94.3%).

Conclusions
In this study, NEP of the HKH region from 2001 to 2018 was estimated based on the CASA model and geostatistical soil respiration model. We found that our results were within the range of the results of previous studies [52][53][54][55][56][57][58]. Although the uncertainty of the input data and the geostatistical soil respiration model resulted in deviations to the NEP estimation, our results reflected the temporal and spatial changes of the carbon

Conclusions
In this study, NEP of the HKH region from 2001 to 2018 was estimated based on the CASA model and geostatistical soil respiration model. We found that our results were within the range of the results of previous studies [52][53][54][55][56][57][58]. Although the uncertainty of the input data and the geostatistical soil respiration model resulted in deviations to the NEP estimation, our results reflected the temporal and spatial changes of the carbon sources/sinks in the HKH region. The grassland ecosystems had the potential to sequester carbon and showed an increasing trend from west to east. The grassland ecosystems were carbon sinks at 3000-5000 m, while carbon sources occurred at both lower than 3000 m and higher than 5000 m. Grassland had the strongest net carbon sequestration ability at 3000-4000 m, while grassland above 5000 m had the weakest net carbon sequestration capacity. Moreover, grassland NEP exhibited a stronger positive correlation with precipitation of the previous year than with precipitation of the current year. NEP was negatively correlated with temperature and not significantly correlated with solar radiation, and there was no obvious lag relationship. The above results based on large-scale analysis can help us understand the temporal and spatial changes of grassland ecosystem NEP in the HKH region and its response to climate change, which is of great significance to the development of grassland policies in the HKH region.
Author Contributions: D.G. designed this research experiment, performed analysis and verification, and drafted the manuscript. X.S. administrated the study, gave some constructive suggestions, and applied for funding for this study. R.H. participated in the design of the study and revised the draft to improve the presentation of this study. X.Z. offered some suggestions on visualization and the draft. S.C. participated in the realization and improvement of algorithm programming. Y.J. polished the pictures in the article. Y.Z. participated in editing the paper and checking for grammatical errors. X.C. provided guidance on thinking and logic. All authors have read and agreed to the published version of the manuscript.