Inter-annual Climate Variability and Vegetation Dynamic in the Upper Amur (Heilongjiang) River Basin in Northeast Asia

Long-term (1982–2013) datasets of climate variables and Normalized Difference Vegetation Index (NDVI) were collected from Climate Research Union (CRU) and GIMMS NDVI3g. By setting the NDVI values below the threshold of 0.2 as 0, NDVI_0.2 was created to eliminate the noise caused by changes of surface albedo during non-growing period. TimeSat was employed to estimate the growing season length (GSL) from the seasonal variation of NDVI. Statistical analyses were conducted to reveal the mechanisms of climate-vegetation interactions in the cold and semi-arid Upper Amur River Basin of Northeast Asia. The results showed that the regional climate change can be summarized as warming and drying. Annual mean air temperature (T) increased at a rate of 0.13 °C per decade. Annual precipitation (P) declined at a rate of 18.22 mm per decade. NDVI had an insignificantly negative trend, whereas, NDVI_0.2 displayed a significantly positive trend (MK test, p < 0.05) over the past three decades. GSL had a significantly positive rate of approximately 2.9 days per decade. Correlation analysis revealed that, NDVI was significantly correlated with amount of P, whereas, GSL was highly correlated with warmth index (WMI), accumulation of monthly T above the threshold of 5°C. Principal regression analysis revealed that the inter-annual variations of NDVI, NDVI_0.2 and GSL were mostly contributed by WMI. Spatially, NDVI in grassland was more sensitive to P, whereas, T was more important in areas of high elevation. GSL in most of the areas displayed high sensitivity to T. This study examined the different roles of climate variables in controlling the vegetation activities. Further studies are needed to reveal the impact of extended GSL on the regional water balance and the water level of regional lakes, providing the habitats for the migratory birds and endangered species.


