Spatiotemporal Dynamics of Vegetation Net Primary Productivity and Its Response to Climate Change in Inner Mongolia from 2002 to 2019

Abstract

Understanding vegetation dynamics and their responses to climate change are essential to enhance the carbon sequestration of the terrestrial ecosystem under global warming. Although some studies have identified that there is a close relationship between vegetation net primary productivity and climate change, it is unclear whether this response exists in ecologically fragile areas, especially in Inner Mongolia, in which multiple ecological ecotones are related to vegetation types. This study uses the Carnegie–Ames–Stanford Approach (CASA) model to estimate vegetation NPP in Inner Mongolia from 2002 to 2019 and focuses on the spatial and temporal changes of NPP of different vegetation types and their responses to three typical climate factors: precipitation, temperature, and solar radiation. The results show that the NPP estimated by the CASA model agrees well with the observed NPP (R² = 0.66, p < 0.001). The vegetation NPP in Inner Mongolia decreases gradually from northeast to southwest, and the average NPP is 223.50 gC ∙ m⁻². From 2002 to 2019, the NPP of all vegetation types trended upward, but exhibiting different rates. The vegetation types, ranked in order of decreasing NPP, are forest, cropland, grassland, and desert. The NPP response of different vegetation types to climate factors possesses significant differences. The cropland NPP and grassland NPP are mainly affected by precipitation, the desert NPP is controlled by both precipitation and solar radiation, and the forest NPP is determined by all three climate factors.

1. Introduction

At present, the carbon cycle is an important index to evaluate the effects of global climate change [1,2]. As an essential part of the terrestrial carbon cycle, vegetation net primary productivity (NPP) can effectively quantify the production capacity of the terrestrial ecosystem, which has become a hot topic in global carbon cycle research [3]. The net primary productivity (NPP) of vegetation refers to the rate at which vegetation accumulates organic matter [4]. NPP is the difference between carbon absorbed by plant photosynthesis and carbon released by plant respiration [5,6]. NPP is the key factor in judging the carbon source and sink function of the ecosystem, and in evaluating the sustainable development of the terrestrial ecosystem. With the in-depth study of global change and the carbon cycle, the estimation of regional NPP has attracted more and more attention from scholars all over the world [7]. The research on regional NPP has been identified as one of the core contents in the International Geosphere-Biosphere Programme (IGBP), Kyoto Protocol, and the Paris Agreement [8,9,10]. The carbon cycle in ecologically vulnerable areas is very sensitive to global climate change. Quantitative assessment of regional NPP can provide theoretical support for accelerating “carbon peak” and “carbon neutralization” in ecologically vulnerable areas and provide important guidance for coping with climate change and regional ecological environment protection.

Inner Mongolia (IM) is an important ecological barrier in northern China. Before the year 2000 and the growth of population and lifestyles being transformed from traditional nomadic to modern settled, land degradation caused by overgrazing was a very serious issue [11,12]. From 1985 to 1999, the degraded grassland area of the Xilin River Basin in IM reached 72% of the total area of the basin [13]. After the year 2000, to solve the problem of land degradation, a variety of ecological compensation measures were implemented, such as direct cash compensation for rest, forbidden grazing grassland, and ecological projects such as returning farmland to forest or grassland [11,14]. Up to now, the current ecosystem function has been improved to a certain extent. However, there are also some negative effects, such as the decrease in herdsmen’s income and the rise of the unemployment rate [15]. This has raised concerns about the sustainability of ecological restoration. Therefore, it is necessary to explore the evolution of different vegetation ecosystems in IM and form a reasonable land management and utilization model. This will provide a reference for the ecological restoration of other ecologically fragile areas.

Due to the coexistence of a variety of ecological ecotones, such as forest-steppe ecotone, farming-pastoral ecotone, and desert-oasis ecotone, IM has become a typical ecologically fragile area and is one of the areas most sensitive to global climate change [16,17,18]. In recent years, there have been many abnormal climate changes in Inner Mongolia, such as extreme drought, which makes NPP fluctuate very considerably in this region [19]. Many researchers have analyzed the contribution of climate factors to NPP change in IM [20,21,22,23]. However, most of them have only focused on grassland types and ignored the response of other vegetation types to climate factors. In the study of global change, the response of vegetation to climate change shows great temporal and spatial heterogeneity, and different vegetation types have different sensitivity and response characteristics to climate change [24,25,26]. Therefore, only analyzing the grassland ecosystem is not enough to understand the carbon cycle process of the regional ecosystem. The quantitative estimation and the evolution of temporal and spatial characteristics of vegetation NPP and its response analysis to different climate factors are of great significance to evaluate the carrying capacity of the terrestrial ecosystem, and formulate sustainable development policies in IM.

