Analyzing NPP Response of Different Rangeland Types to Climatic Parameters over Mongolia

: Global warming threatens ecosystem functions, biodiversity, and rangeland productivity in Mongolia. The study analyzes the spatial and temporal distributions of the Net Primary Production (NPP) and its response to climatic parameters. The study also highlights how various land cover types respond to climatic ﬂuctuations from 2003 to 2018. The Boreal Ecosystem Productivity Simulator (BEPS) model was used to simulate the rangeland NPP of the last 16 years. Satellite remote sensing data products were mainly used as input for the model, where ground-based and MODIS NPP were used to validate the model result. The results indicated that the BEPS model was moderately effective (R 2 = 0.59, the Root Mean Square Error (RMSE) = 13.22 g C m − 2 ) to estimate NPP for Mongolian rangelands (e.g., grassland and sparse vegetation). The validation results also showed good agreement between the BEPS and MODIS estimates for all vegetation types, including forest, shrubland, and wetland (R 2 = 0.65). The annual total NPP of Mongolia showed a slight increment with an annual increase of 0.0007 Pg (0.68 g C per meter square) from 2003 to 2018 ( p = 0.82) due to the changes in climatic parameters and land cover change. Likewise, high increments per unit area found in forest NPP, while decreased NPP trend was observed in the shrubland. In conclusion, among the three climatic parameters, temperature was the factor with the largest inﬂuence on NPP variations ( r = 0.917) followed precipitation ( r = 0.825), and net radiation ( r = 0.787). Forest and wetland NPP had a low response to precipitation, while inter-annual NPP variation shows grassland, shrubland, and sparse vegetation were highly sensitive rangeland types to climate ﬂuctuations.


Introduction
Rangeland ecosystems are vital for farming and the environments that regulate the global carbon cycle [1]. Rangeland includes several vegetation types such as grasslands, shrublands, savannas, and marshes, which together cover about 40% of the world's land area [2,3]. Mongolia has 2.6% of the world's rangelands, constituting around 83% of the country's territory [4,5]. Although the livestock industry is the primary engine of the country's economy, Mongolian rangelands face unprecedented threats from climate change and overgrazing. Both of them have almost an equal impact on Mongolian rangelands [6,7], Bao et al. [47] have quantified spatiotemporal patterns of terrestrial NPP over Mongolia from 1982 to 2011, mainly focusing climatic constraints on NPP during the growing season. Bat-oyun et al. [48] focused on climatic and grazing effects on grassland productivity for 2005 to 2007; Lin et al. [49] calculated monthly NPP from 2000 to 2004. Although their study covered different periods, they focused only on growing season at a monthly time scale and relied on station data for meteorological input. However, it is still unclear how different rangeland types in Mongolia respond to climatic parameter changes. Considering this background, this study focuses on daily and annual spatial and temporal variation of rangeland NPP in Mongolia as a response to the fluctuation in three climatic parameters (precipitation, temperature, and solar radiation).
This study first examines the spatiotemporal variation of NPP over five different vegetation types due to the climatic parameters over the last 16 years (2003-2018) using the BEPS model. The second objective was to highlight the most vulnerable rangeland types to the climatic parameters, which is very urgent to intervene the adverse impacts and ensure environmental sustainability.

Study Area
The study has focused on Mongolia, which is located in Central Asia. In terms of climate, the country is located in the arid and semi-arid regions between the Great Siberian boreal forest and the Central Asian deserts. As a result, Mongolian rangelands are highly sensitive to environmental changes [60]. Mongolian rangelands have diverse ecosystems ranging from the Gobi Desert in the south to the forest-steppe in the north (Figure 1), with high inter-annual precipitation variations and fluctuations in annual and seasonal rangeland productivity.
how different rangeland types in Mongolia respond to climatic parameter changes. Considering this background, this study focuses on daily and annual spatial and temporal variation of rangeland NPP in Mongolia as a response to the fluctuation in three climatic parameters (precipitation, temperature, and solar radiation).
This study first examines the spatiotemporal variation of NPP over five different vegetation types due to the climatic parameters over the last 16 years (2003-2018) using the BEPS model. The second objective was to highlight the most vulnerable rangeland types to the climatic parameters, which is very urgent to intervene the adverse impacts and ensure environmental sustainability.

