Spatial heterogeneity of climate variation and vegetation response for Arctic and high-elevation regions from 2001–2018

In recent decades, an amplification of warming in Arctic and high-elevation regions has been widely observed, along with a general enhancement of vegetation growth. However, driven by variability in controlling factors and complex mechanisms, climate and vegetation changes can be highly heterogeneous in space and time. In this study, an analysis is performed separating Arctic and Tibetan Plateau (TP) vegetated areas into various units according to a map of terrestrial ecoregions. The most recent variations of heat, moisture, and vegetation growth (MODIS Normalized Difference Vegetation Index) are evaluated over 2001–2018. Relationships among the climate and vegetation variables are assessed. Six distinct change patterns are identified: (1) synchronized increase of day and night temperature and precipitation during April to October coinciding with strong vegetation greening, (2) profound warming with no change in precipitation and vegetation, (3) an increase of summer temperature and vegetation with a negative latitudinal gradient, (4) greening under increased precipitation without warming, (5) browning not likely being driven by climate, (6) warming only during the nighttime and moderately enhanced vegetation growth. It is demonstrated that vegetation growth in the Arctic and TP is largely controlled by nighttime temperature and precipitation, as opposed to daytime temperature. The exception is the Canadian Arctic, where greening is directly related to summer daytime warming, and a contrasting relationship is observed on the TP. The underlying causes of these patterns are discussed, relating them to multiple mechanisms reported in the literature. These findings may help to further understand the changing Arctic and high-elevation climate and its effects on vegetation growth.