Although there are several global NPP products, they are not always suitable for regional-scale research. This is because global NPP products use global-specific parameters and ignore inter- and intra-regional variability [27,28]. Therefore, for regional-scale research, it is necessary to use localization parameters to accurately estimate regional NPP based on existing NPP models. The CASA model is the most representative productivity estimation model [29]. The CASA model works by considering the process of plant photosynthesis and using the concept of light energy utilization (LUE). The CASA model has been widely applied because it can well evaluate the evolution trend and spatio-temporal variability of NPP on global and regional scales [30,31]. It should be noted that the application of the original maximum LUE rate (0.389 gC ∙ MJ⁻¹) in the CASA model to different vegetation types at the regional scale will bring large deviations [32]. Therefore, we localized the maximum LUE for different vegetation types to better describe the temporal and spatial characteristics of regional NPP. The main objectives of this study are: (1) Spatial distribution pattern and influencing factors of vegetation NPP in IM; (2) Temporal variation characteristics and environmental driving forces of vegetation NPP in IM; and (3) Temporal and spatial variation of NPP of different vegetation types and its response to climate change in IM.

2. Materials and Methods

2.1. Study Area

Inner Mongolia (IM) is in the plateau region of northern China, with longitude from 97°12′ E to 126°04′ E and latitude from 37°24′ N to 53°23′ N. IM covers a total area of 1180,000 km² and has a typical temperate continental climate. The average annual precipitation and the average annual temperature varies spatially, falling from 450 mm in the northeast to 50 mm in the southwest, and rising from −1 °C in the northeast to 10 °C in the southwest, respectively. This results in the gradual transition of IM from humid and semi-humid areas to arid and semi-arid areas from east to west. Accordingly, the vegetation types are distributed from east to west as eastern forest, central grassland, and western desert. IM is one of the most important animal husbandry and agricultural production bases in China, but its terrestrial ecosystem is fragile and vulnerable to climate change and human activities. Therefore, it is vital to understand the spatiotemporal dynamic of vegetation productivity for the sustainable development of IM. The location of IM in China and the land cover map in 2019 are shown in Figure 1. The vegetation types in IM mainly include forest, grassland, and cropland. It should be noted that the desert area in the west of IM is covered with a specific proportion of shrubs and has a specific carbon sequestration capacity. Therefore, the desert is also regarded as a type of vegetation in this study. According to the 2019 MODIS land cover products (MCD12Q1, see Section 2.2.1 details), the proportions of four main vegetation types including forest, grassland, cropland, and desert in the total area of IM are 10.7%, 59.9%, 6.3%, and 22.5%, respectively.

2.2. Data

2.2.1. Remote Sensing Data

