Climate Dynamics of the Spatiotemporal Changes of Vegetation NDVI in Northern China from 1982 to 2015

: As an important part of a terrestrial ecosystem, vegetation plays an important role in the global carbon-water cycle and energy ﬂow. Based on the Global Inventory Monitoring and Modeling System (GIMMS) third generation of Normalized Difference Vegetation Index (NDVI3g), meteorological station data, climate reanalysis data, and land cover data, this study analyzed the climate dynamics of the spatiotemporal variations of vegetation NDVI in northern China from 1982 to 2015. The results showed that growth season NDVI (NDVI gs ) increased signiﬁcantly at 0.006/10a ( p < 0.01) in 1982–2015 on the regional scale. The period from 1982 to 2015 was divided into three periods: the NDVI gs increased by 0.026/10a ( p < 0.01) in 1982–1990, decreased by − 0.002/10a ( p > 0.1) in 1990–2006, and then increased by 0.021/10a ( p < 0.01) during 2006–2015. On the pixel scale, the increases in NDVI gs during 1982–2015, 1982–1990, 1990–2006, and 2006–2015 accounted for 74.64%, 85.34%, 48.14%, and 68.78% of the total area, respectively. In general, the dominant climate drivers of vegetation growth had gradually switched from solar radiation, temperature, and precipitation (1982–1990) to precipitation and temperature (1990–2015). For woodland, high coverage grassland, medium coverage grassland, low coverage grassland, the dominant climate drivers had changed from temperature and solar radiation, solar radiation and precipitation, precipitation and solar radiation, solar radiation to precipitation and solar radiation, precipitation, precipitation and temperature, temperature and precipitation. The areas controlled by precipitation increased signiﬁcantly, mainly distributed in arid, sub-arid, and sub-humid areas. The dominant climate drivers for vegetation growth in the plateau climate zone or high-altitude area changed from solar radiation to temperature and precipitation, and then to temperature, while in cold temperate zone, changed from temperature to solar radiation. These results are helpful to understand the climate dynamics of vegetation growth, and have important guiding signiﬁcance for vegetation protection and restoration in the context of global climate change.


Introduction
Climate change has a significant impact on the structure and function of global ecosystem [1,2]. The increase in temperature and the decrease in precipitation lead to serious drought, which will be more serious and frequent in the 21st century [3], especially in the middle latitude and arid areas [4].
As an important component of a terrestrial ecosystem [5], the dynamic changes of vegetation can reshape the global carbon and water cycle and energy flow [6]. The climate dynamics of vegetation change has become a hot topic in the field of ecology, geography, and climatology [4,7]. It is also helpful to reveal the response of vegetation biophysics to climate change [1].
The rapid development of remote sensing technology provides massive data freely available for studying global long-term vegetation change [8],vegetation dynamic monitoring [9], ecological environment assessment [10], and land cover classification [11]. The Normalized Difference Vegetation Index (NDVI) is usually used as the evaluation index in monitoring vegetation growth using remote sensing satellite data [12], because NDVI is closely related to the photosynthetic active radiation absorbed by photosynthesizing tissues [13], and it is a good representative of photosynthetic capacity [14]. Generally, the increase or decrease in NDVI indicates the enhancement or weakening of vegetation growth activity [14]. Therefore, NDVI can be used as an effective index to study the response of vegetation to climate change [9,15].
Temperature and precipitation are closely related to vegetation growth [16], and are usually selected as the dominating climate factors to analyze the response of vegetation to climate change [17]. Generally, solar radiation is also an important climate factor, because it can also control vegetation growth through spectral composition, light intensity, sunshine time, etc. [18,19]; however, related studies are mainly focused on some small-scale areas [20,21] or global scale [22]. Many studies on vegetation dynamic changes only consider the monotonic trend, assuming that the vegetation change trend remains unchanged throughout the study period [14], and the typical simple linear regression method is used to obtain the linear trend [23], which may lead to a wrong interpretation of the vegetation growth. As a common nonlinear analysis method, piecewise linear regression can better explain the changing process of vegetation growth [24]. The correlation coefficient or partial correlation coefficient is usually used to determine the degree of influence for climate factors on vegetation growth over the study period. However, the complex response of vegetation growth to the change of climate factors in different spatiotemporal scales [25] may lead to change different dominant climatic factors to vegetation growth. Moreover relevant researches are mainly concentrated on the Qinghai Tibet Plateau [26], the Loess Plateau [27], or the global scale [14].
Northern China covers a complex ecological environment composed of four climatic regions and five temperature zones, which is the key gathering place of soil and water conservation, vegetation protection, and restoration [28]. Therefore, many ecological projects have been implemented to improve vegetation coverage in China [29], such as "Grain-for-Green Project" [30], "Natural Forest Protection Project", "Three North Shelter Forest Region" [17], etc. Climate dynamics of vegetation change in northern China are mainly concentrated in Inner Mongolia [31], Loess Plateau [32], Northwestern China [33,34] and Qinghai Plateau [35,36]. There are few reports about the effect of solar radiation on vegetation growth in northern China, and it is not clear whether the relationship between vegetation dynamics and climate change has changed. Therefore, it is necessary to analyze the dominant climate factors of vegetation growth in different geographical regions, and investigate whether the dominant climate factors have changed from 1982 to 2015, and, finally, obtain the spatial pattern of dominant climate factors.
Some prior studies chose the land cover data in a certain period to analyze the relationship between vegetation and climate change in the whole period [37][38][39][40]. However, this ignores the possibility that land cover types may have changed during the study period. Some studies have shown that with the implementation of ecological engineering, the spatiotemporal characteristics of land cover have changed significantly [41,42]. Constant land cover data in the whole period will bring uncertainty to the study of the relationship between vegetation and climate change. Therefore, it is necessary to obtain unchanged land cover types from 1982 to 2015. Compared with grassland and woodland, cropland is more actively affected by human activities, which cannot reflect the response of natural vegetation to climate change. Therefore, the purpose of this study is to: (1) characterize the temporal trends and spatial patterns of climate factors and vegetation (especially grassland and woodland) in different periods of northern China from 1982 to 2015; (2)