Study Area
The study has focused on Mongolia, which is located in Central Asia. In terms of climate, the country is located in the arid and semi-arid regions between the Great Siberian boreal forest and the Central Asian deserts. As a result, Mongolian rangelands are highly sensitive to environmental changes [60]. Mongolian rangelands have diverse ecosystems ranging from the Gobi Desert in the south to the forest-steppe in the north (Figure 1), with high inter-annual precipitation variations and fluctuations in annual and seasonal rangeland productivity.
The country has a total area of 1.565 million km 2 ; the border extends between the latitudes of 41°35′ N to 52°09′ N and the longitudes of 87°44′ E to 119°56′ E ( Figure 1). This country's elevation decreases gradually from west to east, with an average elevation of 1580m above sea level. The region has a wide variety of climatic conditions with four distinct seasons; it has high annual and diurnal temperature fluctuations and low rainfall. Annual total precipitation varies from 50 mm in the Gobi Desert region to about 400 mm in the northern mountain region, whereas the annual temperature ranges in the study area from −8 °C to 6 °C [61]. The majority of the annual rainfall (about 85%) falls in the growing season (May to September), where 50-60% of rain is concentrated in June, July, and August. Meanwhile, less than 20% of precipitation falls in the form of snow in winter [62]. There are 230-260 sunny days in a year [61], and the number of rainy days increases from south to north direction. The average total sunshine duration in Mongolia is about 2600-3300 h per year [63]. The country has a total area of 1.565 million km 2 ; the border extends between the latitudes of 41 • 35 N to 52 • 09 N and the longitudes of 87 • 44 E to 119 • 56 E (Figure 1). This country's elevation decreases gradually from west to east, with an average elevation of 1580m above sea level. The region has a wide variety of climatic conditions with four distinct seasons; it has high annual and diurnal temperature fluctuations and low rainfall. Annual total precipitation varies from 50 mm in the Gobi Desert region to about 400 mm in the northern mountain region, whereas the annual temperature ranges in the study area from −8 • C to 6 • C [61]. The majority of the annual rainfall (about 85%) falls in the growing season (May to September), where 50-60% of rain is concentrated in June, July, and August.
Meanwhile, less than 20% of precipitation falls in the form of snow in winter [62]. There are 230-260 sunny days in a year [61], and the number of rainy days increases from south to north direction. The average total sunshine duration in Mongolia is about 2600-3300 h per year [63].

NPP Estimation from BEPS Model
The study applied BEPS model simulation to examine the spatial-temporal of NPP variation over Mongolian rangelands. The model was operated in ArcGIS environments and continuously run for years from 2003 to 2018. The annual NPP at the spatial resolution of 500m was calculated as the sum of the daily NPP value between 1 January and 31 December, respectively (day of the year). At last, the statistical analysis of NPP at different land cover types was performed, and the spatial and temporal trend analysis of NPP was investigated. The general flowchart of the overall methodological process and the primary inputs is presented in Figure 2. The BEPS calculates Gross Primary Productivity (GPP), NPP, and Evapotranspiration (ET) for each pixel based on; (1) remotely sensed vegetation parameters such as Leaf Area Index (LAI), vegetation cover types, field capacity, tree cover, and wilting point, (2) topographic parameters such as elevation and latitude, and (3) daily meteorological data. Likewise, land cover data used to determine vegetation structure and type. Spatially distributed LAI is needed to calculate canopy conductance, photosynthesis, and leaf respiration. The model requires daily meteorological data to specify the environmental conditions of plant growth [64]. The DEM and soil data are for defining hydrological parameters for the hydrological module. Finally, the model needs annual tree cover, latitude, and plant type data (C 3 ).