Introduction
With global climate change, cold areas such as the Arctic and high-elevation regions have experienced amplified warming, known as the Arctic Amplification (AA) and High-elevation Amplification (HA) respectively (Wang et al 2016). These two phenomena are characterized by faster warming rates than the rest of the world (Wang et al 2014, Dai et al 2019, although the specific magnitudes can vary with different locations (Boeke and Taylor, 2018). Studies have been conducted to understand the underlying drivers and mechanisms of AA and HA from multiple aspects, including elevated atmospheric CO 2 forcing (Stuecker et al 2018), snow-and icealbedo feedbacks (Screen and Ian, 2010, Taylor et al 2013, Hernández Henríquez et al 2015, and surface, cloud, and water vapor radiative feedbacks (Schweiger et al 2008, Goosse et al 2018. However, the specific consequences of amplified warming still remain uncertain. Along with the rapidly changing climate, vegetation in AA and HA regions generally exhibits greening trends, which gradually alters the landscape and environment of these cold areas (Elmendorf et al 2012a). Evidence for a greening of tundra ecosystems in the Arctic have been widely observed as increases in satellite-derived vegetation indices (e.g. Jia et al 2009, Myers-Smith et al 2011, Miles and Esau, 2016, found to be consistent with biogeochemical processes (Rocha et al 2018). Tundra greening at more local scales has been detected with in situ observations (Elmendorf et al 2012b, Emmerton et al 2016. Relationships between vegetation greenness in high-latitude ecosystems and potential drivers such as temperature (Keenan and Riley, 2018), sea ice decline (Bhatt et al 2010, Dutrieux et al 2012, snow cover dynamics (Pedersen et al 2018), extreme warming events (Pedersen et al 2018), summer sea level pressure (Bhatt et al. 2013), and anthropogenic factors (Zeng et al 2013) have been addressed for different regions and at diverse scales. However, a comprehensive understanding of the drivers of the spatial heterogeneity of vegetation change in the Arctic is still elusive. As a representative HA region, known as the 'third pole' of the Earth, the Tibetan Plateau (TP) has been reported to have enhanced vegetation growth, yet clear conclusions are lacking (Song et al 2011, Shen et al 2013, Wang et al 2017. Climate controls of vegetation dynamics on the TP have recently begun to be addressed (Hua and Wang, 2018).
As mentioned above, specific changes of climate and vegetation in AA and HA regions can be driven by different mechanisms and therefore are highly heterogeneous (Elmendorf et al 2012b). Greening trends and phenology changes were compared between North America and Eurasia by Beck and Goetz (2011) and Zeng et al (2011). They found that North America exhibited a stronger trend of enhanced vegetation growth with an earlier start of the growing season than Eurasia during the first decade of the 21st century. The spatial heterogeneity of vegetation response to climate change was investigated along two latitudinal gradients within the North Slope of Alaska (Jia et al 2006), along environmental gradients in Greenland (Karami et al 2017), and over the full latitudinal extent of Arctic tundra in North America (Alaskan North Slope and western Canadian Archipelago ; contrasting change patterns with respect to latitude were identified. Additional studies were conducted for the circumpolar Arctic (Bhatt et al 2010, Rocha et al 2018, the Siberian Arctic tundra (Dutrieux et al 2012), the Yamal Peninsula , and across the TP (Hua and Wang, 2018). However, the spatio-temporal patterns of variability and trends in vegetation and climate can be complex, and our understanding is still limited (Miles andEsau 2016, Myers-Smith et al (2019)).
To obtain a general pattern of climate and vegetation variation over AA and HA areas, and attribute controls to some of the heterogeneities of change, we conducted a study in Arctic and TP regions based on a terrestrial ecoregion map. The Arctic and TP are both underlain by permafrost. They have similar temperature, humidity, and landscapes, and have been widely reported as getting warmer and greener in the past decades. By incorporating the TP into the analyses, we intend to compare how the climate and ecosystem in these two cold regions of different causes (i.e. high latitude and high altitude) change under the background of global warming. The similarities and differences between their change patterns may help us locate the driving factors of the amplification phenomenon. The most recent version of the Moderate Resolution Imaging Spectroradiometer (MODIS) and Global Land Data Assimilation System (GLDAS) products were used to derive variables from 2001-2018 with reasonable spatial resolutions for this analysis (500 m for the vegetation, 1 km for the temperature variables, and 0.25°for precipitation variables). Compared with the widely used AVHRR GIMMS3g dataset, MODIS provides more resolved data, albeit spanning a shorter time period (Bhatt et al 2010, Dutrieux et al 2012; however, its product archives enable a robust trend analysis of almost two decades. Coarseresolution analyses are suitable for regions such as the Arctic and TP, as they have generally not been extensively affected by large-scale disturbances (Raynolds et al 2011). By separating the two areas into numerous subregions (ecoregions) of similar natural terrestrial communities, the heterogeneity of change patterns and distinct relationships between vegetation and climate variables can be evaluated. The potential driving factors of each pattern can then be discussed. These patterns may help to reveal additional clues about the mechanisms of environmental change in AA and HA regions, and contribute to better forecasts for these systems.

Study area and data sources
Our study focuses on characteristic cold region ecosystems of the Arctic and TP including tundra, and alpine grass and shrub landscapes. To identify the boundaries of these areas, we employ the terrestrial ecoregion map published by the World Wildlife Fund (WWF), in which 867 ecoregions are defined by biomes and biogeography. The map is produced to accurately reflect the complex distribution of eco-climates and the Earth's natural vegetation communities (Olson et al 2001). Twenty-six ecoregions that cover typical vegetation areas of the Arctic and TP are selected, with details listed in table A1 and their distributions presented in figure 1. Note that this selection excludes ecoregions on small islands and peninsulas, which are not necessarily representative of the broader geographic scope of the study. Twenty-one ecoregions are chosen for the Arctic, which are all in the tundra biome. Ecoregions 1-2 represent the upper and lower parts of Kalaallit Nunaat Arctic tundra located in the coastal area of Greenland. Ecoregions 3-5 denote the low, middle, and high Arctic tundra respectively in Canada. Ecoregions 6-14 are distributed mainly throughout Alaska, while ecoregions 15-21 are in Siberia. On the TP, five ecoregions (22-26) are chosen, which encompass the extent of the plateau and are all in the montane grasslands and shrublands biome.
To analyze the change of climate and vegetation growth, six parameters which are easy to acquire through remote sensing techniques were selected, including April to October land surface temperatures during the local daytime and nighttime (AO-LSTD and AO-LSTN respectively), the summer warmth index (SWI), April to October precipitation (AO-P), annual precipitation (AP), and the annual maximum Normalized Difference Vegetation Index (Max-NDVI). The five climatic parameters are widely used to evaluate the heat and moisture status of vegetation growth. The period from April to October is commonly used to study northern vegetation activities (Piao et al 2014). The summer here refers to months where mean temperatures are above freezing (>0°C). These parameters are computed and aggregated for each chosen ecoregion from 2001 to 2018 using four datasets with relatively high spatial and temporal resolutions (table 1). The processing is performed within the Google Earth Engine platform (GEE, code.earthengine.google.com). For more details about the datasets, please refer to the GEE explorer (explorer.earthengine.google.com) and the corresponding ImageCollection IDs.

Parameter computation
Parameters were chosen to reflect the primary controlling variables of heat and moisture conditions. Temperature indicators including AO-LSTD, AO-LSTN, and SWI are derived from the MODIS 8-day combined land surface temperature (LST) product with a resolution of 1 km (http://doi.org/10.5067/MODIS/ MOD11A2.006). The product is generated using the generalized split-window algorithm based on thermal infrared bands, with LST_Day and LST_Night bands retrieved from the day and night algorithms respectively based on MODIS observations. Mean monthly temperatures in April-October are summed to generate AO-LSTD and AO-LSTN, and mean monthly temperatures greater than 0°C are summed to calculate the SWI, for each pixel in each year. Pixels with insufficient quality per the corresponding quality indicator (QC) band are excluded. Finally, the AO-LSTD, AO-LSTN, and SWI are averaged over the remaining pixels for each ecoregion.
Moisture indicators including AO-P and AP are obtained using total precipitation simulated in the GLDAS-2.1 system (https://ldas.gsfc.nasa.gov/gldas/) with a resolution of 0.25°. The system is designed to generate optimal fields of land surface states and fluxes, through advanced modeling and data assimilation with satellite and ground-based observational data. The monthly mean precipitation is calculated and then summed over the entire year for AP, and summed during April-October for AO-P. The spatially averaged values of both parameters are computed for each ecoregion.  The Max-NDVI is selected to denote the peak in vegetation aboveground biomass and photosynthetic capacity. The variable is robust to seasonal data gaps, which is particularly important due to frequent cloud and snow cover in polar and high-elevation areas (Raynolds et al 2012). The parameter is derived from the MODIS 16-day combined NDVI and Normalized Difference Water Index (NDWI) datasets with a resolution of 0.5 km. Both products are generated from the MODIS version 6 MCD43A4 surface reflectance composites (https:// lpdaac.usgs.gov/products/mcd43a4v006/). NDWI is obtained with the Near-IR band and a second IR band. Similar to NDVI, it ranges from −1.0 to 1.0, with high values indicating larger likelihoods of water. First, the maximum NDVI value is computed in each year for each pixel. Then, a threshold of <0.1 is used to remove bare areas and those covered with snow and ice for the entire year. NDWI values are then checked for the remaining pixels, and a threshold of >0.7 is used to mask out surface water. Finally, the mean Max-NDVI is obtained using the remainder of the pixels for each ecoregion.
Note that the Max-NDVI only represents the peak status of vegetation, whereas changes such as those to growing season length are not considered. This could lead to an incomplete assessment of vegetation change. Similar problems might also exist for the climatic data, although five parameters were used to capture diverse aspects of heat and water status. The data processing was implemented with the original spatial resolution of each dataset, with different quality control strategies to maximize data preservation yet avoid mixed and contaminated pixels. Because values of all the parameters were aggregated for each ecoregion in the final analysis, resampling was not performed despite the datasets of different resolutions.

Data analysis
We employed the Mann-Kendall test, a non-parametric method that does not assume any particular distribution to the dataset, to evaluate if a linear trend exists for the change of selected parameters during 2001-2018 (Mann 1945). Mann-Kendall is not sensitive to missing data or outliers, not affected by the length of the time series, and robust to irregular spacing of time points. Therefore, it is widely applied in remote sensing and field measurement data analyses, which often suffer from data gaps (Wang et al 2017). All Mann-Kendall tests in this study assumed a significance level of 0.1.
When a monotonic trend is detected, the Theil-Sen estimator (Sen 1968), a nonparametric technique for linear trend estimation (El-Shaarawi and Piegorsch, 2002), is used to calculate the slope. This method is designed for robust linear regression that uses the median of all the slopes between paired values as the final slope. By using this strategy, the methodology is insensitive to outliers and can significantly outperform simple linear regression methods, such as the least square fit for skewed and heteroskedastic data. This will largely ensure a reliable slope estimation even when the data quality is low.
To assess the relationships between vegetation growth and multiple climatic parameters, the partial correlation method (Baba et al 2004) was used. This method measures the degree of association between two variables, while removing the effect of a set of controlling variables to avoid misleading information caused by correlations among independent variables.

Climate trends
Trends for all selected parameters are shown in figure 2. We observed that for both heat and moisture parameters, there were no significant negative trends, which indicates that the major trend of climate change in AA and HA areas is warming with more precipitation from 2001-2018. Only two ecoregions exhibit significantly increasing trends in SWI, which are the low Arctic tundra (slope=0.73°C months/year) and middle Arctic tundra (slope=0.40°C months/year). However, neither of these two ecoregions show significant increases in AO-LSTD and AO-LSTN, suggesting that warming in these ecoregions may actually be a rearrangement of the seasonal distribution of temperatures (i.e. the warmer summers may co-occur with cooler springs or autumns). Compared with SWI, more ecoregions are found to have increased AO-LSTD and AO-LSTN with slopes ranging from 0.40°C/year to 1.52°C/year, which may be attributed to warming in spring and/or fall. All of ecoregions in Alaska and Siberia display positive trends in both AO-LSTD and AO-LSTN, although only half of these trends are significant. The slopes for day and night tend to be similar for the same ecoregion. The Northeast Siberian coastal tundra (Ecoregion 18) only displays a significant increase in AO-LSTN, while the Northwest Russian-Novaya Zemlya tundra (Ecoregion 21) only exhibits increased AO-LSTD; however, we can see from table A2 that trends for their other AO temperature indicator are similar but less significant (with p-values just over 0.1). In contrast, Central and east TP (Ecoregion 22-24) only show significant increasing trends in AO-LSTN, a particularly interesting distinction. A small region in the western TP experienced an increase in precipitation, which may have offset the trend of warming with less solar radiation and more evapotranspiration, as both day and night temperatures were unchanged. Unlike temperature, the patterns of AP and AO-P are similar. Most regions in Alaska and Greenland show increasing trends for both parameters, while slight differences exist for a few costal ecoregions in Alaska, which may have active sea-land moisture interactions.

The trend of vegetation growth and its relationships with climate variation
Strong increasing trends in Max-NDVI ranging from 0.0016 to 0.0032 per year can be found on both Alaska and the east Siberia, approximately between 60°N and 70°N (figure 4). The Canadian Arctic tundra and the northwest TP show moderate levels of greening. It is interesting to notice that the slopes for low, middle, and high Arctic tundra are 0.0015, 0.0011, and 0.0008 per year respectively, essentially a decrease in response with increasing latitude. On Greenland, the high and low segments of the Kalaallit Nunaat Arctic tundra show opposite trends of Max-NDVI variation. The low Arctic ecoregion in Greenland is the only one that displays a significant browning trend (slope=−0.0016/year). In west Siberia and southeast TP, the changes of Max-NDVI are not significant. Greening was the trend for more than half of ecoregions in AA and HA areas over the past 18 years, however other situations still exist.
The relationships between vegetation growth and climate conditions were analyzed with two groups of parameters using a partial correlation analysis. The first group included Max-NDVI, SWI and AP reflecting summer warmth and annual moisture conditions respectively. In comparison, the second group contains Max-NDVI, AO-LSTD, AO-LSTN, and AO-P. The climatic factors denote heat and moisture status from April through October which is the period widely considered in the literatures to study the relationships between vegetation and climate for the northern hemisphere. Since the time above freezing can vary significantly for ecoregions of Arctic and TP, this group using the uniform period can be important to compare different AA and  Correlations are marked to be nonsignificant when the corresponding p-value0.1. Partial correlation coefficients were also computed using all six selected parameters. However, because patterns are hard to interpret with too many factors, the results are not displayed.
HA areas as well. Results are presented in figure 3, in which the correlation is marked to be nonsignificant when the corresponding p-value0.1. Detailed statistics can be found in tables A3 and A4.
For most selected ecoregions, except those in the Canadian Arctic and west Siberia, the Max-NDVI is negatively correlated with AO-LSTD, but positively correlated with AO-LSTN and precipitation ( figure 3). Therefore, when greening and warming occur concurrently in AA and HA areas, the increased temperature during the daytime may not be a driving factor for the enhancement of vegetation growth. Instead, the greening trend may be mainly controlled by warming during nighttime and increased precipitation. Only the middle Arctic tundra shows a positive relationship between Max-NDVI and SWI (r=0.54, p<0.05), meaning that the greening trend in this area is closely related with summer warming. The situation is opposite for the majority of the TP (r=−0.57to −0.63, p<0.05). Vegetation in the central TP is sensitive to AP, but not to AO-P, denoting that winter precipitation could play a dominant role in vegetation growth, possibly related to the dormancy requirements of some vegetation in the winter. In Kalaallit Nunaat high Arctic tundra, the vegetation growth correlates positively with both AP and AO-P. In the only browning ecoregion (Kalaallit Nunaat low Arctic tundra), vegetation does not seem to respond to any of the selected climate parameters. This suggests that the decrease of vegetation growth in this area may relate to other factors, such as wild fire and human activities. In west Siberia, the Max-NDVI also presents no significant correlation with any climatic parameter.

Different patterns of climate change and the vegetation responses
Dynamics of the selected Arctic and TP ecoregions are displayed in figure 4. From the analyses above, we can essentially divide these ecoregions into six groups with different change patterns (A to F in table 2).
Siberia can be separated into two regions (east and west), where the eastern part has similar dynamics to those of Alaska (Group A); both are coastal regions around the Bering Strait. Warming in this group occurs concurrently during the day and night from April through October, but not specifically during the summer; this could be related to spring or autumn warming. Precipitation was enhanced in most of the ecoregions during the same period, and the vegetation growth was enhanced mainly in the eastern part to the largest extent among all groups.  West Siberia (Group B) shows a strong increase in temperature (April-October), whereas precipitation did not increase concomitantly. This group exhibits no significant change of vegetation growth. The absence of greening in this area may be the major reason for the low overall greening trend in Eurasian tundra compared with North America (Bhatt et al 2010, Beck and Goetz 2011, Zeng et al. 2011. In the Canadian Arctic tundra (Group C), warming occurs along a negative latitudinal gradient (i.e. greatest warming in lowest latitudes), which is consistent with the pattern of snow cover variation reported by Hernández Henríquez et al 2015. Greening trends can be seen along a similar gradient (i.e. greater greening with decreasing latitude) which is also mentioned by Epstein et al. (2012). The precipitation of this group has not changed significantly and does not correlate with vegetation growth, suggesting that the vegetation in this area is mainly responding to summer temperatures.
The TP can be split into two parts as well. The northwest TP and north Greenland (Group D) present no warming but a significant increase in precipitation. The Southern Greenland (Group E) presents a similar trend of climate change, but reduced vegetation growth which seems to be an exception among all the selected ecoregions.
The central and eastern parts of the TP (Group F) only show increased LST at night, which is in contrast to Groups A and B. The vegetation growth is negatively correlated with SWI, which contrasts with Group C. This demonstrates different interactive mechanisms between vegetation and climate between AA and HA areas. Greening in this group only occurs in the most northern ecoregion, which has the most profound increase in AO-LSTN among the three ecoregions in the group.

Causes and consequences of different change patterns 4.2.1. Climate effects on vegetation
For Groups A and B, the significant increase of temperature in April-October may relate to severe sea ice loss (Bhatt et al 2010, Dutrieux et al 2012, which potentially led to amplified warming especially in the spring through ice-albedo feedback mechanisms. The synchronized increase of temperature and precipitation of Group A produced favorable conditions for vegetation growth. In contrast, the vegetation growth of Group B may be limited by water stress due to profound warming and no change in precipitation. For Group C, greening may relate to enhanced vegetation activity in response to the increase in extreme warm days in the Arctic (Piao et al 2014), since this area only exhibited increased summer temperatures. The observed negative latitudinal gradient of change may relate to mechanisms discussed in Räisänen (2008), which suggests that an increase in air temperature over regions that are not extremely cold can convert some component of what would otherwise be snowfall into rain, and also increase the rate of snow melt. These mechanisms led to faster and earlier snow melt in summer for areas at lower latitudes, which could also accentuate warming through snow-albedo feedbacks.
For Group D, the atmosphere may get warmer and hold more water vapor which causes a significant increase of precipitation in this area (Callaghan et al 2011). This leads to a greening trend. However, the dramatically increased precipitation can provide a cooling effect on the land surface, and consequently the land surface temperature did not change significantly.
For Group E, the south Greenland (Kalaallit Nunaat Low Arctic tundra) is browning with no clear climate drivers, which may imply other influences (e.g. disturbances). Hua and Wang (2018) reported that NDVI and precipitation has become less correlated, whereas an enhanced influence of temperature on vegetation growth was observed in TP, especially for Group F. According to the results of our study, only the northwest TP (in Group D) shows a significant correlation between Max-NDVI and AO-P, whereas almost the entire TP shows a positive correlation between Max-NDVI and AP (correlation coefficients ranging from 0.45 to 0.76, similar to the coefficients of the relationships between Max-NDVI and temperature parameters). This suggests that precipitation plays an important role in the vegetation growth on the TP; however, we cannot affirm whether this relationship is weakening or not due to the limited length of our time series.

Vegetation response feedbacks to climate
With regard to vegetation effects on amplified warming, Shen et al (2015) illustrated two opposing feedbacks of enhanced vegetation growth on atmospheric warming. On one hand, greening can decrease albedo and lead to a further increase of near-surface temperatures, mainly in the spring. On the other hand, greening enhances evapotranspiration (ET), an atmospheric cooling process, especially during the daytime in summer. The relative strengths of these opposing effects may determine the overall temperature feedback. We can see that TP shows possible ET-induced cooling in summer, as its Max-NDVI is significantly negatively correlated with SWI (figure 3i), which is consistent with the theory that the negative feedback dominates on the TP. Although an increase of Max-NDVI is only shown in the northwest part of TP (western TP in Group D and the northern part of Group F), an extended growing season length has been observed for the rest of TP areas . The underlying cooling effect during daytime for Group F (dampened increasing trend of AO-LSTD when AO-LSTN increases monotonically) implies that vegetation growth in central and east TP could be enhanced as well, but in different ways such as the aforementioned longer growing season.

Conclusions
Over the past few decades, the Arctic and the TP have both experienced profound warming and vegetation greening, which is widely considered as a characteristic of the amplification phenomenon of climate change. However, the specific dynamics can vary significantly across different areas driven by various factors and mechanisms. In this paper, we analyzed the most recent variation of climate and vegetation in typical Arctic and High-elevation Amplification areas at the ecoregion scale, using the latest version of MODIS and GLDAS products over 2001-2018. Heat, moisture, and vegetation growth were considered. Six distinct groups of dynamics were discovered. Alaska and east Siberia are characterized with an increase of day and night temperature and precipitation along with significant greening mainly in the eastern part. West Siberia exhibits profound warming without a significant change in either precipitation or vegetation growth. The warming of these two groups is possibly dominated by the sea-ice-albedo feedback, and the different vegetation conditions may relate to their precipitation variation. The Canadian Arctic shows greening and warming in summer with a negative latitudinal gradient (i.e. greatest warming and greening at lowest latitudes). The temperature is potentially affected by the snow-albedo feedback, and an increase of extreme hot days in summer could be the major cause of greening in this group especially around 60°−75°N. Northern Greenland and the northwest TP present no warming but a significant increase in precipitation with enhanced vegetation growth. Potential warming of this group is possibly dampened with the more frequent precipitation. The browning trend in the southern part of Greenland has no clear climatic drivers and may be attributed to disturbance and human activities. Central and eastern TP only show increased LST at night, which can be a result of the cooling effect induced by enhanced vegetation evaporation during the daytime. However, in this group, a greening trend is only observed for the northernmost ecoregion, which indicates that vegetation growth in the other ecoregions may be controlled by different aspects of temperature. Although additional work is still needed to fully understand these patterns, this research provides an improved understanding of the heterogeneity of climate and vegetation variation in AA and HA areas, providing new insights into change mechanisms in these regions that are highly sensitive to global change.  Table A2. The p-values and slopes of trends for the variation of SWI, AO-LATD, AO-LSTN, AP, AO-P, and Max-NDVI. The significant trends with p<0.1 are shaded. Table A3. The partial correlation coefficients and p-values between Max-NDVI, SWI, and AP for selected ecoregions. The significant correlations with p<0.1 are shaded.