Study Area
The study area (73 • 37 -125 • 57 E and 31 • 41 -53 • 18 N) is located in north China, including the administrative provinces of Shanxi, Hebei, Inner Mongolia, Gansu, Shaanxi, Qinghai, Ningxia, and Xinjiang. The altitude range of the study area is −156 to 7347 m (Figure 1a), and the total area is about 4,497,787 km 2 . The complex climate conditions can be divided into humid, sub-humid, sub-arid, and arid regions according to precipitation. According to the latitudes of the earth, it can be divided into subtropical, warm temperate, middle temperate, cold temperate, and plateau climate zones [43]. So, the study area is summarized as 11 different climate zones ( Figure 1a and Table 1). Based on multiphase land cover data, the spatial distribution of grassland and woodland unchanged in 1980-2015 was obtained (in Section 2.2.3) (Figure 1b).  The third generation of NDVI product created by the Global Inventory Monitoring and Modeling System (GIMMS NDVI3g) is the latest version of the long-time series NDVI data set (https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/). It covers the period from 1982 to 2015, with 8 km spatial resolution and 15 days temporal resolution, and is widely used to monitor vegetation changes at global or regional scales [44]. The 15-day temporal resolution GIMMS NDVI3g were aggregated into monthly NDVI data by Maximum Value Composite (MVC) [45]. Finally, the NDVI in growing season (NDVI gs ) was obtained by averaging the monthly NDVI values from April to October [46,47].