NPP Estimation from BEPS Model
The study applied BEPS model simulation to examine the spatial-temporal of NPP variation over Mongolian rangelands. The model was operated in ArcGIS environments and continuously run for years from 2003 to 2018. The annual NPP at the spatial resolution of 500m was calculated as the sum of the daily NPP value between 1 January and 31 December, respectively (day of the year). At last, the statistical analysis of NPP at different land cover types was performed, and the spatial and temporal trend analysis of NPP was investigated. The general flowchart of the overall methodological process and the primary inputs is presented in Figure 2. The BEPS calculates Gross Primary Productivity (GPP), NPP, and Evapotranspiration (ET) for each pixel based on; (1) remotely sensed vegetation parameters such as Leaf Area Index (LAI), vegetation cover types, field capacity, tree cover, and wilting point, (2) topographic parameters such as elevation and latitude, and (3) daily meteorological data. Likewise, land cover data used to determine vegetation structure and type. Spatially distributed LAI is needed to calculate canopy conductance, photosynthesis, and leaf respiration. The model requires daily meteorological data to specify the environmental conditions of plant growth [64]. The DEM and soil data are for defining hydrological parameters for the hydrological module. Finally, the model needs annual tree cover, latitude, and plant type data (C3). The model consists of a soil water balance model, stomatal and mesophyll conductance model, photosynthesis, and respiration model, which eventually estimate daily and yearly accumulated NPP. In the BEPS, the NPP is computed as the difference between the The model consists of a soil water balance model, stomatal and mesophyll conductance model, photosynthesis, and respiration model, which eventually estimate daily and yearly accumulated NPP. In the BEPS, the NPP is computed as the difference between the total Agronomy 2021, 11, 647 5 of 19 carbon uptake through photosynthesis (GPP) and the carbon release due to respiration (R a ) [55,65]. The following Equations (1)-(3) were used to run the NPP [66].
where NPP is net primary productivity, GPP (gross primary productivity) is the carbon flux during photosynthesis. The subscripts 'sun' and 'shade' denote the sunlight and shaded components of daily mean photosynthesis rate (A) and LAI. The factor converts the GPP unit into g C m −2 day −1 . Daily R a is the autotrophic respiration, while R m and R g are the maintenance respiration and growth respiration by all living parts including fine roots and leaves, respectively [67-69].

Validation Analysis
The NPP estimated from the BEPS model was compared with the corresponding NPP measured from field observations at the exact time and geographical location. We conducted two statistical analyses to evaluate BEPS model accuracy. Firstly, the R 2 (Coefficient of determination) represents the BEPS model's fitness to simulate NPP compared to the ground observations (Equation (4)). The second statistical analysis to evaluate the model efficiency is the Root Mean Square Error (RMSE) from Equation (5).
where x and y represent NPP values estimated from the observed data and BEPS model, respectively, in the ith year; x represents the mean annual amount of NPP estimated from the ground, y represents the annual average NPP estimated from the model; n = 206 represents the number of biomass samples; and R x,y is the correlation coefficient between x and y. The R 2 value is the coefficient of determination which ranges from 0 to 1. The high R 2 value reflects how the model fitness in predicting NPP over the study area.
where x i is the ground-measured NPP, y i is the simulated NPP, and n is the number of samples. The smaller the RMSE value, the higher the accuracy of the model in estimating NPP.

Trend Analysis
The non-parametric Mann-Kendall (MK) test was used to assess trends in the time series of the variables [69]. The significance of the annual NPP and its driving three climatic parameters (temperature, precipitation, and net radiation) were tested at the 95% level, and MK statistics for time series was calculated as Equation (6): where Z i and Z j are mean NPP and climatic parameters value in the year i and j, respectively, i > j, n = 16 years, representing the length of time series from 2003 to 2018 and sgn (Z i − Z j ) is sign function defined by Equation (7): Agronomy 2021, 11, 647 6 of 19 The null hypothesis (H o ) is that there is no trend in the series, whereas the alternative hypothesis (H 1 ) is that an increasing or decreasing monotonic trend exists in the series. The presence of a statistically insignificant trend was evaluated based on the p-value. H o (no trend) is rejected if the p-value is less than a predefined significance level of 0.05. In addition to the linear regression method, the Sen's slope method was also used to estimate the slope of NPP and its driving parameters. If the time series data presents a linear trend, the true slope (change per unit time) of a trend can be estimated by the non-parametric index developed by Sen [70], which is based on the assumption of a linear trend [71]: where x i and x j are changing NPP and climate parameters value at the time i and j, respectively. The annual NPP and climate parameters trend is computed as an average change: a negative slope indicates a decreased trend and a positive slope denotes an increasing trend.