The MODIS NDVI product MOD13A1 from 2002 to 2019 were used in this study (https://lpdaac.usgs.gov/products/mod13a1v061/. Accessed on 20 November 2021). The temporal resolution of MOD13A1 is 16 days and the spatial resolution is 500 m. The monthly NDVI data for NPP calculation were obtained by the maximum values composite (MVC) method for 16-day NDVI [33]. The WGS84/Albers equal-area conic projection was adopted for all the data used in this study.

Yearly land cover data from 2002 to 2019 is derived from MODIS land cover product MCD12Q1 (https://lpdaac.usgs.gov/products/mcd12q1v006/. Accessed on 20 November 2021). MCD12Q1 has a spatial resolution of 500 m and a temporal resolution of 1 year. It uses the land cover classification system defined by the International Geosphere Biosphere Program (IGBP), utilizing a total of 17 land cover classes [34]. According to the needs of this study, we reclassified original the 17 land cover classes into six main classes, as follows: (1) evergreen needleleaf forests, evergreen broadleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, mixed forests, closed shrublands, and woody savannas are reclassified as forest; (2) open shrublands, savannas, grasslands and permanent wetlands are reclassified as grassland; (3) croplands, and cropland and natural vegetation mosaics are reclassified as cropland; (4) urban and built-up lands are reclassified as developed; (5) barren is reclassified as desert; and (6) permanent snow and ice, and water bodies are reclassified as water.

2.2.2. Meteorological Data

The meteorological data from 2002 to 2019, including monthly average temperature, monthly cumulative precipitation, and monthly total solar radiation, are derived from 46 meteorological stations in or surrounding IM, which can be downloaded from the National Meteorological Science Data Center (http://data.cma.cn. Accessed on 20 November 2021). To match the above remote sensing data both spatially and temporally, the meteorological grid images are generated by Kriging interpolation for the original meteorological data. Kriging method is a widely used interpolation method, which is characterized by low deviation and high precision [35].

2.2.3. Field Sampling Data

To validate the NPP estimated by the CASA model, we measured the biomass data of 42 grassland sample plots in IM in August 2016 and August 2018. Grassland rather than other vegetation types was selected because the above ground part of grassland can be completely collected and easily measured, while other vegetation types, such as forest and desert, need destructive logging for analysis to take place. Considering that the CASA model has been widely verified [36], we did not use other methods to obtain the observed data for the forest, cropland, and desert. Only grassland was used to verify the feasibility of the CASA model in IM. All sample plots were randomly selected and were expected to be widely distributed in IM. However, some areas, such as eastern grassland, are not allowed to be sampled due to local ecological protection policies. The locations of sampling sites are shown in Figure 1. The area of each sample plot is 1 m × 1 m The vegetation aboveground was collected in its entirety and dried to constant weight in an oven at 70 °C, and the dry weight was recorded. Then the corresponding underground biomass was obtained by multiplying the dry weight of above ground biomass by a factor of 2.8 according to the recommendation of the Intergovernmental Panel on Climate Change [37]. Finally, the carbon conversion rate was set to 0.475 to obtain the observed NPP [38].

2.3. Methods

2.3.1. NPP Estimation

In this study, the CASA model was utilized to estimate the monthly NPP of vegetation in IM. Following this, the yearly NPP from 2002 to 2019 was aggregated by monthly NPP in one year. The CASA model is a typical light use efficiency (LUE) model driven by remote sensing and meteorological data and has a robust performance in estimating long time-series NPP [29]. The NPP is determined by absorbed photosynthetically active radiation (APAR) and LUE factor (ε) in CASA model using Equation (1):

$$NPP\left(a,m\right)=APAR\left(a,m\right)\times \epsilon \left(a,m\right)$$

(1)

where $NPP\left(a,m\right)$ [gC ∙ m⁻²] is NPP of pixel $a$ in month $m$. $APAR\left(a,t\right)$ [MJ ∙ m⁻²] is photosynthetically active radiation absorbed by pixel $a$ in month $m$, and $\epsilon \left(a,m\right)$ [gC ∙ MJ⁻¹] is the actual LUE of pixel $a$ in month $m$. APAR is determined by the total solar radiation (SOL) and the fraction of absorbed photosynthetically active radiation (fPAR) using Equation (2):

$$APAR\left(a,m\right)=SOL\left(a,m\right)\times fPAR\left(a,m\right)\times 0.5$$

(2)

where $SOL\left(a,m\right)$ [MJ ∙ m⁻²] is the total solar radiation of pixel $a$ in month $m$, the constant 0.5 is the proportion of solar effective radiation (400 nm–700 nm) that can be used by vegetation in SOL. $fPAR\left(a,m\right)$ is the absorption ratio of vegetation to incident photosynthetically active radiation (PAR). fPAR has a good linear relationship with NDVI and thus can be calculated by Equation (3):

$$fPAR\left(a,m\right)=\frac{\left(NDVI\left(a,m\right)-NDV{I}_{k,min}\right)\times \left(fPA{R}_{max}-fPA{R}_{min}\right)}{NDV{I}_{k,max}-NDV{I}_{k,min}}+fPA{R}_{min}$$

(3)

where $NDVI\left(a,m\right)$ is the NDVI value of pixel $a$ in month $m$, $NDV{I}_{k,max}$ and $NDV{I}_{k,min}$ are the maximum and minimum NDVI value for the k^th land cover class, respectively. $fPA{R}_{max}$ and $fPA{R}_{min}$ are constants and are set to 0.95 and 0.001, respectively.

LUE factor (ε) is the efficiency of vegetation converting APAR into organic carbon for vegetation. It is mainly affected by temperature and precipitation and is calculated by Equation (4):

$$\epsilon \left(a,m\right)=T_1\left(a,m\right)\times T_2\left(a,m\right)\times W\left(a,m\right)\times {\epsilon }_{max}$$

(4)

where $T_1\left(a,m\right)$ and $T_2\left(a,m\right)$ are temperature stress factors, $W\left(a,m\right)$ is water stress factor, ${\epsilon }_{max}$ [gC ∙ MJ⁻¹] is the maximum LUE. Because the maximum LUE has a great impact on NPP estimation, it should be set carefully. The monthly maximum LUE of global vegetation used in the traditional CASA model is 0.389 gC ∙ MJ⁻¹. Currently, many researchers modify this value according to the specific vegetation types. Here, we set ${\epsilon }_{max}$ for each vegetation type in IM based on the research results of reference [39]. The details for calculating $T_1\left(a,m\right)$, $T_2\left(a,m\right)$, and $W\left(a,m\right)$ refer to the relevant literature [5,29,40].

2.3.2. Trend Analysis

In this study, the Theil–Sen estimator [41] was used to analyze the temporal variations of vegetation NPP in IM from 2002 to 2019. The Theil–Sen estimator was selected because it has good robustness to outliers. The inter-annual change rate of a single pixel location is the slope of the trend line in the Theil–Sen estimator equation, which is calculated by Equation (5):

$$slop{e}_a=Median \left(\frac{NP{P}_j-NP{P}_i}{j-i}\right),i,j=1, 2,\dots , N$$

(5)

where $slop{e}_a$ is the slope of the trend line in Theil–Sen estimator, $NP{P}_i$ and $NP{P}_j$ are the NPP at year $i$ and $j$ ($j$ > $i$), respectively. N is the number of years during the study period (N = 18 in this study). The slope with a negative value indicates a decreasing trend of NPP, while the slope with a positive value indicates an increasing trend of NPP. F test was used to conduct the significance test of inter-annual trend. This significance only represents the confidence level of trend change. The equation of F test is as follows:

$$F=U\times \frac{N-2}{Q}$$

(6)

$$U={\displaystyle \sum }_{i=1}^N{\left(\widehat{NP{P}_i}-\overline{NPP}\right)}^2$$

(7)

$$Q={\displaystyle \sum }_{i=1}^N{\left(NP{P}_i-\widehat{NP{P}_i}\right)}^2$$

(8)

where $NP{P}_i$ are the NPP at year $i$, $\widehat{NP{P}_i}$ is the predicted value of $NP{P}_i$, $\overline{NPP}$ is the average value of NPP during the study period, and N is the number of years. According to the test results, the change trend is divided into the following five levels: extremely significant decrease ($slop{e}_a<0, P<0.01$), significant decrease ($slop{e}_a<0, 0.01<P<0.05$), no significant change ($P>0.05$), significant increase ($slop{e}_a>0, 0.01<P<0.05$), extremely significant increase ($slop{e}_a>0, P<0.01$).

2.3.3. Partial Correlation Analysis

To analyze the impact of climate factors on NPP in IM, partial correlation coefficients [42] between annual NPP and three climate factors, e.g., annual average temperature, annual precipitation, and annual cumulative solar radiation for each pixel location from 2002 to 2019, are calculated. Partial correlation analysis is chosen because it can exclude the influence of other related variables when studying the correlation relationship between two variables. A second-order partial correlation analysis is required because four variables (NPP and three climate factors) are considered in this study. The second-order partial correlation coefficient is calculated by Equation (9):

$$R_{ya·bc}=\frac{R_{ya·b}-R_{yc·b}·R_{ac·b}}{\sqrt{\left(1-R_{yc·b}^2\right)\left(1-R_{ac·b}^2\right)}}$$

(9)

where $R_{ya·bc}$ is the second-order partial correlation coefficient between variable y and a after fixing variables b and c. $R_{ya·b}$, $R_{yc·b}$ and $R_{ac·b}$ are the first-order partial correlation coefficient between variable y and a, variable y and c, and variable a and c, respectively. The first-order partial correlation coefficient is calculated by Equation (10):

$$R_{pq·s}=\frac{R_{pq}-R_{ps}·R_{qs}}{\sqrt{\left(1-R_{ps}^2\right)\left(1-R_{qs}^2\right)}}$$

(10)

where $R_{pq·s}$ is the first-order partial correlation coefficient between variable p and q after fixing variable s. $R_{pq}$, $R_{ps}$, and $R_{qs}$ are the Pearson correlation coefficient between variable p and q, variable p and s, and variable q and s. Among the three climate factors of each pixel, the one with the largest partial correlation coefficient is regarded as the main factor. This is a reasonable result, as the larger the coefficient, the stronger the correlation [43].

3. Results

3.1. NPP Validation

We used 42 biomass data from grassland sample plots to evaluate the estimation accuracy of the CASA model. All biomass data were measured in August 2016 and August 2018 and were then converted into observed NPP. Figure 2 shows the result of correlation analysis between observed NPP and estimated NPP. We see that the estimated NPP is significantly correlated with the observed NPP (R² = 0.66, p < 0.001). This indicates that the CASA model is suitable for multi-year NPP estimation in IM.

3.2. Spatiotemporal Dynamics of Vegetation NPP

3.2.1. Spatial Distribution Pattern of Vegetation NPP

The average annual NPP and its standard deviation in IM from 2002 to 2019 was shown in Figure 3a,b, respectively. Overall, the average annual NPP shows obvious spatial heterogeneity. This is because there are significant differences in vegetation types and climate in different subregions. The average annual NPP of all vegetation types in IM from 2002 to 2019 is 223.50 gC ∙ m⁻². It is characterized by decreasing spatially from northeast to southwest. The average annual NPP of the forest in the northeast is mostly above 400 gC ∙ m⁻², while that of the desert in the west is mostly less than 100 gC ∙ m⁻². The NPP of the grassland and the cropland in central IM is mainly between 100–400 gC ∙ m⁻², and the grassland NPP is lower than the cropland NPP.

The standard deviation of the average annual NPP in IM from 2002 to 2019 also has spatial heterogeneity, but its pattern is different from the average annual NPP. The large standard deviation (>55 gC ∙ m⁻²) mainly appears in the eastern regions of IM. This indicates that the inter-annual NPP fluctuated greatly in these regions during the study period. The possible reason is that these regions are close to being of a grassland type and have experienced noticeable degradation and restoration during the study period. In contrast, the standard deviation of the western desert type is very low (<15 gC ∙ m⁻²). This can be expected because the vegetation in the desert is mostly sparse shrubs that have low but stable productivity.

3.2.2. Inter-Annual Variation of Vegetation NPP

The spatial pattern of inter-annual change rate and change trend significance of NPP in IM from 2002 to 2019 was shown in Figure 4a,b, respectively. For the inter-annual change rate of NPP, the value of most pixels is positive, which means that NPP increases in these pixel locations over time. On the contrary, some pixels with negative values indicate that NPP decreases in these pixel locations over time. The area percentage corresponding to the different change rate ranges is shown in

Table 1. The area percentage of the inter-annual change rate and change trend significance of NPP in IM from 2002 to 2019.

Change Rate (gC∙m−2∙a−2)	Percentage (%)	Trend	Percentage (%)
<−3	0.6	ESD	11.0
−3–0	15.3	SD	1.7
0–3	53.1	NSC	12.5
3–6	23.7	SI	1.9
>6	7.3	ESI	72.9

. The area with a positive change rate accounted for 84.1% of IM, and the area with a negative change rate only accounted for 15.9% of the IM. Obviously, the area with increased NPP is significantly greater than the area with decreased NPP. The region with a large NPP increase mainly appears in the east of IM, while the regions with sporadic NPP decrease appear in the western desert and middle grassland. In addition, from the significance of the inter-annual change trend, most pixels are extremely significantly increased (ESI), and only a few pixels are extremely significantly decreased (ESD). The corresponding area percentage is also shown in Table 1. The ESI accounted for 72.9%, and the ESD was only 11%. It is worth noting that the pixels with no significant change (NSC) accounted for 12.5% and are mainly concentrated in the western desert area. This indicates that the change of NPP in the desert type has no linear relationship with time. The area percentage of significant increase (SI) and significant decrease (SD) is very close, as they are 1.9% and 1.7%, respectively.

Figure 5a shows the inter-annual change of NPP for four vegetation types, e.g., forest, grassland, cropland, and desert in IM from 2002 to 2019. The range of inter-annual change for four vegetation types is completely different. The average annual NPP from high to low is forest (496.17 gC ∙ m⁻²), cropland (348.64 gC ∙ m⁻²), grassland (222.25 gC ∙ m⁻²), and desert (64.49 gC ∙ m⁻²). The annual NPP of forest, cropland, and grassland have some fluctuations, while the desert is relatively stable. The fluctuations of the forest, cropland, and grassland increased considerably after 2010. This is mainly because there were more extreme weather events, such as drought, in IM after 2010 [44]. The NPP of cropland and grassland decreased significantly in 2007 and 2017, which is due to extreme drought events in these two years. Low precipitation inhibited vegetation growth. It is worth noting that the forest NPP did not decrease significantly in 2007; this is because the drought event in 2007 was related to continuous high temperatures [45,46]. To a certain extent, the increase in temperature was beneficial to the increase in the forest NPP, as it offset the decrease in the forest NPP caused by drought. In addition, the NPP of forest and cropland peaked in 2014, which may be related to solar radiation; the solar radiation in 2014 was the highest throughout the study period [44]. The trend fitting parameters of inter-annual NPP change are listed in

Table 2. The linear fitting parameters of inter-annual NPP change trend for different vegetation types in IM from 2002 to 2019.

Types	Slope	R2	p
Forest	3.6988	0.6784	<0.01
Cropland	2.1864	0.3305	<0.05
Grassland	1.7651	0.1442	>0.05
Desert	0.0806	0.0051	>0.05

. All the change rates of the four vegetation types are positive, which indicates that NPP had an increasing trend for four vegetation types from 2002 to 2019. Among the types, the change trend of grassland and desert type are NSC (p > 0.05), while forest type and cropland type are ESI (p < 0.01) and SI (p < 0.05), respectively.

Figure 5b shows the total NPP for four vegetation types, e.g., forest, grassland, cropland, and desert in IM from 2002 to 2019. The inter-annual change of the total NPP for all vegetation in IM is significant. The highest total NPP for all vegetation was 405.24 TgC ∙ a⁻¹ in 2018 and the lowest was 311.57 TgC ∙ a⁻¹ in 2007. The inter-annual change rate of the total NPP for all vegetation was 3.27 TgC ∙ a⁻¹. Although the inter-annual change of total NPP for all vegetation fluctuated greatly, it is still ESI (R² = 0.3520, p < 0.01). Grassland contributes the most (the average percentage is 58.9%) to the total NPP for all vegetation, and the inter-annual change trend of total NPP for grassland is consistent with that for all vegetation types. This can be expected because the total area of grassland is greater than the other three types. The contribution (the average percentage is 25.1%) of forest types to the total NPP for all vegetation is lower than that of grassland, but higher than that of cropland (the average percentage is 8.8%) and desert (the average percentage is 6.8%) types. It is worth noting that the total NPP for forest, cropland, and desert types does not fluctuate significantly from year to year and was relatively stable during the study period.

3.3. The Impact of Climate Factors on Vegetation NPP

Figure 6 illustrates the second-order partial correlation coefficients between annual NPP and three climate factors, e.g., annual precipitation, annual average temperature, and annual cumulative solar radiation, at each pixel location in IM from 2002 to 2019. Moreover, the main driving factor of each pixel location is also shown. Overall, the response of annual NPP to the three climate factors has obvious spatial heterogeneity, but it can still be observed that the annual NPP at most pixel locations has a strong correlation with precipitation and solar radiation, while only a few pixel locations have a strong correlation with temperature. Specifically, Figure 6a illustrates the partial correlation coefficient between annual NPP and annual precipitation. Among them, the pixel locations with positive correlations account for 95.73% of IM, and most of them have a high value of correlation coefficients (e.g., red color indicates a positive correlation coefficient with high value). The pixel locations with negative correlations only account for 4.27% of IM. Spatially, the pixel locations with positive correlations are mainly distributed in the semi-arid, semi-humid and humid areas from the middle to the east, while the pixel locations with negative correlation are mainly located in the arid areas in the west. For the annual temperature shown in Figure 6b, the pixel locations with positive correlation account for 37.81% of IM and are mainly distributed in the northeast forest area. The pixel locations with negative correlations account for 62.19% of IM and are mainly located in the central grassland area. For annual cumulative solar radiation shown in Figure 6c, the pixel locations with positive correlations account for 82.74% of IM and are mainly distributed in the northeast forest area and the western desert area. The pixel locations with negative correlations account for 17.26% and are mainly located in the central grassland area. It is worth noting that although the proportion of positively correlated pixel location is high (82.74%), the value of the correlation coefficient is generally low (yellow-green color indicates a positive but relatively low correlation coefficient). In addition, the comprehensive analysis of the influence of the three climate factors on annual NPP is shown in Figure 6d. It can be found that precipitation is the main driving factor of the inter-annual change of NPP in IM. The pixel locations are mainly driven by precipitation, which accounts for 82% of IM. The inter-annual change of NPP in a small part of the west and northeast (accounting for 14.85 of IM) is mainly driven by solar radiation. Only the sporadic pixel locations (about 3.15% of IM) are driven by temperature.

3.4. The Response of NPP of Different Vegetation Types to Climate Factors

Table 3. Mean second-order partial correlation coefficients between annual NPP and three climate factors for the four vegetation classes.

	Forest	Grassland	Cropland	Desert
NPP_Precipitation	0.4490	0.5887	0.5212	0.2888
NPP_Temperature	0.1557	−0.1462	−0.1229	−0.0153
NPP_Solar radiation	0.3198	0.1654	0.1995	0.2884

lists the average partial correlation coefficients between the annual NPP of the four vegetation types, e.g., forest, grassland, cropland, and desert, and three climate factors. In general, all four vegetation types have the strongest correlation with precipitation, followed by solar radiation and temperature. However, the response of each vegetation type to climate factors is significantly different. There is a positive partial correlation coefficient between forest and all three climate factors, which means the forest NPP is driven by all three climate factors. The response of grassland and cropland to the three climate factors has a similar pattern, in which there are positive coefficients for precipitation and solar radiation, and a negative coefficient for temperature. It is worth noting that the partial correlation coefficient between grassland or cropland and precipitation is significantly greater than that with solar radiation. This indicates that the NPP of grassland or cropland is mainly driven by precipitation. Desert NPP has similar partial correlation coefficients with precipitation and solar radiation, which means that desert NPP is driven by both precipitation and solar radiation, but is independent of temperature due to the partial correlation coefficient being close to zero.

4. Discussion

4.1. Vegetation NPP Estimation in IM

The vegetation NPP in IM estimated by different models differs greatly. The average NPP of vegetation in IM from 2002 to 2019 estimated in this study is 223.50 gC ∙ m⁻², of which the NPP of the forest, grassland, cropland, and desert are 496.17 gC ∙ m⁻², 222.25 gC ∙ m⁻², 348.64 gC ∙ m⁻² and 64.49 gC ∙ m⁻², respectively. Compared with the existing studies, these results are reasonable. For example, for the grassland NPP, the result of this study is similar to those estimated by the same model in reference [46] and in reference [47], but lower than those estimated by the EcoC-G (ecological consumption of grassland) model in reference [48] and greater than those estimated by Synthetic model in reference [20] and by ORCHIDEE model in reference [49]. The main reasons for this diversity are as follows: (1) different study periods; (2) differences in structure, mechanism, and key parameters between models; and (3) different input data source and spatiotemporal resolution.

4.2. Climate Change and Vegetation Response in IM

Vegetation activity has changed significantly due to global warming in recent years. Many studies have shown that vegetation activities in the middle and high latitudes of the northern hemisphere, including most areas of China, have increased significantly [50,51]. The results of this study are consistent with the above results. This study indicates that the vegetation NPP in IM shows an increasing trend from 2002 to 2019. The NPP of different vegetation types has the same increasing trend but different increasing rates. The NPP of forests increases significantly, while deserts’ increases only slightly.

Most previous studies only focused on one vegetation type, that is, grassland, due to the area of grassland being the largest type [52,53]. This may lead to bias in the evaluation of regional NPP. This study found that there was no significant change trend of the inter-annual grassland NPP (although the slope of the linear fitting is positive, p > 0.05 means that it failed the significance test, see Table 2). However, the total NPP of IM has increased significantly (the slope of the linear fitting is 3.27 TgC ∙ a⁻¹, and p < 0.01, see Figure 5b), which indicates that other vegetation types, such as forest and cropland, may have a more important impact on the regional NPP than grassland. Therefore, it is necessary to consider all vegetation types to deeply understand the dynamics of regional NPP. The zonal distribution of vegetation is determined by precipitation, temperature, and solar radiation. Accordingly, the sensitivity of different vegetation types to the changes in precipitation, temperature, and solar radiation is also quite different. A large amount of water vapor evaporated from the ocean is transported to the land through atmospheric circulation, which is the main source of atmospheric precipitation on the land. The precipitation on land decreases gradually from coastal to inland, and the vegetation types also change as you go from the forest, to grassland, to arid desert. IM spans 19 longitude zones, and its vegetation types show regular longitude zonal distribution. In general, the NPP of all vegetation types in IM is positively correlated with precipitation. In other words, precipitation is the main factor affecting the vegetation NPP in IM. This conclusion of this study is consistent with previous similar studies [20,52]. Furthermore, this study also found that the sensitivity of different vegetation types to precipitation was different, and the partial correlation coefficients were grassland > cropland > forest > desert.

There is a positive correlation between the forest NPP and all three climate factors, e.g., precipitation, temperature, and solar radiation. The correlation coefficient is precipitation > solar radiation > temperature, which is consistent with the similar study in northeast China in reference [54]. The forest in IM is mainly distributed in the northeast region where there is abundant precipitation and very low temperatures. Therefore, the type of forest is mainly the evergreen coniferous forest, which is insensitive to temperature, leading to the low correlation between the forest NPP and temperature. In contrast, the study in Guangxi province, China, in reference [55], found that temperature has a more important impact on the forest NPP than precipitation. It is necessary to distinguish the study area and vegetation types when analyzing the response of the forest NPP to climate change. In addition, most studies conclude that the change of the grassland NPP is mainly affected by precipitation, which is contrary to the change of temperature [56,57]. Our study also confirmed this point. This is because grasslands in IM are mainly distributed in arid and semi-arid areas with less precipitation. When the temperature rises, the soil evaporation increases, resulting in the decrease in soil water content. This means that the soil cannot provide sufficient water for grassland growth, resulting in a decline in productivity. In addition to climatic factors, the growth status of grassland is related to grassland management measures. For example, reference [58] found that precipitation could not improve the degradation of the plant community in desert grassland of IM caused by heavy grazing. The results of reference [16] show that the implementation of environmental projects such as returning cropland to grassland and the control of overgrazing play a vital role in the increase in the grassland NPP in recent years. The response of the cropland NPP to climate change has the same pattern with grassland, but there is only a slight decrease in correlation coefficient with precipitation. This is because some croplands depend on artificial irrigation, which reduces the dependence on precipitation, such as the Hetao irrigation area in southwestern IM [59]. The desert NPP in IM is mainly affected by precipitation and solar radiation, but the inter-annual fluctuation of NPP is very small due to low vegetation coverage and poor growth status.

4.3. Limitations and Future Research Directions

The CASA model is driven by remote sensing and meteorological data. In practical application, the spatial and temporal resolution of remote sensing and meteorological data are always different. For example, the spatial resolution of MODIS NDVI (MOD13A1 product) used in this study is 500 m, while the spatial resolution of meteorological data may reach 10 km. The transformation between data with different resolutions will bring some uncertainty [60,61] to the estimation results of the CASA model. How to obtain or generate remote sensing and meteorological data with the same resolution to reduce uncertainty is one possible research direction for the future.

In addition, from 2002 to 2019, vegetation NPP in IM shows an increasing trend, which may be due mainly to two aspects. On the one hand, the increase in precipitation in IM may be the main factor leading to the increase in vegetation NPP. On the other hand, human protection of the environment may also be a reason for the increase in vegetation NPP. This study discussed the correlation between vegetation NPP and climate factors in IM but did not analyze the impact of human activities on vegetation NPP. Although relevant studies have proved that the change of NPP in IM is mainly affected by climate change [62], human activities should not be ignored as an important factor [63]. Therefore, how to quantify the impact of various human activities on vegetation NPP is another important research direction in the future.

5. Conclusions

In this study, the spatiotemporal characteristics of NPP for different vegetation types in IM from 2002 to 2019 were obtained by the CASA model, and its response to climate factors was analyzed. The results show that the vegetation NPP in IM gradually increases from south to north and from west to east. The NPP of different vegetation types has a significant difference, and the NPP of the forest, cropland, grassland, and desert are 496.17 gC ∙ m⁻², 348.64 gC ∙ m⁻², 222.25 gC ∙ m⁻², and 64.49 gC ∙ m⁻², respectively. Annual NPP of vegetation in IM has an obvious increasing trend but with some fluctuations from 2002 to 2019. A total of 84.1% of the study area increased significantly, and only 15.9% of study area decreased. The annual NPP of all vegetation types has the same increasing trend but different increasing rates. The increasing rate from high to low is 3.6988 gC ∙ m⁻²∙ a⁻¹, 2.1864 gC ∙ m⁻²∙ a⁻¹, 1.7651 gC ∙ m⁻²∙ a⁻¹, and 0.0806 gC ∙ m⁻²∙ a⁻¹ for the forest, cropland, grassland, and desert, respectively. The partial correlation between NPP and three climate factors, e.g., precipitation, temperature, and solar radiation is calculated for each pixel location. We found that precipitation is the most important impact factor for the change of vegetation NPP in IM, followed by solar radiation and temperature. The responses of different vegetation types to climate factors are significantly different. The forest NPP are jointly affected by three climate factors. The grassland NPP and cropland NPP mainly depend on precipitation, while the desert NPP is affected by both precipitation and solar radiation. The above research results are of great significance for understanding the spatiotemporal characteristics of vegetation in IM and formulating relevant vegetation protection and restoration policies.