Meteorological Datasets
The daily 829 climate station data provided by China Meteorological Data Service Center (http://data.cma.cn/) was used. Combined with the values of several days before and after, linear interpolation was used to fill in the missing or outlier values. If there were too many consecutive blank values, the station would be removed. Finally, daily data from each station was merged into monthly data. Among the stations, 38 meteorological stations were randomly selected to verify the accuracy. The monthly raster data at a 8 km spatial resolution was obtained by interpolating the data from neighboring weather stations with ANUSPLIN 4.2 software [48]. The accuracy evaluation showed that the monthly raster data was significantly correlated with the corresponding validating weather data (p < 0.01). We choose it for interpolation because this method has good performance and has been widely used [32,49].
Monthly solar radiation data with a spatial resolution of 0.25 degrees were obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) (https://www. ecmwf.int/) ERA5. ERA5 was the fifth generation ECMWF reanalysis data for the global climate and weather covering the past 4 to 7 decades, which replaces the ERA-Interim reanalysis. To match the spatial resolution of NDVI, the thin-plate spline method was used to interpolate solar radiation data into the 8 km spatial resolution with Digital Elevation Model(DEM), which was collected from the NASA Shuttle Radar Topographic Mission, as the concomitant variation [50].
Finally, we obtained the temperature (refers to mean temperature), precipitation (refers to cumulative precipitation) and solar radiation (refers to mean solar radiation) in the growing season of vegetation.

Other Geospatial Ancillary Data
The land cover data with a 1 km spatial resolution in 1980,1990,1995,2000,2005,2010, and 2015 were provided by Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC, http://www.resdc.cn), which were widely used in regional scale related research [20,51]. Based on remote sensing images (e.g., Landsat and HJ-charge coupled device (CCD)) and ground survey data, the land cover data of China were generated by visual interpretation. The data includes six first-class classification and 25 s-class classification, and the classification accuracy was above 94.3% and 91.2%, respectively [51]. Using the nearest neighbor resampling method, the spatial resolution of seven phases of land cover data was resampled to 8 km to match that of GIMMS NDVI3g. The types of land cover mainly include woodland, cropland, grassland, water, built-up land, and unused land. The description of the land cover classification system can be found on the website mentioned above. The sub-classification of woodland (WL) includes forestland, shrubbery, sparse woodland, and other woodland, and grassland is composed of high coverage grassland (HCG), medium coverage grassland (MCG), and low coverage grassland (LCG). Using seven phases of land cover data, we derived the spatial distribution of unchanged grassland and woodland areas during 1980-2015 (Figure 1b), in order to minimize the potential impact of land cover change on the relationship between NDVI and climate factors. The WL only accounted for 19.30% of the study area (referring to the total area of woodland and grassland), the sub-classification of WL was merged into a single land cover type. HCG, MCG, and LCG accounted for 21.98%, 26.82%, and 31.90% of the total area respectively. Finally, four land cover types were used to study the response of vegetation to climate change.
Climate zone data was derived from RESDC, which was based on the division of temperature zone and precipitation.

Sen Median Trend Analysis and Mann-Kendall Test
Sen [52] median trend analysis and the Mann-Kendall test [53,54] were widely used in meteorological and hydrological research [55,56]. Sen median trend analysis has good ability to avoid outliers or measurement errors, and is not interfered by outliers [57,58]. Meanwhile, it is more suitable than the linear regression method for studying the trend of vegetation change [59]. The calculation formula is as follows: where ρ is the slope of time series data; when ]ρ > 0, it reflects that NDVI shows an increasing trend; otherwise, a decreasing trend occurs. x i and x j are the time series data. The Mann-Kendall test was used to determine the significance of the NDVI trend [30]. When |Z| > Z 1−a/2 it indicates that the trend of NDVI changes significantly. Z 1−a/2 is the corresponding value of the standard normal distribution function at the level of α. The trend of NDVI was divided into six levels according to ρ and |Z|: very significant increase (ρ ≥ 0, |Z| > 2.58, p < 0.01), very significant decrease (ρ < 0, |Z| > 2.58, p < 0.01); significant increase (ρ ≥ 0, |Z| > 1.96, p < 0.05), significant decrease (ρ < 0, |Z| > 1.96, p < 0.05), nonsignificant increase (ρ ≥ 0, |Z| ≤ 1.96, p > 0.05), and non-significant decrease (ρ < 0, |Z| ≤ 1.96, p > 0.05).

Piecewise Linear Regression Analysis
Due to the influence of climate change and human activities, the trend of NDVI usually has obvious stage (i.e., turning point) [60]. The overall change trend of vegetation in the study period cannot truly reflect the changing characteristics of NDVI in different stages. Piecewise linear regression model is very powerful in extracting the turning point (TP) of NDVI time series, and has been widely used to evaluate the change trend of hydrology [61], NDVI, and climate over time [31,62]. In a single TP case, the calculation formula is as follows: When long time series data contain more than one TP [63]. The formula is as follows: where y is the dependent variable; t is the independent variable; β 0 is the intercept; β 0 and β 1 + β 2 are the slopes of before and after the TP, respectively; α i are the estimated ith TPs of time series data, and ε is the residual random error; n is the number of TP in the time series data. Least squares linear regression was used to estimate the above parameters to minimize the sum of squares of fitting residuals, and p-value < 0.05 was considered significant.