Introduction
At the background of climate change, the contributions of climatic factors, such as radiation, precipitation, and temperature, on vegetation dynamics remains unclear (Craine et al 2012, Seddon et al 2016. For example, radiation was found to be the dominant factor in the rainforest and Antarctica , Andrade et al 2018, precipitation was found to be dominant in arid and semiarid areas (Choat et al 2012), while temperature was dominant in areas of cold, high elevation and high degree of latitude (Seddon et al 2016).
The dominating factors for vegetation activities are complicated in cold and semi-arid mid-latitude regions such as northeastern Asia, where plant life cycle is likely limited by temperature and also water shortage. Previous studies of vegetation dynamics in the mid-latitude regions had concluded that enhanced vegetation activity were related to regional climate warming (Fensholt et al 2012, Cong et al 2013, Zelikova et al 2015. However, the positive effect of climate warming on vegetation dynamics was insignificant despite the pronounced increase in temperature (Garonna et al 2014). In addition, there have been weakening relationships between vegetation dynamics and temperature at continental scales (Piao et al 2014), which complicates our understanding of vegetation responses to climate change in regions characterized as both cold and arid/ semi-arid.
Furthermore, the co-linearity among climate variables can mislead the conclusions of statistical analyses and results in diversified trends and pseudo driving factors , Seddon et al 2016. For example, studies on temporal variations of vegetation activities in cold and arid/semi-arid region of northeastern Asia provided controversial conclusions. Shinoda and Nandintsetseg (2011) investigated the soil moisture and plant phenology in Mongolian grassland and highlighted the availability of water in controlling vegetation activities. In contrast, Lin et al (2015) reported the insignificant relationship between precipitation and vegetation dynamics. To provide further controversy, Mohammat et al (2013) concluded temperature variability was the dominant factor in controlling the vegetation growth in northeast Asia, which complicates our understanding of the relationship between climate change and vegetation response.
The upper stream of Amur River is a typical region in north-eastern Asia characterized as both cold and dry (Simonov and Egidarev 2017). Vegetation growth and phenological change are vulnerable to climate variability (Shinoda and Nandintsetseg 2011, Mohammat et al 2013, Lin et al 2015. Besides the controversy about the driving factors, previous studies chose political borders as their study boundaries disregarding the natural watershed and basin boundaries. As the climate change will pose a threat to both ecosystem function and the provision of ecosystem services, a detailed analysis of the relationship between climate change and vegetation response over the past decades will be helpful to further our understanding of natural processes in this region. In this study, three decades (1982 to 2013) of climate variables and vegetation dynamics in the upper Amur River basin in northeastern Asia were examined to address the following questions: (1) what are the regional patterns of climate change and vegetation dynamics in this region?(2) What are the roles of climate variables and human activities, if any, in controlling the vegetation dynamics?(3) What is the contribution of each climate variable on influencing the vegetation dynamics?

Study site
The upper stream of the Amur River is located in northeastern Asia and the watershed includes territories of three countries: Russia, Mongolia, and China (figure 1). The climate of Amur basin is influenced by East Asian monsoon and characterized as warm and moderate humid summer alternated with cold and dry winter. Herbaceous plants, such as Stipa krylovii, Stipa baicalensi, and Leymus chinensis are dominant in this area (Simonov and Egidarev 2017). Based on the long-term climatic observation from Manzhouli station (49°36′0′N, 117°25′59.99′E) located approximately in the center of the study area, the annual mean temperature in the watershed is −0.6°C with −23.3°C in coldest month and 19.6°C in warmest month, and the mean total annual precipitation is 303.4 mm with 88.2% delivered during the summer of the monsoon season (June-September). The aridity index, the ratio of mean annual precipitation to mean annual potential evapotranspiration proposed by United Nations Environment Programme (UNEP), is 0.49 in Manzhouli station. According to Köppen classification of climates (Köppen 1936), the study area is situated between cold semi-arid climate (Bsk) and Monsoon-influenced warm-summer humid continental climate (Dwb) and Monsoon-influenced subarctic climate (Dwc).
This region was characterized as scarce of water but high ecological values (Simonov and Egidarev 2017). It possesses tremendous biodiversity with at least 600 species of vascular plants, 30 fish species, at least 303 bird species and about 35 mammal species. Grassland, rivers, lakes and wetlands provide a variety of habitats and niches for endangered species. As a result, the Amur River basin was identified as one of the world's 200 most valuable wilderness places. Hulun Lake (covering approximately 2,339 km 2 ) located in the center of upper stream of Amur, is enrolled in the Wetlands of International Importance (United Nations' Ramsar List).
For this study, we obtained the watershed boundary of Amur River (Heilongjiang) from HydroSHEDS (https://hydrosheds.cr.usgs.gov/) and intercepted the watershed boundary with approximately 1000 m elevation contour extracted from ETOPO1 model (Amante and Eakins 2009). This incepted sub-basin area is 687,675 km 2 with low population density of<20 individuals per km 2 (Global Human Settlement Map, European Commission, http://ghsl.jrc.ec.europa.eu/index.php). There were no change in pixels corresponding to build-up area from 2001 to 2013 (supplementary table S1 is available online at stacks.iop.org/ ERC/2/061003/mmedia) and change in goat density from 1990-2008 in Mongolia territory was considered to be negligible (Liu et al 2013).

Data source 2.2.1. NDVI data
The Normalized Difference Vegetation Index (NDVI) has been used extensively for reflections of vegetation dynamics and plant phenology (Pinzon and Tucker 2014). The Global Vegetation Inventory Modeling and Mapping Studies (GIMMS) NDVI3g dataset  was used to reveal the spatial and temporal distribution of vegetation dynamics (https://ecocast.arc.nasa.gov/). The GIMMS NDVI3g dataset was normalized for sensor calibration loss, orbital drift, and atmospheric effects (e.g. volcanic eruptions). We followed He et al (2017) method (removing data that was flagged 3 to 7-corresponding to bad data quality) to obtain the long-term NDVI images with a spatial resolution of 1/12°at 15-day intervals.
As the threshold of 0.2 in NDVI is widely considered representative for the onset/offset of plant growing period , Liang et al 2016, we developed an index of NDVI_0.2 by setting the pixel value below the threshold of 0.2 to zero. NDVI_0.2 can be more effective in representing the vegetation growth and dynamic as a result of the influence of noise, such as snow cover during winter, is removed.

Climate data
The Climatic Research Unit (CRU: University of East Anglia Climatic Research Unit, 2008, http://cru.uea.ac. uk/data) processed scheme harmonized station-based observation data sources to generate reliable estimates of global climate with spatial interpolation at 0.5°grid resolution (Mitchell and Jones 2005). The dataset CRU TS was taken for climate analysis as it provided time series of several climate variables with high accuracy (Harris et al 2014). We extracted 7 climate variables, including the Palmer Drought Severity Index (PDSI) (Palmer 1965) calculated by a physical water balance model, from grids corresponding to the study area at monthly interval from 1982 to 2013. The CRU-extracted data at Manzhouli site (49°36′0′N, 117°25′59.99′E) were validated by ground meteorological observations collected from China Meteorological Data Service Center (http:// nmic.cn/).
The PDSI is the most prominent index of meteorological drought (Dai 2011). It incorporates antecedent precipitation, moisture supply, and moisture demand into a hydrological accounting system, with the intent to measure the cumulative departure (relative to local mean conditions) in atmospheric moisture supply and demand at the surface.
Warmth index (WMI) and cold index (CDI) are considered as the key bio-climatic factors that influence plant growth, species distribution and terrestrial net primary production (Kira 1945, Wen et al 2018. In this study, WMI and CDI were prepared for each pixel using equations (1) and (2), which count the annual sum of positive and negative differences between monthly means and 5°C. As the vegetation activity is mostly influenced by coldness before the growing season (You et al 2019), the CDI in this study is counting the negative temperature difference (T−5°C) between two consecutive summers.

Growing season length
The TIMESAT program was employed in this study to smooth the raw data and extract seasonality parameters (Eklundh and Jönsson 2015). For each pixel, the median filtering method was applied to remove spikes and outliers of NDVI. Then, Savitzky-Golay method was applied to generate a smoothed NDVI time series, from which seasonality parameters were extracted (Savitzky andGolay 1964, Eklundh andJönsson 2015).
For each pixel, time series of 768 anomalies (32 year×12 month×2 images/month) was extracted. All of the settings of parameters followed the previous studies (Chen et al 2004, He et al 2017. The start of growing season (SOS) is the time for which the left edge has increased to 50% of the seasonal amplitude measured from the left minimum level. The time for the end of the season (EOS) is the time for which the right edge has decreased to 50% of the seasonal amplitude measured from the right minimum level. The growing season length (GSL) is the time from the start to the end of the growing season.
To validate the GSL extracted by TIMESAT, we conducted a comparison between the TIMESAT extracted GSL and the available phenological records reported at Hulunbeier station (49°19′N, 119°55′E) (Zhou and Liu 2017). Good agreement suggests the TIMESAT is applicable in this study (table S2; figure S1).

Human influence index
As the regional climate change was combined with human activity, a comparison between the spatial distribution of NDVI trend and human activity is helpful in understanding the role of human activity in the regional climate-vegetation relationship. The Global Human Influence Index (HII) Dataset was collected from Socioeconomic Data and Applications Center, NASA's Earth Observing System Data and Information System (EOSDIS)-Hosted by CIESIN at Columbia University (http://sedac.ciesin.columbia.edu/wildareas/) (Columbia Wildlife Conservation Society and Center for International Earth Science Information Network 2005). The HII was created from nine global data layers covering human population pressure (population density), human land use and infrastructure (built-up areas, nighttime lights, land use/land cover), and human access (coastlines, roads, railroads, navigable rivers).
Atmospheric CO 2 concentrations at Mauna Loa observatory was obtained from National Oceanic and Atmospheric Administration Earth System Research Laboratory (NOAA ESRL; http://esrl.noaa.gov/gmd/ ccgg/trends). We ignored latitudinal differences in CO 2 concentration as these are small compared to the signal of inter annual variation. Table 1 displayed all of the data sources and abbreviations.

Statistical analysis 2.3.1. Trend analysis
The trend of each variable was analyzed using the Mann-Kendall test (Mann 1945, Kendall 1975, McLeod et al 2012. A pre-whitening procedure (MK-TFPW) was applied to reduce the influence of autocorrelation on the significance of Mann-Kendall test results (Mann 1945, Yue andWang 2002). Positive tau value from MK-TFPW test indicates an increasing trend, whereas, negative tau represents a decreasing trend.

Change point detection
To detect change points in these time series, the Mann-Whitney-Pettit test was conducted (Pettitt 1979). Each variable was standardized to z-score with long-term mean and standard deviation. A low pass filter (10 years) of the standardized time series was calculated to reveal the temporal pattern of time series fluctuation.

Correlation analysis
Correlation analyses between climatic variables and vegetation activities were conducted to reveal the driving factors. Then partial correlation analyses were also conducted to reveal the relative importance among the comparable climate variables.

Principal regression analysis
The climate variables were de-aggregated to match GIMMS NDVI3g both spatially and temporally. In each pixel across the study area, all time series were transformed to z-score anomalies using long-term mean and standard deviation, and then principal regression analyses were conducted following the procedure reported by Seddon et al (2016). To present the results spatially, the relative importances of the climate variables were plotted using color composition. Figure 2 and table 2 display the 32-year trends of climatic variables and vegetation activities from 1982 to 2013. Annual T increased insignificantly with a slope of 0.13°C per decade (p=0.36), TMN increased insignificantly with a slope of 0.09°C per decade (p=0.50) and TMX also increased insignificantly with a slope of 0.17°C per decade (p=0.28), indicating the increased diurnal range of air temperature at the background of climate warming. The increased seasonal range of temperature was revealed by the integrated effects of increased WMI and decreased CDI.

Climate change and trend analyses
The seasonal fluctuation of PET follows the of variation of T. The regional average annual PET is 702.50±32.75 mm which is about 2 times higher than annual P (349.58±65.11 mm). From 1982 to 2013, annual P insignificantly declined at a rate of 18.22 mm per decade (p=0.14), whereas, annual PET significantly increased at a rate of 17.35 mm per decade (p=0.04). CLD and PDSI displayed significant negative trend, suggesting the increased solar radiation and aridity. Overall, the regional climate change was characterized as warming and drying.
Results of the Mann-Whitney-Pettit test of climate variables and atmospheric CO 2 concentration are displayed in table 2. A significant change point in P occurred around the year 1998, when P shifted from a moderate (insignificant) trend of 4.961 mm per decade (p=0.87) to a great (significant) slope of 78.52 mm per decade (p=0.02). Similarly, this significant change point (around the year 1998) was also detected in humidityrelated climate variables (PET and PDSI). Temperature variation (T, TMN and TMX) displayed an insignificant change point (around the year 1988) over the past decades.

Vegetation dynamics
The long term average NDVI is 0.35 and has an insignificantly negative linear trend over the past three decades. Whereas, NDVI_0.2 has a significantly positive MK trend (table 2). The average GSL was 135.2 days and significantly increased with a slope of 2.94 days per decade (p=0.00). Figure 3 displayed the spatial distribution of averaged NDVI and GSL. Spatially, NDVI in the eastern region (forest) displayed low inter-annual stand deviation, whereas, NDVI in grassland of the middle and western subregions displayed considerable inter-annual variation ( figure 3(b)). The pixels of positive NDVI trend were mostly distributed in the middle of the region. Spatially, the average annual GSL in high elevated (forest) areas were higher than the low land (steppe) located in the middle of the region. The positive GSL trend was mostly distributed in the eastern sub-region, whereas, the negative GSL trend was mostly distributed in the western part.
The spatial distribution of human influence did not overlap with negative NDVI trend (figure 4), suggesting the negative NDVI trend can't be explained by human activities exclusively.

Correlation analyses
Correlation analyses revealed that T, TMX, TMN PET, CLD, CDI and WMI were of insignificant relationship with NDVI (table 3), indicating the vegetation growth was insensitive to climate warming. As climate warming has an effect of drying the air, NDVI_0.2 displayed significantly negative correlations with T, TMX, TMN, WMI and PET. Contrast to the insignificant and negative roles of temperature regimes, both NDVI and NDVI_0.2 had significantly positive correlations with P (p=0.04) and PDSI (p=0.02), suggesting the vegetation growth was mostly sensitive to variation of moisture availability. The GSL displayed significantly positive correlations with PET (p=0.04), TMX (p=0.04, Kendall method) and WMI (p=0.00) but significantly negative correlations with P (p=0.02), PDSI (p=0.01), suggesting the GSL was driven differently from the patterns of NDVI and NDVI_0.2. Partial correlation revealed that GSL was mostly sensitive to WMI, whereas, NDVI and NDVI_0.2 displayed high sensitivities to both of PDSI and WMI.
The atmospheric CO 2 concentration displayed insignificant correlations with NDVI and NDVI_0.2 (table 3), suggesting the limited effect of carbon fertilization on plant growth. However, the significant correlation between CO 2 and GSL suggests a considerable influence of increased atmospheric CO 2 concentration on plant phenological changes. The atmospheric CO 2 concentration displayed a persistent increase that differs from the inter-annual variation of climate variables. Partial correlation analysis was conducted to compare the importances of increased atmospheric CO 2 and WMI, which were significantly correlated with GSL. The result showed that, the partial correlation between CO 2 and GSL is 0.150 (controlling WMI, p=0.420), whereas, the partial correlation between WMI and GSL is 0.379 (controlling CO 2 , p=0.036). Therefore, WMI weigh higher importance than CO 2 in controlling the inter-annual variation of GSL.
These results of correlation analyses were confirmed by low pass filters (10 years) of their standardized time series (figure 2). The fluctuation of standardized P had a similar temporal pattern with that of NDVI. And the temporal fluctuation of standardized WMI had a similar pattern with that of GSL.

Principal regression analyses
Principal analysis of the climate variables revealed that the first 3 principals could explain 84.4% of the total variance. Table 3 displayed the Pearson correlation coefficients between climate variables and the first 3 principal components (PC1, PC2 and PC3). The first principal (PC1) was highly correlated with temperature regime (T, TMN, TMX), the second principal (PC2) was highly correlated with moisture conditions (P, PDSI), and the third principal was highly correlated with CLD, suggesting the regional climate change can be summarized as the changes of T, P and CLD. These three principal components represented the temperature, moisture and radiation conditions which were widely considered as driving factors for vegetation activities.
Principal regression analysis of NDVI and GSL against T, P and CLD were plotted in figure 5. Spatially, NDVI in the low land (middle and southern part) was mostly sensitive to P, whereas, NDVI in higher elevation was mostly sensitive to T. GSL in this region was mostly influenced by T.
In order to reveal the importance of each climate variable, 6 principal components (totally explain 98% of inter annual variation) was extracted and then acted as independent variables. The principal regression analysis showed that WMI played a primary role in determining vegetation activities (NDVI, NDVI_0.2 and GSL), whereas, PDSI ranked the 2nd (table 3).

Characteristics of climate change
In this study, air temperature showed a moderate-positive trend of 0.13°C per decade, which was lower than the global average level (more than 0.7°C for the period 1986-2016) (Team et al 2014). The increased diurnal range of temperature revealed in this study differs from the widely reported decreasing of diurnal temperature range (Team et al 2014). The decreased P and PDSI combined with the increased PET over the past decades suggests the regional climate drying, which has been reported in central and northeastern Asia (Huang et al 2012, Chuai et al 2013, Xia et al 2014, Eckert et al 2015, Deng and Chen 2017. As a result of the climate variables were highly correlated, the significantly decreased CLD might be caused by the decreased P. The significantly increased PET might be influenced by the combined effects of the increased T and the decreased CLD (increased solar insolation). And the increased WMI coincided with the recently reported significant summer warming (Chu et al 2019, Wu et al 2020).  Table 3. Correlation (partial) coefficients between climatic variables (T, TMX, TMN, CDI, WMI, P, PET, PDSI and CLD) and vegetation characteristics (NDVI, NDVI_0.2 and GSL). The first 3 principal components (PC1, PC2 and PC3), explained 85% of total variance, were indicators of temperature, aridity and radiation. The PCA Coefficients are the multiple linear regression coefficients multiplied by the loadings of climate variables.

Pearson Coefficients
Spearman's rho Partial Correlation ( Contrary to the climate variables and vegetation dynamics that displayed high inter annual-variations, atmospheric CO 2 displayed persist increasing trend over the past decades. As a result, the correlation between the NDVI and atmospheric CO 2 concentration was weak (table 3). And the GSL displayed lower partial correlation coefficient than WMI, suggesting the limited effects of carbon fertilization on the regional vegetation dynamics.
Time series analyses revealed a significant change point around the year 1998, which was also reported by studies in the North Korea

Trends and controlling factors of vegetation dynamics
This study revealed that NDVI had no obvious trend over the entire study period . The insignificant and slightly negative linear trend of NDVI (table 2)  In term of the dominance of climate factors on controlling vegetation dynamics, the correlation analyses revealed that temperature changes were less important than precipitation. These findings coincided with the conclusion that vegetation dynamics of grassland were mostly related to changes in precipitation, rather than temperature  of accumulated temperature (above the threshold of 5°C) rather than instantaneous temperature (Jiang et al 2015, Piao et al 2015. As a result of CO 2 is the substrate for leaves photosynthesis, elevated atmospheric CO 2 have a significant influence on vegetation productivity and carbon assimilation (Keenan et al 2013). Contrary to the CO 2 fertilization effect that facilitates biomass accumulation and greening of vegetation cover (Cheng et al 2017), this study revealed the weak correlation between atmospheric CO 2 concentration and variations of NDVI and NDVI_0.2, suggesting the insignificant role of increased CO 2 on vegetation dynamics. However, increased CO 2 decreases stomatal conductance, reducing water loss through leaves photosynthesis and increasing water use efficiency particularly for water-limited environments, where drought stress is the main limitation on plant growth (Swann et al 2016). As a result, the increased atmospheric CO 2 concentration facilitated the vegetation phenological changes.
Therefore, it is concluded that the positive GSL trend might be the combined effects of the increased WMI and atmospheric CO 2 concentration. And the significantly positive trend of NDVI_0.2 was a consequence of the significantly extended GSL.

Conclusion
This study investigated the relationship between climate variables and vegetation dynamics. Time series analyses revealed that CO 2 , PET, and WMI displayed significant positive trends, whereas, P, CLD, AI and PDSI displayed significantly negative trends. As a result, the regional climate change can be characterized as climate warming and drying.
Correlation analyses revealed that drought stress weighs high importance on the inter-annual variations of NDVI and NDVI_0.2. A positive trend of temperature did not result in significant greening of vegetation over the past three decades. However, climate warming did have effects in extending the GSL and facilitating the plant growth. The increased atmospheric CO 2 concentration played an negligible role in facilitating plant growth. However, it played an active role in extending the vegetation's growing period.
Spatial distribution of controlling factors coincides with the spatial distribution of vegetation types. Forested land is sensitive to temperature, whereas, meadow grassland is sensitive to drought stress. Overall, climate warming, combined with the increasing atmospheric CO 2 concentration, has an effect of extending the plant growing period. Further studies are needed to reveal the impacts of the asymmetric warming and extended plant growing period on the regional carbon/water flux and carbon/water balance.