Statistical and Correlation Analysis
We extracted each climate parameter and NPP of five different vegetation types from 2003 to 2018 within the study area using the zonal statistical analysis. Then, Pearson correlation analysis was used to study the correlation between NPP and climatic parameters (temperature, precipitation, and net radiation). Their correlation level in different vegetation types was analyzed separately (9).
where x and y represent NPP values and climatic parameters (temperature, precipitation, and net radiation) in the ith year; x represents the mean annual NPP, yrepresents the mean annual value of climatic parameters (temperature, precipitation, and net radiation) from 2003 to 2018; n = 16 years, representing the length of time series. R x,y < 0 means a negative correlation coefficient between x and y and vice versa [72].

Data Acquisition and Preparation
All the input datasets were processed beforehand NPP estimation in the same coordinate system (latitude/longitude) and unified spatial resolution. The following section explains the input data characteristics and their preparation.

LAI Data
The Leaf Area Index (LAI) is one of the main input variables to most distributed ecosystem process models. In this study, the LAI product was derived from Moderate Resolution Imaging Spectroradiometer (MODIS) reflectance data. The MODIS LAI product (MCD15A3H) can be obtained from the Earth data (https://doi.org/10.5067/MODIS/ MCD15A3H.006, accessed on 1 January 2019), with 4-day temporal resolution by composing the pixel value of each four days from NASA's Terra and Aqua satellites in 500-m spatial resolution. In this study, we downloaded the 4-days temporal products of MCD15A3H for the years 2003-2018. Then, LAI gaps were interpolated to daily intervals from the original four days raster using the direct linear interpolation method.

Climate Data
The daily climate data (precipitation, air temperature, net surface radiation, wind speed, snow depth, and dew point temperature) for 2003-2018 were derived from the European Centre for Medium-Range Weather Forecast (ECMWF) data products. The ECMWF provides global gridded three-hour time-varying climate products (http://www. ecmwf.int/research/era, accessed on 10 December 2018) with a spatial resolution of 0.125 • × 0.125 • . The daily input data was estimated from 3 hourly gridded climate data, while the daily precipitation was derived as the sum of the daily rainfall. The net surface radiation was calculated as the sum of the surface net thermal radiation and surface net solar radiation. ECMWF meteorological data used in this study was successfully applied in modeling [73], including the BEPS model [74]. Besides that, data reliability was achieved through the correlation analysis between ECMWF climate data and ground data (70 stations), which showed R 2 values 0.64 and 0.70 for precipitation and temperature, respectively.

Land Cover and Tree Cover Data
Land cover plays a crucial role in terrestrial carbon biogeochemical processes, and thus, the land cover data is a significant input variable to the BEPS model. The global land cover dataset was downloaded from the European Space Agency (ESA) Climate Change Initiative Land Cover (CCI-LC) product through the project website (http://maps.elie.ucl. ac.be/CCI/viewer/download.php, accessed on 17 June 2019). The maps are available on an annual basis from 1992 to 2015 at 300 m spatial resolution. The CCI-LC classifies land cover into 22 types based on the Land Cover Classification System (LCCS) method. We transformed the original 22 classes of CCI-LC products into two separate codes to suit the model and vegetation analysis ( Table 1). The first was to estimate NPP through the BEPS model, while the second was for analyzing NPP variation for different vegetation types. Also, we clipped the global maps to derive the LC dataset for Mongolia between 2003 and 2015. Then, the 22 CCI-LC classes were aggregated into 18 categories as input to the BEPS model. We transformed the CCI-LC classes into nine general classes to estimate the NPP of each vegetation type for vegetation analysis. Besides the LC maps, the BEPS model requires an additional layer of percent tree cover. The MODIS has a percent tree cover product with a spatial resolution of 250 m from vegetation continuous vegetation fields (MOD44B) product at the yearly interval [75]. We derived the percent tree cover layers over Mongolia to cover the whole study period.  The soil-related variable required by BEPS is the Available soil Water-holding Capacity (AWC). AWC is defined as the water content difference between field capacity and permanent wilting point in the soil. Both, wilting point and field capacity at a spatial resolution of 0.083 • × 0.083 • were downloaded from the International Geosphere-Biosphere Program, Data, and Information System (IGBP-DIS; https://daac.ornl.gov/cgi-bin/dsviewer.pl? ds_id=569, accessed on 15 January 2019). Besides that, DEM and latitude data are also important input variables to the BEPS model. Hence, the 30 m DEM was downloaded from the Shuttle Radar Topography Mission (SRTM) (https://dwtkns.com/srtm30m/, accessed on 15 June 2019). Then, the latitude data was derived from the DEM data based on its y-coordinates.

Field Sampling Data and Data Collection
For model validation, we used the aboveground biomass (AGB) data derived from 47 stations across the country (Figure 1), which was provided by the Information and Research Institute of Meteorology, Hydrology, and Environment (IRIMHE) of Mongolia. The dominant plant community in the Mongolian rangeland is characterized by C 3 , coolseason species [77]. The AGB of grass was measured once per month: on the twenty-fifth of each months during the growing season (from May to September of the year) for 2003-2018. Forty-seven plots (1 × 1 m) were hand-harvested (clipped) at 1 cm above the ground surface. Later, it was dried in the lab until the weight remained constant and weighted using a balance with an accuracy of 0.1 g. The dried AGB values were expressed as gram per meter square area (g/m 2 ), then converted to NPP using a carbon conversion factor of 0.475 (for grass and foliage components) as C per gram dry mass [39]. Finally, the NPP model results were compared with the corresponding observed NPP from the same time and location.

The Accuracy of Estimating NPP from the BEPS Model
The correlation analysis examined between BEPS-simulated NPP with AGB gathered from May to September for the whole study period showed a moderate positive correlation (R 2 = 0.59) with root-mean-square-error (RMSE) of 13.22 g C m −2 (Figure 3a). The mean annual NPP over Mongolian rangelands was 0.1839 Pg C y −1 (1 Pg C = 10 15 g C) from 2003 to 2018. However, the derived NPP of this study is smaller than previous studies by Bao et al. [47] and Li et al. [49] (NPP correspond 0.3 Pg C y −1 and 0.71 Pg C y −1 , respectively. The results show discrepancies, possibly because different data input and methods were used to estimate vegetation production. For example, Bao et al. [47] have used NDVI products of 8 km spatial resolution as vegetation index, while we have used LAI of 500-m spatial resolution. In addition to that, previous studies have estimated monthly NPP using station-based monthly meteorological data. The BEPS model also considers the foliage clumping effect, which reduces sunlit leaves and increases shaded leaves. Therefore, the BEPS NPP value might be less than another model result. shrubland, forest and wetland. The comparison was made over randomly selected 461 locations covering all vegetation types. A good agreement (R 2 = 0.65) was found between BEPS and MODIS NPP over the study area (Figure 3b). Overall, the analyses indicated that the good efficiency of the BEPS model to estimate the NPP over Mongolian rangelands. However, the BEPS estimated NPP result was higher compared to the observed NPP, while MODIS NPP predicted much higher NPP than BEPS did in the study area ( Figure 3).

Inter-Annual Variations in NPP
The mean annual NPP's spatial distribution over Mongolia from 2003 to 2018 was estimated ( Figure 4). NPP spatially increases from the south to the north. While the minimum NPP was estimated over sparse vegetation in the south, the maximum NPP was concentrated over forested areas in the north. The spatial and temporal NPP trends based on the Mann-Kendall test and the Sen's slope have presented in Figures 4b and 5a. During the study period, the annual NPP across Mongolia has increased in 59% of total vegetated areas (up to 13.74 g C m −2 ) where the Mann-Kendall test indicates increment was not significant (Figure 4b). A decreasing trend was found in about 41% of the vegetated areas (up to -14.75 g C m −2 ), and 7.66% of them were significantly decreased for the spatial dimension. The increasing NPP trend in the region indicates that the vegetation has played an increasing role in carbon sequestration. This is consistent with the study of Li et al. [79], who reported that the NPP of middle and high latitude regions generally increased from 1961 to 2010. Li et al. [79] also noted that climate change was the second most crucial factor that controls global terrestrial NPP's inter-annual variability after atmospheric carbon dioxide. By contrast, Nemani et al. [80] suggested that climatic parameters mainly determine the global terrestrial NPP: temperature, precipitation, and net radiation; are the main affecting parameters of the inter-annual variability of NPP. Our research considered the 16- Although collected field sample plots for model validation were evenly distributed across the study area, most of them have located over the grazed rangelands (e.g., grassland and sparse vegetation). Hence we used MODIS/Terra Annual NPP (MOD17A3HGF, 500 m) product [78] to validate the BEPS model result for all vegetation types, including shrubland, forest and wetland. The comparison was made over randomly selected 461 locations covering all vegetation types. A good agreement (R 2 = 0.65) was found between BEPS and MODIS NPP over the study area (Figure 3b). Overall, the analyses indicated that the good efficiency of the BEPS model to estimate the NPP over Mongolian rangelands. However, the BEPS estimated NPP result was higher compared to the observed NPP, while MODIS NPP predicted much higher NPP than BEPS did in the study area (Figure 3).

Inter-Annual Variations in NPP
The mean annual NPP's spatial distribution over Mongolia from 2003 to 2018 was estimated ( Figure 4). NPP spatially increases from the south to the north. While the minimum NPP was estimated over sparse vegetation in the south, the maximum NPP was concentrated over forested areas in the north. The spatial and temporal NPP trends based on the Mann-Kendall test and the Sen's slope have presented in Figures 4b and 5a. During the study period, the annual NPP across Mongolia has increased in 59% of total vegetated areas (up to 13.74 g C m −2 ) where the Mann-Kendall test indicates increment was not significant (Figure 4b). A decreasing trend was found in about 41% of the vegetated areas (up to -14.75 g C m −2 ), and 7.66% of them were significantly decreased for the spatial dimension. The increasing NPP trend in the region indicates that the vegetation has played an increasing role in carbon sequestration. This is consistent with the study of Li et al. [79], who reported that the NPP of middle and high latitude regions generally increased from 1961 to 2010. Li et al. [79] also noted that climate change was the second most crucial factor that controls global terrestrial NPP's inter-annual variability after atmospheric carbon dioxide. By contrast, Nemani et al. [80] suggested that climatic parameters mainly determine the global terrestrial NPP: temperature, precipitation, and net radiation; are the main affecting parameters of the inter-annual variability of NPP. Our research considered the 16-y datasets of climatic parameters: temperature, precipitation, and net radiation were used to explore the patterns of inter-annual variability of NPP (Figure 5b-d). As shown in Figure 5a, the annual NPP has slightly increased from 0.169 Pg C y −1 in 2003 to 0.20 Pg C y −1 in 2018, with an annual increase of 0.0007 Pg C (p-value = 0.82). This can be attributed to the annual total precipitation and mean temperature slightly increased over the study period (with an average of 0.081 mm per year, 0.066 • C per year, respectively), while net radiation decreased by −0.018 W/m −2 per year from 2003 to 2018 (Figure 5b-d). This finding is also consistent with that of a previous study, which showed that overall trends of average annual temperature and precipitation have been increasing worldwide from 2000 to 2014 [81]. In addition, the ESA land cover map showed barren land has decreased over the study period, also supported by [82], and the vegetated area has increased. This land-use change could be a reason that annual NPP has risen slightly in the study area. Moreover, the annual NPP variation is probably due to changes in the vegetation growing season length in Mongolia [83].
an average of 0.081 mm per year, 0.066 °C per year, respectively), while net radiation decreased by -0.018 W/m −2 per year from 2003 to 2018 (Figure 5b-d). This finding is also consistent with that of a previous study, which showed that overall trends of average annual temperature and precipitation have been increasing worldwide from 2000 to 2014 [81]. In addition, the ESA land cover map showed barren land has decreased over the study period, also supported by [82], and the vegetated area has increased. This land-use change could be a reason that annual NPP has risen slightly in the study area. Moreover, the annual NPP variation is probably due to changes in the vegetation growing season length in Mongolia [83].
Specifically, lower NPP was found in 2007, 2009, 2010, and 2017 which were attributed to the combined effect of temperature, solar radiation, and precipitation. For example, high temperature and low rainfall were observed in 2007, 2009, 2010, and 2017, which inevitably brought about low NPP values for Mongolian rangelands (Figure 5b-d). Some studies have mentioned that droughts' influences certainly led to a decrement in NPP [84][85][86][87]. Recent research has also confirmed that Mongolia has experienced frequent droughts [24] and warming over the last two decades, which has reduced aboveground biomass by 20% to 65% [87]. At the vegetation type level, the NPP exhibited large variations across five major vegetation types ( Table 2). The forest has the highest NPP with an average value of 195.23 g C m −2 y −1 than other vegetation types. The accumulated annual NPP of the forest was less than grassland (0.05268 Pg C) because the forest only covers 10.83% of the vegetated lands. A large area of the country covered by grassland (41.74%) has the highest annual NPP (0.07693 Pg C). While the dominant rangeland type is the sparse vegetation (45.76%), its annual NPP was only 0.02717 Pg C (Table 2).   To further differentiate the NPP trends of different vegetation types, we performed a statistical analysis to derive annual NPP values per unit area. Annual NPP per square meter for different vegetation types differed greatly, and forest showed the highest annual NPP with an average value of 635.97 g C m −2 y −1 during the last 16 years ( Figure 6). Wetland vegetation had a slightly lower NPP than the forest, with an average value of 525.07 g C m −2 y −1 . The grassland NPP was about 241 g C m −2 y −1 , while shrubland and sparse vegetation were almost lower than 100 g C m −2 y −1 from 2003 to 2018.   (Figure 5b-d). Some studies have mentioned that droughts' influences certainly led to a decrement in NPP [84][85][86][87]. Recent research has also confirmed that Mongolia has experienced frequent droughts [24] and warming over the last two decades, which has reduced aboveground biomass by 20% to 65% [87].
At the vegetation type level, the NPP exhibited large variations across five major vegetation types ( Table 2). The forest has the highest NPP with an average value of 195.23 g C m −2 y −1 than other vegetation types. The accumulated annual NPP of the forest was less than grassland (0.05268 Pg C) because the forest only covers 10.83% of the vegetated lands. A large area of the country covered by grassland (41.74%) has the highest annual NPP (0.07693 Pg C). While the dominant rangeland type is the sparse vegetation (45.76%), its annual NPP was only 0.02717 Pg C ( Table 2). To further differentiate the NPP trends of different vegetation types, we performed a statistical analysis to derive annual NPP values per unit area. Annual NPP per square meter for different vegetation types differed greatly, and forest showed the highest annual NPP with an average value of 635.97 g C m −2 y −1 during the last 16 years ( Figure 6). Wetland vegetation had a slightly lower NPP than the forest, with an average value of 525.07 g C m −2 y −1 . The grassland NPP was about 241 g C m −2 y −1 , while shrubland and sparse vegetation were almost lower than 100 g C m −2 y −1 from 2003 to 2018. Normalized NPP by land cover areas showed increment by 0.27, 5.34, 1.57, and 0. g C m −2 y −1 in grassland, forest, wetland, and sparse vegetation, respectively; while N in shrubland vegetation decreased by -0.48 g C m −2 y −1 (Figure 6). However these N trends were not statistically significant, and the NPP fluctuated due to the changes in c Normalized NPP by land cover areas showed increment by 0.27, 5.34, 1.57, and 0.11 g C m −2 y −1 in grassland, forest, wetland, and sparse vegetation, respectively; while NPP in shrubland vegetation decreased by -0.48 g C m −2 y −1 (Figure 6). However these NPP trends were not statistically significant, and the NPP fluctuated due to the changes in climatic parameters. For example, grassland and sparse vegetation NPP were lower during the dry years (2007/2009/2010/2017), reflecting these two vegetation types are more sensitive to climatic parameters. All vegetation types show depressed NPP in 2009 and 2010, which mainly correspond with less precipitation and higher temperature, while net radiation was average. In contrast, high NPP values were observed in 2014 for all vegetation types, which was also the wettest year of the study period.

Seasonal NPP Variations by Vegetation Types
Furthermore, the statistical analysis was performed to understand seasonal variation in two contrasting (dry and wet) years for the five land cover types because seasonality is a dominant characteristic of most ecosystems. There was an apparent seasonal vegetation NPP pattern depicted in Figure 7: low and stable from October to April of each year, and high in summer, which coincided with the growing season from May to September. A large variation in the growing season NPP was found between dry and wet years, particularly in grassland, sparse vegetation and shrubland. This further verifies that these vegetation types are susceptible to climatic parameters than others. On the other hand, the wetland has less NPP difference among dry and wet years followed by forest, which can be attributed to soil and vegetation water content. Among other vegetation types, sparse vegetation was not following the general NPP trend of the growing season in the dry year, while it soared when NPP in other land cover types started declining (Figure 7). Since the growing season lasts from May to September, we focused only on the growing season (121st to 273rd day of the year) for the relationships between daily NPP and climatic parameters. At a regional scale, temperature showed the highest correlation with daily NPP (r = 0.904), followed by precipitation (r = 0.770) and net radiation (r = 0.745).

The Spatiotemporal Correlation between NPP and the Driving Parameters
Since the growing season lasts from May to September, we focused only on the growing season (121st to 273rd day of the year) for the relationships between daily NPP and climatic parameters. At a regional scale, temperature showed the highest correlation with daily NPP (r = 0.904), followed by precipitation (r = 0.770) and net radiation (r = 0.745). Some studies revealed that precipitation was an important parameter influencing NPP than temperature [48,49], but we found the temperature was the dominant parameter on NPP in the entire Mongolian rangelands compared to precipitation and net radiation. Different temporal scales of previous studies (mean growing season NPP) and the present study (daily NPP) might explain this discrepant concern.
In terms of NPP over all vegetation types, the temperature was also a dominant climatic parameter, followed by precipitation except in forest and grassland, where net radiation remains the second dominant parameter (Table 3). Though different vegetation NPP responses to precipitation diverse, precipitation played a crucial role after temperature to determine NPP in grassland, shrubland, and sparse vegetation. In the forest and wetland, net radiation was the second dominant climatic parameter after temperature. Because water stress is less significant in the forest and wetland since they have enough water content, and it needs a suitable temperature for grow vegetation. Although the temporal correlation between daily NPP and climatic parameters was positive, the spatial correlation analysis indicated negative correlations for temperature, precipitation, and net radiation, which account for 28.4%, 24.5%, and 67.32% of the vegetated lands, respectively. Annual NPP was positively correlated with temperature over central parts of Mongolia, while it was moderately correlated with net radiation in the northeast (Figure 8a,c). This could be attributed to the increased precipitation and decreased temperature in particular regions (Appendix A Figure A1a,c), where annual NPP was showed an increasing trend (Figure 4b). This is primarily attributed to the fact that these regions experience a lower mean temperature and high precipitation, which was also found by previous studies [24,47]. A similar finding was obtained over northern China, where a positive correlation between NPP and net radiation and temperature was observed [88].
Meanwhile, the correlation between NPP and precipitation was negative over northwestern parts (indicated by red color in Figure 8b), where temperature highly tended to increase with a decrement of the net radiation ( Figure A1b,c). These conditions could lead to plant water scarcity caused by accelerated evaporation, which would work against vegetation growth. This existing negative correlation was possibly either due to previous-year precipitation or available annual plants [89]. On the other hand, the area with a positive correlation in east Mongolia can be attributed to perennial plants.
Meanwhile, the correlation between NPP and precipitation was negative over northwestern parts (indicated by red color in Figure 8b), where temperature highly tended to increase with a decrement of the net radiation ( Figure A1c,b). These conditions could lead to plant water scarcity caused by accelerated evaporation, which would work against vegetation growth. This existing negative correlation was possibly either due to previousyear precipitation or available annual plants [89]. On the other hand, the area with a positive correlation in east Mongolia can be attributed to perennial plants.   Figure A1. Trends of (a) temperature and its significance level (p-value); (b) precipitation and its significance level (p-value); and (c) net radiation and its significance level (p-value) from 2003 to 2018.