Partial Correlation Analysis and Coefficient of Variation
Precipitation, temperature and solar radiation are the main climate factors that directly affect vegetation growth [19]; however, the three climatic factors interact each other. Through the second-order partial correlation coefficient [64][65][66], the influence of single climate factor on vegetation growth was calculated, and the other two climate factors were used as control variables. The formula is as follows: where R 12,34 refers to the partial correlation coefficient (PCC) of variables 1 and 2, and variables 3 and 4 are the control variables; R 12,3 refers to the partial correlation coefficient of variables 1 and 2, and variable 3 is the control variable; R 12 refers to the correlation coefficient of variables 1 and 2. Similarly, other variables have similar meanings. Variable 1 refers to NDVI, and variables 2, 3, and 4 represent three climatic factors, respectively. The PCC between NDVI and each climatic factor is calculated by the above formula. The t-test was used to estimate the significance of the calculated partial correlation coefficient in the case of p < 0.05 or 0.01. The coefficient of variation (CV) can reflect the low or high stability of vegetation NDVI time series with higher or lower CV value [67]. The formula is: where CV, σ, NDV I refer to the coefficient of variation, the standard deviation, the mean value of NDVI time series. According to previous studies [33,[68][69][70], the stability of vegetation was divided into three levels: high stability (CV ≤ 0.08), medium stability (0.08 < CV ≤ 0.12), and low stability (CV > 0.12).
To better explain the overall process of this study, the flow chart is shown in Figure 2.  Figure 3 shows the spatial distribution of the multi-year mean NDVI gs and CV of NDVI gs from 1982 to 2015. The NDVI gs of the study area ranged from 0.01 to 0.86 with a spatial average of 0.37. Overall, the spatial pattern of NDVI gs is that it gradually decreased from southeast to northwest, and had obvious spatial heterogeneity, which was basically consistent with the spatial distribution of land cover types ( Figure 1b). The NDVI gs of WL distributed in the southeast of the study area was the largest with a spatial average of 0.59. The coverage area of grassland was the largest, and the NDVI gs average values of HCG, MCG, and LCG were 0.42, 0.34, and 0.22, respectively. Overall, the NDVI gs of vegetation had good stability, indicated by the spatial average of CV being 0.07. According to the corresponding percentages of different CV grades in Figure 3b, high stability region accounted for 63.94%, medium stability region accounted for 26.26%, low stability region accounted for 9.80%, especially in the central (e.g., the Loess plateau) and western parts of the study area. According to Figures 1b and 3b, the inter-annual stability of WL was the best, with the average CV of 0.05, followed by HCG (CV = 0.07). The CVs of MCG and LCG were similar, slightly higher than 0.08.

Temporal Variations of NDVI gs and Climatic Factors in Different Periods
At the regional scale, the TPs of NDVI gs in the study area from 1982 to 2015 were identified by piecewise linear regression, which occurred in 1990 and 2006 ( Figure 4). From 1982 to 2015, NDVI gs showed an obviously significant increasing trend, with a slope of 0.006/10a (p < 0.01), indicating that the vegetation in the past 34 years was constantly improving. During 1982-2015, the temperature also showed a very significant increasing trend, while the precipitation and solar radiation showed a non-significant increasing trend, indicating that the climate conditions were gradually warm and humid. According to the two TPs of 1990 and 2006, the NDVI gs change in 1982-2015 could be divided into three periods: 1982-1990, 1990-2006, and 2006-2015 (Figure 4a). In 1982-1990, the slope of NDVI gs was 0.026/10a (p < 0.01), which was higher than that of NDVI gs in 1982-2015; temperature and precipitation showed an increasing trend. However, solar radiation showed a decreasing trend. In 1990-2006, the slope of NDVI gs suddenly decreased with a nonsignificant decreasing trend of −0.002/10a (p > 0.1). The temperature trend was higher than that in 1982-1990; precipitation had a decreasing trend, with a slope of −12.199 mm/10a (p > 0.1); the slope of solar radiation was 0.467 MJ/m 2 ·10a (p < 0.1). In 2006-2015, the slope of NDVI gs became positive, indicating that the vegetation growth began to improve. At the same time, the slope of precipitation also changed to positive, showing a significant increase trend. On the contrary, the slopes of temperature and solar radiation changed from positive to negative, showing a decreasing trend. From Figure 4, the NDVI gs slope is consistent with that of precipitation, positive-negative-positive change, which is opposite to that of the solar radiation slope. From 1982 to 1990, with the increase of temperature and precipitation, climate conditions become warm and humid, which was conducive to vegetation growth, and then NDVI gs increased accordingly. In 1990-2006, with the increase of temperature and solar radiation, precipitation decreased, resulting in dry climate, which was not conducive to vegetation growth, so NDVI gs decreased accordingly. In 2006-2015, with the decrease of temperature and solar radiation and the increase of precipitation, the climate tends to be cold and humid, and NDVI gs increased rapidly.

Spatial Pattern and Trend of NDVI gs and Climatic Factors in Different Periods
Sen median trend analysis and Mann-Kendall tests were used to analyze the trends and significance levels of vegetation and climate factors in different periods on the pixel scale. Figures 5 and 6, respectively, show the slope and significance level for NDVI gs and climatic factors in 1982-2015. Table 2 shows the area statistics of NDVI gs and climate factors at different significance levels. From 1982 to 2015, the trends of NDVI gs , temperature, precipitation and solar radiation were mainly positive ( Figures 5 and 6, Table 2). The area with extremely significant increase of NDVI gs (p < 0.01) was close to that with non-significant increase (p > 0.05), while the decreasing trend was mainly non-significant (p > 0.05). According to the area statistics of NDVI gs of different land covers at different significance levels (Table S1), from 1982 to 2015, there was obvious improvement in vegetation; LCG had the largest area improvement percentage of 78.08%, followed by MCG (76.28%) and WL (73.45%), and HCG (68.66%). The area of WL, MCG or LCG NDVI gs with significant increase (including p < 0.05 and p < 0.01), which was much larger than that with non-significant increase (p > 0.05), while HCG NDVI gs was opposite. The increase of temperature, precipitation, and solar radiation in the middle of the study area was good for vegetation growth. On the contrary, the decrease of precipitation and the increase of temperature and solar radiation led to the decrease of NDVI gs in the northeast of the study area, which indicated that precipitation might play a decisive role in the dynamic change of vegetation.
The spatial patterns of NDVI gs and climatic factors were significantly different in different periods (Table 2). During 1982-1990, NDVI gs , temperature and precipitation showed an increasing trend, which was mainly composed of non-significant increase (p > 0.05). However, the slope of solar radiation was mainly negative, and the area of decreasing trend accounted for 76.12% of the total area ( Table 2, Figures S1 and S2). The area with increasing trend of WL, HCG, MCG, and LCG NDVI gs was much higher than that with decreasing trend (Table S1). The area with an increasing trend of WL NDVI gs was the highest (21.39%, including p < 0.05 and 0.01), while that of LCG NDVI gs was the lowest (11.30%) (Table S1). From 1982 to 1990, 85.34% of the regions showed an increasing trend of NDVI gs , indicating that vegetation growth continued to improve, and vegetation degradation were mainly distributed in arid and sub-arid regions. The decrease of temperature was mainly located in the southeastern and western parts of the study area. The decrease of precipitation was mainly distributed in the western, southwestern, and southeastern parts of the study area. The increase of solar radiation was mainly concentrated in the northeastern B4 and the southwestern B3, while the other areas showed a decreasing trend.
From 1990 to 2006, the results clearly showed that more than half of the vegetation had been degraded ( Table 2, Figures S3 and S4). The area of NDVI gs increase decreased from 85.34% in 1982-1990 to 48.14% in 1990-2006. Over 84% of the study area, the temperature and solar radiation showed an increasing trend, and the area with increasing or decreasing precipitation was basically the same. (Table 2, Figures S3 and S4). Compared with 1982-1990, the NDVI gs changes of WL, HCG, MCG, and LCG were significantly different, and the area with increasing trend decreased sharply. The area with increasing trend of WL and HCG NDVI gs decreased from 93.24% and 86.20% in 1982-1990 to 41.05% and 40.85% in 1990-2006, both of which decreased by more than 45%. The area with increasing trend of MCG and LCG NDVI gs was slightly higher than that with decreasing trend (Table S1). Vegetation degradation was serious in the northeastern, southwestern, and northwestern parts of the study area. The temperature in most of the study area showed an increasing trend. The precipitation increased in the southwestern, western, and southern parts of the study area, and decreased in the northeastern part of the study area. Except for the plateau climate zone, solar radiation was increasing, especially in the northeast of the study area ( Figures S3 and S4).
In 2006-2015, NDVI gs and precipitation increased in more than 68% of the regions, while temperature and solar radiation decreased in more than 71% of the region, and the significant level was mainly non-significant ( Table 2, Figures S5 and S6). Compared with 1990-2006, the area proportion of the NDVI gs and precipitation with positive trends increased, while that of temperature and solar radiation decreased sharply. Compared with 1990-2006, WL, HCG, MCG, and LCG had significantly improved, and the area with an increasing trend had increased to some extent, but was still less than that in 1982-1990 (Table S1) Table 2).
The trend of NDVI was not a simple linearity, and had obvious differences before and after the TP. During 1982During -1990During , 1990During -2006

Spatial Pattern of the Dominant Climate Factors of NDVI gs
The PCC of NDVI gs with temperature, precipitation, and solar radiation were denoted by R Tem , R Pre and R SR , respectively. The impact of the three climate driving factors was obtained by linearly stretching of the absolute values of R Tem , R Pre and R SR to 0-255, and then by RGB synthesis (Figure 7, Figures S7-S9). Table 3 shows the statistics result of PCC between NDVI gs and climatic factors in different periods. According to the square of the PCC, the spatial distribution of the dominant climate factors was divided ( Figure S10). Table S2 shows the statistics results of climatic factors for different vegetation in different climatic zones.
From 1982 to 2015, the average PCC between vegetation and precipitation was the largest, followed by temperature, both were much higher than solar radiation, which was the same as the area ranking through significance level test, but the area with positive correlation between vegetation and temperature was larger than that with precipitation ( Figure 7, Table 3). The regions with R Tem less than 0 were mainly distributed in sub-humid regions, sub-arid regions, and parts of arid regions. The values of R Pre in the western B3, the northeastern B4, the northern C2 and the parts of C4 were higher (R Pre > 0.4), while regions with R Pre less than 0 were mainly distributed in A1, D2, and D3. R SR was mainly negative in plateau climatic zones, arid regions and some sub-arid regions, while was mainly positive in the humid zones, sub-humid areas. Figure 7d shows the spatial distribution of dominant climate drivers. The vegetation growth in the humid areas, plateau climate zones, and high-altitude areas were mainly affected by temperature and solar radiation. The vegetation growth was primarily affected by precipitation in arid region, sub-arid region, and part of sub-humid region (Figure 7d and Figure S10a, Table S2).
However, the dominant climate factors of vegetation growth changed significantly in different periods. From 1982 to 1990, NDVI gs was positively correlated with temperature and precipitation, with an area of more than 55%, but negatively correlated with solar radiation, with an area of more than 60%, with spatial average values of 0.098, 0.061, and −0.148, respectively ( Figure S7, Table 3). The regions with R Tem less than 0 were mainly concentrated in the eastern plateau climatic zone, B4, northwestern C2 and southern B3. The negative PCC in the central B4 was lower than −0.4, while the positive PCC in the northeast of the study area was higher (R Tem > 0.6). R Pre was mainly negative in the southern and western of the plateau climate zone, the humid and sub-humid zones. The area with R SR less than 0 was widely distributed, and the area with R SR larger than 0 was mainly distributed in the northeast of the study area. The vegetation growth in the eastern of the plateau climate zone, the sub-humid zone, and the high-altitude area was mainly affected by solar radiation. In the northeastern part of the study area, vegetation growth was mainly affected by temperature. The vegetation growth affected by precipitation was mainly distributed in the northeastern B4 and the northwestern B3, while the remaining area was affected by two or three climatic factors (Figures S7d and S10b, Table S2).
Compared with 1982-1990, the correlation between vegetation and temperature and solar radiation decreased significantly, but increased with precipitation during 1990-2006. The area with R Tem larger than 0 was basically the same, and the area with of R Pre and R SR larger than 0 significantly increased ( Figure S8, Table 3). The regions with R Tem < 0 were mainly distributed in the plateau climate zone, some arid, sub-arid, and sub-humid zones. R Pre was mainly negative in the northeast of the study area, southeast and northwest of the plateau climate zone. R SR was mainly positive in E1, the southwestern C2. The growth of vegetation in the eastern of the plateau climate zone and E1 was mainly affected by solar radiation. The vegetation growth was mainly affected by temperature in the northern part of the study area and the southern of the plateau climate zone. The vegetation growth affected by precipitation was mainly distributed in the northeastern B4 and the southwestern B3 ( Figure S8 and S10c).
Compared with 1982-1990 and 1990-2006, the average value of R Tem continued to decrease, and the sign of R SR changed from negative to positive, and R Pre continued to increase. The areas with R Pre > 0 and R SR > 0 increased gradually, while R Tem continued to decrease ( Figure S9, Table 3). R Tem was mainly positive in the plateau climate area and the northeast of the study area. The areas with R Pre < 0 were mainly distributed in the southeastern of the plateau climate zone, southeast, and northeast of the study area. The areas with R SR > 0 were mainly distributed in the northeast of the study area, southwestern of the plateau climate area. Vegetation growth was mainly controlled by temperature in the southern of the plateau climate area, some high-altitude areas, while in the northern parts of the study area, central B3, vegetation growth was mainly limited by solar radiation. The spatial distribution of vegetation growth affected by precipitation was widespread, mainly in the eastern B2-4 and C2. Other regions were mostly affected by two or three climatic factors.
According to the Figures S7-S9, the positive and negative response of vegetation to climate factors in the same region had changed significantly with time. From 1982From -1990From to 1990From -2006 and then to 2006-2015, the area with R Tem > 0 in the northeast of the study area gradually decreased, while in the plateau climate zone, it first decreased and then increased. The area with R SR > 0 increased gradually, especially in the plateau climate zone. The effects of temperature on vegetation growth in plateau climate zone were gradually intensified, and the dominant climate factor change from solar radiation to temperature. In arid, sub-arid, and sub-humid zones, the area with R Pre > 0 gradually increased, and the impact of precipitation on vegetation growth gradually strengthened ( Figures S7-S10, Table S2). Table 4 shows the area percentage of vegetation growth controlled by the three climate factors in different land cover types. Overall, vegetation growth was mainly controlled by precipitation, followed by temperature in the study area during 1982 to 2015. However, the dominant climate factors of vegetation growth also changed in different periods. During 1982-1990, solar radiation was the main climate factor, and the influence of temperature and precipitation was similar. In 1990-2006 and 2006-2015, the impact of precipitation and temperature on vegetation growth increased and the impact of solar radiation weakened, precipitation became the main influencing factor, followed by temperature. Different type vegetation had different responses to climate factors. From 1982 to 2015, HCG and MCG were firstly affected by precipitation, followed by temperature, while WL and LCG were mainly affected by temperature, followed by precipitation, and the solar radiation had the lowest impact on the growth of four vegetation. The dominant climate factors affecting the growth of four vegetation had difference in 1982-1990, 1990-2006, 2006-2015: the main climate factors were temperature and solar radiation, precipitation and temperature, precipitation and solar radiation (by the magnitude of the impact) for WL; and were solar radiation and precipitation, precipitation and temperature, precipitation and solar radiation for HCG. For MCG, they were precipitation and solar radiation, precipitation and temperature, precipitation and temperature. For LCG, they were solar radiation and temperature, precipitation and temperature, temperature and precipitation. Figure 7. Spatial distributions of the partial correlation coefficients (PCCs) between NDVI gs and temperature (a), precipitation (b), and solar radiation (c) and spatial pattern of dominant climate factors (d) in the study area, from 1982 to 2015. Note: Tem, Pre, and SR refer to temperature, precipitation, and solar radiation, respectively.

Discussion
Using GIMMS NDVI3g, meteorological station and climate reanalysis data, we evaluated vegetation changes in the northern China in recent 34 years, and analyzed the impact of climate factors on vegetation NDVI before and after the TP. Generally, the driving factors affecting vegetation growth include climate and non-climate factors, which are mainly human activities and natural disturbances [22], such as agricultural irrigation, afforestation, urbanization, pests and human disturbances [71][72][73][74]. Climate factors can provide necessary growth conditions for vegetation [22], and dominate vegetation change in large-scale areas, and restrict the geographical distribution of vegetation in a terrestrial ecosystems [32].
However, human activities can change vegetation growth in small areas in a short time, which has a significant impact on regional ecosystem functions [75,76]. At present, we only pay attention to the effects of three climatic factors on the growth of vegetation, and do not quantitatively evaluate the potential effects of human activities or natural disturbances on vegetation growth in different periods before and after the TP.
It is generally believed that climate warming has a positive impact on vegetation growth in high latitude areas of northern hemisphere, which is conducive to improving photosynthesis of vegetation and prolonging the length of growth season [77]. Similarly, Peng [78] also considered that temperature rise was the main reason for the increase of vegetation in northern China. However, the results of Wu [22] indicated that the vegetation growth at high northern latitudes were mainly affected by temperature and solar radiation from 1982 to 2008, which further confirmed the reliability of our results. On the long time series scale , the increase of vegetation greenness in high-latitude areas of northern China was mainly affected by temperature and solar radiation. However, the dominant climate factors that affect vegetation growth usually change over time, and this fact will be ignored in the overall analysis during the study period. The results show that the influence of solar radiation on vegetation growth increases gradually in the northern high latitude region, while the effect of temperature decreases gradually, which is opposite in plateau climate region. Climate factors can promote or limit vegetation growth. In the plateau climate zone and humid area, excessive precipitation is not conducive to vegetation growth, because it caused lower temperature and insufficient sunshine time to inhibit vegetation photosynthesis [79]. When solar radiation exceeds the optimal amount for vegetation growth, the extra amount of solar radiation leads to a decrease of evapotranspiration, which has a negative impact on the growth of vegetation. In the northeastern part of the study area, temperature rise and precipitation increase triggered the significantly increase of NDVI in the growing season before the 1990s. However, after the 1990s, the increase of temperature and solar radiation, and the decrease of precipitation, matched the decrease of NDVI, which is consistent with the studies of Peng [78]. After 2006, with the sufficient precipitation and decrease of temperature and solar radiation, vegetation growth has been effectively improved, because precipitation can regulate soil moisture, thus affecting the root activity and water status of vegetation [78]. In addition, when the temperature is lower than the optimum temperature of photosynthesis, temperature rise will promote vegetation growth, but when the temperature is higher than the optimum temperature, vegetation photosynthesis will be inhibited, and the consumption of respiration and nutrients will be accelerated to limit vegetation growth [80]. From 1982 to 1990, the temperature increased and the solar radiation decreased in the plateau climate zone promoted vegetation growth ( Figures S2 and S7), while the continuous increase of temperature in 1990-2006 exceeded the optimum temperature of photosynthesis and then restricted vegetation growth. The excessive temperature accelerated vegetation transpiration and the loss of soil water, partly offset the effect of precipitation increase in promoting vegetation growth. Many studies have shown that the solar radiation decreased obviously after the 1960s, and began to increase in the 1990s [81][82][83][84]. During 1982During -1990, the possible reason why the decrease of solar radiation promoted the more active growth of vegetation may be that the solar radiation is higher than what is needed to support the growth of vegetation. There is a physiological reason for such a phenomenon. With the decrease of solar radiation, the stomatal conductance and photosynthesis rate of leaves continued to increase, and NDVI would increase accordingly [85]. This negative correlation between solar radiation and vegetation growth has also been found in other studies [21,86,87].
The response of vegetation growth to climate factors shows obvious geographical heterogeneity. The response of the same vegetation to climate factors will also change in different geographical locations. During 1982-2015, the woodland growth was mainly controlled by temperature, followed by precipitation in the study area (Table 4), but those in humid areas were more affected by temperature and solar radiation. The woodland in the sub-humid zone or sub-arid was mainly controlled by temperature and precipitation, but the difference of temperature or humidity zones can change the order of the dominant climate factors (Table S2). Generally, precipitation may play a dominant role in controlling the growth of grassland vegetation [75]. The response of grassland to climate factors was slightly different in different climatic zones. Grassland in B3-4 accounted for 38.32% of the total grassland area including HCG, MCG, and LCG, and its growth was mainly controlled by precipitation. The grassland in the plateau climate zones (D2-4) accounts for 30.02% of the total grassland area, and the dominant climate factor for its growth was temperature. However, the grassland in the warm temperate zone (C2, C4) accounted for 27.77% of the total grassland area, and its growth was mainly limited by temperature and precipitation (Table S2). These results describe the spatial pattern of climate driving factors of vegetation growth in northern China, and provide a reference for understanding the response of vegetation change to climate change.
In this study, we only used the single NDVI gs as the evaluation index, and did not analyze the response of vegetation to climate factors on annual average or seasonal average NDVI. In addition, the lag effect of vegetation response to climate change was not considered. In the study, the land cover data of multiple periods are used to obtain the unchanged land cover area in the study period, which can avoid the uncertainty of the relationship between vegetation and climate change caused by a single period land cover data. We focused on the change of vegetation growth in WL, HCG, MCG, and LCG, and did not analyze the transfer of four vegetation types in different periods. The unchanged land cover data reduced the spatial distribution of grassland and woodland to a certain extent, which may cover the fact that the LCG was transferred to the MCG and HCG, or the MCG and HCG degenerated into the LCG. In the future, we will focus on the first level land cover data to avoid this situation.
In the future, we will consider the impact of terrain factors (such as altitude and slope) on the relationship between vegetation and climate change, and quantitatively evaluate the relative impact of human activities on vegetation change. At the same time, more climate data (such as soil moisture, evapotranspiration, and drought index) will be used to analyze the relationship between vegetation and climate change.
Supplementary Materials: The following are available online at https://www.mdpi.com/2072-4 292/13/2/187/s1, Figure S1: Spatial distribution of the slopes of NDVIgs (a), temperature (b), precipitation (c), and solar radiation (d) in the growing season of the study area during 1982-1990. Figure S2. Spatial distribution of the significant levels of NDVIgs (a), temperature (b), precipitation (c), and solar radiation (d) in the growing season of the study area during 1982-1990. Figure S3. Spatial distribution of the slopes of NDVIgs (a), temperature (b), precipitation (c), and solar radiation (d) in the growing season of the study area during 1990-2006. Figure S4. Spatial distribution of the significant levels of NDVIgs (a), temperature (b), precipitation (c), and solar radiation (d) in the growing season of the study area during 1990-2006. Figure S5. Spatial distribution of the slopes of NDVIgs (a), temperature (b), precipitation (c), and solar radiation (d) in the growing season of the study area during 2006-2015. Figure S6. Spatial distribution of the significant levels of NDVIgs (a), temperature (b), precipitation (c), and solar radiation (d) in the growing season of the study area during 2006-2015. Figure S7. Spatial distributions of the PCCs between NDVIgs and temperature (a), precipitation (b) and solar radiation (c) and spatial pattern of dominant climate factors (d) in the study area from 1982 to 1990. Figure S8. Spatial distributions of the PCCs between NDVIgs and temperature (a), precipitation (b) and solar radiation (c) and spatial pattern of dominant climate factors (d) in the study area from 1990 to 2006. Figure S9. Spatial distributions of the PCCs between NDVIgs and temperature (a), precipitation (b) and solar radiation (c) and spatial pattern of dominant climate factors (d) in the study area from 2006 to 2015. Figure S10. The spatial distribution of the main control areas of three climate factors in different periods. Table S1. Statistical results of NDVIgs significance level of different land cover in different periods (%). Table S2. Statistics results of climatic factors for different vegetation in different climatic zones (%).
Author Contributions: R.S. processed the data, analyzed the results, and wrote the manuscript. R.S., S.C. and H.S. revised the manuscript. All authors have read and agreed to the published version of the manuscript.