Effect of climate and thaw depth on alpine vegetation variations at different permafrost degrading stages in the Tibetan Plateau, China

ABSTRACT

Understanding the driving forces for alpine vegetation variations at different permafrost degrading stages is important when the Tibetan Plateau is experiencing climate warming. We applied the modified Frost Number model to simulate frozen ground distributions in the Tibetan Plateau and calculated the maximum thawing depth by the Stefan approach. We classified the simulated frozen ground into three subzones: seasonal frozen ground zone, changing zone, and permafrost zone. We evaluated the effects of precipitation, air temperature, and maximum thawing depth on Normalized Difference Vegetation Index (NDVI) in the subzones across five different stages from 1982 to 2012. The results show that permafrost zone, changing zone, and seasonal frozen ground zone account for about 30.6 percent, 23.3 percent, and 46.1 percent of the study area, respectively. Over the five stages, permafrost areas decreased at fast, slow, fastest, and then slowest rate from stage1 to stage 5, and the large continuous permafrost area has been degraded into pieces. Precipitation is strongly correlated with NDVI and contributes most `to the changes of NDVI. Maximum thawing depth and particularly air temperature show a much smaller correlation and contribute less to the variation rate of NDVI. The findings will have broad applications in investigating the impact of climate and environment changes on alpine vegetation variations in the Tibetan Plateau.

KEYWORDS:

• Permafrost degradation
• maximum thawing depth
• NDVI
• Tibetan Plateau

Introduction

Permafrost, the frozen ground at low temperature below 0 °C for ≥2 years (Muller 1947), is highly sensitive to climate change and surface disturbances, especially to air temperature changes (Li et al. 2008; Yang et al. 2010). Permafrost in the Tibetan Plateau, approximately accounting for 53 percent of the plateau area, is the world’s highest altitude and largest frozen ground region in low latitude (Cheng et al. 2012). Compared with high-latitude permafrost, permafrost in the Tibetan Plateau is warmer and shallower (Wu and Zhang 2008; Wu, Zhang, and Liu 2010; Yang et al. 2010). Therefore, permafrost in the Tibetan Plateau is more sensitive to air temperature changes and it is considered the key indicator for climate and environment changes (Wang et al. 2016). During recent decades, a significant climate warming occurred on the Tibetan Plateau and the increase in temperature came earlier and faster than in other parts of China (Pan and Li 1996; Liu and Chen 2000; Wang et al. 2000; Cheng and Wu 2007; Wu and Zhang 2010; Kuang and Jiao 2016; Lu, Zhao, and Wu 2017). The increasing temperature on the Tibetan Plateau resulted in permafrost thawing, which is confirmed by both site observation data (Lemke et al. 2007; Wu and Zhang 2010) and numerical modeling (Li et al. 2008; Cheng et al. 2012; Nan et al. 2012; Nan, Huang, and Zhao 2013; Lu, Zhao, and Wu 2017). The permafrost degradation is mainly characterized by a reduction in areal extent of permafrost occurrence and active layer deepening (Jin et al. 2009; Xue et al. 2009; Wu and Zhang 2010).

Studying the changes of permafrost in the Tibetan Plateau has become an important issue. Although field observations can provide more realistic and accurate data, these data are hard to obtain because of the harsh environment. Numerical simulations make it possible to study the temporal and spatial variations of permafrost. Different models have been applied in the Tibetan Plateau, such as the Altitude model (Li and Cheng 1999), the MAGT (mean annual ground temperature) model (Nan, Shuxun, and Yongzhi 2002), the Frost Number model (Nan et al. 2012; Lu, Zhao, and Wu 2017), the TTOP (temperature at the top of permafrost) model (Smith and Riseborough 1996) and the CLM (community land model) model (Guo, Wang, and Li 2012; Guo and Wang 2013; Guo, Wang, and Wang 2017). The Frost Number model has been widely applied in many regions (Lu, Zhao, and Wu 2017). Though many models have been used to simulate the distribution of permafrost in the Tibetan Plateau, they are mainly focused on its current and future distributions. Few models provided the evolution of permafrost, and none of them outlined a degrading permafrost zone in the Tibetan Plateau. Until now, the distribution and evolution of degrading permafrost on the plateau are still unknown.

Degradation of permafrost and its subsequent effects on the ecological environment have attracted attention from many scientists. Considering that the Tibetan Plateau is a cold and arid region, some studies regard precipitation and temperature as the two main factors controlling vegetation changes (Pang, Wang, and Yang 2016; Zhang et al. 2016). However, the changes in permafrost conditions have significant impacts on ecosystems (Yang et al. 2010). Yi et al. (2011) suggested that climate warming might promote the expansion of alpine grassland in permafrost regions. Other studies suggest that permafrost degradation is responsible for the environmental deterioration and grassland degeneration (Wang et al. 2006; Yang et al. 2010; Zhou et al. 2015). The thawing of permafrost results in a lowering of the water table, thus leading to degradation of alpine grassland (Cheng and Wu 2007; Wang et al. 2008). The increase in the thickness of active layer will affect the alpine ecosystem (Wu et al. 2015), and the response of vegetation to frozen ground degradation vary spatially with the change rate of the frozen ground depth (Qin et al. 2017). In general, permafrost ecosystems are fragile, and the relationship between vegetation and permafrost is complicated (Yang et al. 2010); therefore, studies at different spatial scales or different places may reach different conclusions.

There is a lack of long-term observations in dynamics of the alpine ecosystem during permafrost degradation (Yang et al. 2010). Only a few studies have been conducted to quantitatively identify the main drivers to vegetation changes (Zhang et al. 2016). Until now, it remains unclear whether permafrost and seasonally frozen ground degradation make a positive contribution (Cuo et al. 2015; Qin et al. 2016) or negative contribution (Jin et al. 2009; Iijima et al. 2014) to vegetation growth. It is important to reveal the impacts of permafrost degradation on vegetation variations in the Tibetan Plateau.

The objective of this study is to analyze the contributions of driving factors to alpine vegetation variations in different stages of frozen ground degradation. This study simulated the frozen ground distributions by employing the modified Frost Number model and divided the frozen ground into seasonal frozen ground zone, changing zone, and permafrost zone. The simulation period from 1982 to 2012 was divided into five typical stages based on the Normalized Difference Vegetation Index (NDVI) variations. We quantitatively analyzed the relations between precipitation, temperature, the maximum thawing depth, and NDVI in the subzones. We also investigated the contributions of these factors to the annual variation rate of NDVI in different stages and different zones.

Materials and methods

Study area

The Tibetan Plateau, considered to be Earth’s third pole and located in western China (Figure 1a), has an average altitude of above 4,000 m sea level and covers an area of about 2.57 × 10⁶ km² (Zhong et al. 2010). Climate on the Tibetan Plateau is characterized by strong solar radiation, low air temperature, high daily temperature differences, and low annual temperature differences (Zhong et al. 2010). Mean annual air temperatures range from −10 to 10 °C (Sun, Miao, and Duan 2015; Kuang and Jiao 2016) and mean annual precipitation ranges from <50 mm to >1,000 mm from the northwest to the southeast (Ding et al. 2016; Kuang and Jiao 2016). Precipitation on the Tibetan Plateau occurs mainly in the summer. Alpine permafrost is widely distributed in the central region of the plateau (Cheng and Jin 2013). From southeast to northwest of the plateau, the NDVI values decrease gradually (Figure 1b). Forests and shrubs in the southeast are linked to NDVI values mainly >0.7, while meadows and steppes in the middle are linked to NDVI values mainly >0.5 (Figure1b and Figure 2). In the northwest, the low NDVI values are desert (Figure1b and Figure 2). Most of the plateaus are arid and cold, and covered by alpine vegetation. Temperature and moisture conditions at the biological limitation level make the ecosystem highly sensitive to climate or environment changes (Kato et al. 2004; Piao et al. 2011; Gao et al. 2013).

Figure 1. Location of the Tibetan Plateau (a) and distribution of the average annual maximum NDVI for 31 years during 1982 to 2012 (b). Topography of the Chinese mainland is derived from a digital elevation model (http://srtm.csi.cgiar.org).

Figure 2. Vegetation type maps based on the Vegetation Map of China.

Data and processing

Climatic data

Daily average ground surface temperature from 78 standard meteorological stations on the Tibetan Plateau over 1982–2012 was used. Monthly gridded precipitation and air temperature, which are based on the precipitation and air temperature records of 2,472 national meteorological stations and were interpolated using ANUSPLIN developed by the Australian National University (Hutchinson and Xu 2004), were also used in this study. Data assessment shows that each of the gridded box series is highly correlated with the original observational series with a small error (NMIC 2012). Therefore, the gridded values can be used to reflect the variations of precipitation. All temperature and precipitation data are provided by the China Meteorological Administration (CMA) and the National Meteorological Information Center (NMIC) (https://data.cma.cn).

Remote sensing data

Remote sensing has the potential to detect and monitor changes of alpine vegetation at a variety of spatial and temporal scales (Stow et al. 2004). NDVI, which is the most widely used index, is frequently used to assess spatiotemporal changes in vegetation dynamics (Purevdorj et al. 1998; Nemani et al. 2003; Liang et al. 2012). In this study, the latest version of Global Inventory Modeling and Mapping Studies (GIMMS) NDVI3g from 1981 to 2012 (https://ecocast.arc.nasa.gov) was used. The GIMMS NDVI3g has been widely used to study vegetation at global and regional scales because of long time series, wide extent of observation, and short satellite revisit time (Liang et al. 2012; Piao et al. 2015). The dataset was provided at a spatial resolution of 1/12 degree and a 15-day temporal resolution. The product was carefully assembled from different advanced very high resolution radiometer (AVHRR) sensors, removing several satellite errors such as calibration loss, orbit drift, and some special environmental effects like volcanic eruptions (Holben 1986; Pinzon and Tucker 2014). We also obtained the largest annual NDVI value for each pixel in the region using the maximum-value composite method (Holben 1986) to reduce the influence of clouds and atmospheric constituents (Hope et al. 2003; Stow et al. 2004; Bao et al. 2014). The NDVI used are peak growing season values.

The Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type product (MCD12Q1) was used to screen barren or sparsely vegetated pixels (Friedl et al. 2002, 2010). The International Geosphere Biosphere Program (IGBP) classification scheme was used with MCD12Q1 imagery. Areas that were barren or sparsely vegetated from 2001 have not been analyzed (Figure 3).

Figure 3. The difference map between the simulated frozen ground area and the map of permafrost in the Tibetan Plateau, 1:3,000,000 (Li and Cheng 1996). Provided by http://westdc.westgis.ac.cn.

Methods

Temporal and spatial variation of ground surface temperature

Gridded mean annual ground surface temperature (MAGST) is based on daily average ground surface temperature records from 78 national meteorological stations. As most of the meteorological stations are distributed in the eastern part of the Tibetan Plateau, the sparse distribution of the meteorological stations on the west of the Tibetan Plateau has a large effect on the precision of the MAGST data. To improve accuracy, the gridded MAGST is interpolated using linear multiple regression (Wang, Nan, and Zhao 2011). The regression model between MAGST and topographical factors (altitude, latitude, and longitude) is as follows:

(1) $MAGS{T}_{sim}=d+b_1\left[altitude\right]+b_2\left[latitude\right]+b_3\left[longitude\right]$(1)

(2) $\epsilon =MAGS{T}_{sim}-MAGS{T}_{obs}$(2)

(3) $MAGST=MAGS{T}_{sim}+{\epsilon }_{interp}$(3)

where MAGST_sim is the simulation result, MAGST_obs is the observation data from 78 national meteorological stations; b₁, b₂, and b₃ are the multiple regression coefficients; d is the regression constant; $\epsilon$ is the residual value; and $\epsilon$_interp is the grid residual value, which is interpolated by the Kriging method for the entire plateau.

Frost number model

The Frost Number model (Nelson and Outcalt 1987), a statistical-empirical model, is able to simulate and predict frozen ground distributions and has been applied successfully in many regions (Lu, Zhao, and Wu 2017). In this article, the modified Frost Number model was used, a new variable E is introduced to explicitly reflect the thermal effects of frozen and thawed states (Nelson and Outcalt 1987; Nan et al. 2012). Previous studies found that this model can depict the major characteristics of frozen ground distributions over the Tibetan Plateau (Nan et al. 2012; Ran et al. 2015; Lu, Zhao, and Wu 2017).

(4) $F=\frac{\sqrt{DDF}}{\sqrt{DDT}+E\sqrt{DDT}}$(4)

(5) $E=\frac{E_t}{E_f}=\sqrt{\frac{{\lambda }_{\mathrm{t}}{Q}_f}{{\lambda }_f{Q}_t}}$(5)

where F is the surface Frost Number and a threshold value of 0.5 was used to distinguish between permafrost and seasonal frozen ground (Nan et al. 2012; Ran et al. 2015; Lu, Zhao, and Wu 2017), F ≥0.5 and <0.5 are considered as permafrost and seasonal frozen ground, respectively; DDF and DDT are the ground surface freezing and thawing indices and defined by the degree-day totals above or below 0 °C, respectively, expressed in degree-days; E is the calibrated factor; subscript t and f are thaw soil and frozen soil; λ is thermal conductivity (W·m⁻¹·°C⁻¹); Q is the absorbed or released heat (kJ·m⁻³), Q = L·γ·θ, in which L is the mass-based latent heat of fusion or freeze (J·kg⁻¹); γ is the unit weight of soil (kg·m⁻³); θ is the average liquid soil water content of the frozen ground (m³·m⁻³).

The Stefan solution (Stefan 1891; Jumikis 1977; Lunardini 1981) was applied to estimate maximum thawing depth. It is considered to be maximum thawing depth in the permafrost zone, while in other zones are considered to be maximum freezing depth, which is expressed as follows:

(6) $TD=\left\{{ }^{\hspace{0.17em}}\begin{array}{c}\sqrt{\frac{2{\lambda }_{t}DDT}{Q_t}}\\ \sqrt{\frac{2{\lambda }_{f}DDF}{Q_f}}\end{array}\right.$(6)

Contribution of different factors to NDVI variations

We choose precipitation (P), air temperature (T), and maximum thawing depth (D) as the main factors influencing NDVI variations. The similar approach is used to separate the contribution of different factors to variation of evaporation (Rana and Katerji 1998; Yang and Yang 2012; Zhang et al. 2016).

(7) $\begin{array}{rl}& slope=C(P)+C(T)+C(D)+OF\\ & =\frac{\mathrm{\partial }NDVI}{\mathrm{\partial }P}\frac{\mathrm{\partial }P}{\mathrm{\partial }n}+\frac{\mathrm{\partial }NDVI}{\mathrm{\partial }T}\frac{\mathrm{\partial }T}{\mathrm{\partial }n}+\frac{\mathrm{\partial }NDVI}{\mathrm{\partial }D}\frac{\mathrm{\partial }D}{\mathrm{\partial }n}+OF\end{array}$(7)

where slope is the annual variation rate of NDVI. We calculate slope by a simple linear regression (Edwards 1984). C(P), C(T), and C(D) are the contributions of precipitation, air temperature, and thawing depth to the annual variation rate of NDVI, respectively, n is the number of years. OF is the residual contribution. OF may contain human factors and some other uncertain natural factors such as solar radiation and humidity. ∂NDVI/∂P, ∂NDVI/∂T and ∂NDVI/∂D are the slope of the linear regression line between the annual variation rate of NDVI and precipitation, air temperature, and thawing depth, respectively. ∂P/∂n, ∂T/∂n, and ∂D/∂n are the slope of the linear regression line between n and precipitation, air temperature, and thawing depth, respectively.

Mann–Kendall trend analysis

The Mann–Kendall trend test (Mann 1945; Kendall 1948) is one of the widely used non-parametric tests to detect significant trends in time series. To detect the existence of any step change points in a time series x_t = (x₁, x₂, x₃, …x_n), the Mann–Kendall rank statistic d_k is calculated as following:

(8) $d_k=\sum _{i=1}^k{m}_i{\hspace{0.17em}}^{\stackrel{\hspace{0.17em}}{\hspace{0.17em}}}2\le k\le n$(8)

where m_i is the number of later terms in the series whose values exceed x_i .

The mean and variance of the normally distributed statistic d_k are given by

(9) $E\left(d_k\right)=\frac{k\left(k-1\right)}{4}$(9)

(10) $var\left(d_k\right)=\frac{k\left(k-1\right)\left(2k+5\right)}{72}$(10)

The normalized variable statistic UF(d_k) is estimated by the following formula:

(11) $UF\left(d_k\right)=\frac{d_k-E\left(d_k\right)}{\sqrt{var\left(d_k\right)}}\begin{array}{cc}\hspace{0.17em}& \hspace{0.17em}\end{array}k=1,2,3,...,n$(11)

The retrogressive variable statistic UB(d_k) (backward sequence) is also calculated with Eq. (11) but with a reversed series of the data. The intersection of the forward and backward curves statistics UF(d_k) and UB(d_k) indicates the beginning of the step change point (Partal and Kahya 2006; Ye et al. 2013; Wang 2014). Positive UF values indicate an increasing trend, and negative UF values indicate a decreasing trend. Two critical values were set for significance measurement as U_0.05 = ±1.96. The sequence curve of UF(d_k), UB(d_k), and two horizontal lines indicating the critical values are presented in Figure 5.

Figure 4. Simulated subzones in the Tibetan Plateau. Seasonal frozen ground zone, changing zone, and permafrost zone account for about 46.1 percent, 23.3 percent, and 30.6 percent of the study area, respectively.

Figure 5. Mutation times in the annual-maximum and zone-averaged NDVI series from 1982 to 2012 in the three subzones and the entire study area. The intersection of the solid line and the dashed line is the mutation point. (a) The mutation times in permafrost zone. (b) The mutation times in changing zone. (c) The mutation times in seasonal frozen ground zone. (d) The mutation times in the whole Tibetan Plateau.

Results

Permafrost distribution pattern

The seasonal frozen ground and permafrost modeled by the modified Frost Number model are shown in Figure 4. Permafrost is located in the center and the Qilian Mountain in the Tibetan Plateau. Seasonal frozen ground is distributed in the east and south to the permafrost zone. The area of permafrost and seasonal frozen ground is about 3.5 × 10⁵ km² and 7.31 × 10⁵ km², respectively. In order to test the accuracy of the simulation results, the simulated permafrost map was compared with the 1:3,000,000 Map of Permafrost in the Tibetan Plateau (Li and Cheng 1996). The permafrost map was based on field survey and was considered as the standard permafrost distribution map for the Tibetan Plateau. The map has also been widely used in previous permafrost studies (Cheng et al. 2012; Ran et al. 2012; Lu, Zhao, and Wu 2017). We selected the simulated permafrost map from the simulation period (1995–2000) to compare with the permafrost map by Li and Cheng (1996) because this simulation period is close to 1996. As long as a simulation map is consistent with the permafrost map, we consider that the simulated results are consistent with the field survey data. We calculated the area of the simulated map and verified map, and the simulated map is 90.5 percent similar to the verified map.

According to the simulated results, we classified the permafrost map into three subzones: seasonal frozen ground zone, changing zone, and permafrost zone. Seasonal frozen ground zone, changing zone, and permafrost zone account for 46.1 percent, 23.3 percent, and 30.6 percent of the study area, respectively (Figure 4). Here, we defined the changing zone as the region transforming from permafrost into seasonal frozen ground over the past 31 years.

NDVI mutation time

The abrupt changes of NDVI time series from 1982 to 2012 is identified by the Mann–Kendall method. The abrupt changes occurred at similar times across the whole Tibetan Plateau and the three subzones (Figure 5). The intersection of the dash line and the solid line shown in Figure 5 is the mutation point of NDVI. The mutation times of NDVI series in each region have been counted in Table 1. Abrupt changes affecting the whole Tibetan Plateau occurred at four times: 1984–1985, 1987, 2001–2003, and 2008–2009.

Table 1. Mutation times from 1982 to 2012.

Category	Mutation times
Permafrost zone	1987, 2001–2003, 2008–2009
Changing zone	1984, 2001–2003, 2008–2009
Seasonal frozen ground zone	1984, 2001–2003, 2008–2009
The whole TP	1985, 1987, 2001–2003, 2008–2009

Based on the four mutation times, the three zones can be divided into five stages during the past 31 years: 1982–1987, 1988–1994, 1995–2000, 2001–2007, and 2008–2012. The five stages were named stage 1 to stage 5, respectively.

Spatiotemporal pattern of permafrost

The spatiotemporal distributions of permafrost and seasonal frozen ground in the five stages are shown in Figure 6. The total area of permafrost in each stage is 6.727 × 10⁵ km², 5.548 × 10⁵ km², 5.288 × 10⁵ km², 3.554 × 10⁵ km², and 3.551 × 10⁵ km², respectively (Figure 6f). The permafrost area is decreasing over the five stages. The decreased areas between two stages are 11.79 × 10⁴ km², 2.60 × 10⁴ km², 17.34 × 10⁴ km², and 0.03 × 10⁴ km², respectively. The area of permafrost degradation is the largest from stage 3 to stage 4. The decreased area is the smallest from stage 4 to stage 5, due to the slow warming of the climate. The degradation of permafrost mainly begins in the east and southeast of the Tibetan Plateau. These parts also correspond to the high-temperature areas. Over the past 31 years, the seasonal frozen ground has been embedded into the center of the continuous permafrost zone and the large area of continuous permafrost has been degraded into pieces (Figure 6a–e). Permafrost is no longer a continuous aquiclude for groundwater.

Figure 6. Distribution and area changes of permafrost and seasonal frozen ground during the five stages. (a)–(e) Distributions of permafrost and seasonal frozen ground in stage 1, stage 2, stage 3, stage 4, and stage 5, respectively. (f) Area changes of permafrost and seasonal frozen ground over the five stages.

Figure 7. The relationships between NDVI and precipitation, temperature, and maximum thawing depth in the three subzones. The black spots represent 31-year-averaged data in each pixel. The red dots represent the averaged data every 0.01 NDVI interval.

Figure 8. Contributions of precipitation, temperature, and maximum thawing depth to NDVI over the five stages.

Correlations of the main factors with NDVI

In previous studies, temperature, precipitation, and groundwater were considered as the main factors determining the growth of alpine vegetation in the cold and arid Tibetan Plateau (Guisan, Theurillat, and Kienast 1998; Yang et al. 2010; Zhong et al. 2010; Manel et al. 2012). To reveal the impact of temperature and moisture conditions on alpine vegetation, we calculated the relation of precipitation, temperature, and maximum thawing depth to NDVI variation in the three subzones (Figure 7). There is a good positive correlation between precipitation and NDVI in all the three subzones (Figure 7a–c). The slopes of the relationships between NDVI and temperature in the seasonal frozen ground zone and permafrost zone are positive (Figure 7d and e), whereas the slope in the changing zone is negative when NDVI values are >0.5 and positive when NDVI is <0.5 (Figure 7f). The slopes of the relationships between NDVI and the maximum thawing depth in the permafrost zone were are (Figure 7g). The relationships between NDVI and the maximum thawing depth in the changing zone and seasonal frozen ground zone are different. In the changing zone where NDVI is <0.3 (Figure 7h), and in the seasonal frozen ground zone where NDVI is between 0.44 and 0.7 (Figure 7i), the slopes of the relationships between the maximum thawing depth and NDVI are negative.

The correlation coefficients (R²) can be used to determine the relevance of the three factors to NDVI. Based on the correlation coefficients, the correlation between NDVI and precipitation is the most significant in all the subzones (Figure 7a–c). In each subzone, the maximum thawing depth has more significant relation with vegetation growth than the temperature.

Contribution of factors to NDVI

For further understanding the contribution of each driving factor to vegetation growth, we calculate the contribution rate of precipitation, temperature, and maximum thawing depth at each stage based on Equation (8) using annual data of the three factors (annual precipitation, annual average temperature, and annual maximum thawing depth). The contribution rates of the three factors in each pixel are shown in Figure 8.

The contributions of precipitation to NDVI variations are mainly positive over the five stages except for the 4th stage, as shown in row 1 of Figure 8. The area with positive correlation was increased from 52 percent to 62 percent (Table 2). At the 4th stage, the regions with negative contributions of precipitation to NDVI are mainly located in the Three Rivers Source Region (Yangtze River, Yellow River, and Mekong River) and the region around Lhasa (Figure 8). The positive contributions of precipitation mainly occurred in the lower-NDVI areas in the central and north of the Tibetan Plateau, while the negative ones are in the higher-NDVI areas in the east and northeast of Tibetan Plateau. In summary, the contributions of precipitation to the vegetation in the cold and arid Tibetan Plateau are different in space and time.

Table 2. The area percentage of factors (precipitation, temperature and maximum thawing depth) driving NDVI change over the five stages.

Driving factor	Correlation	1st stage	2nd stage	3rd stage	4th stage	5th stage
Precipitation	Negative	46%	48%	38%	54%	41%
	Positive	54%	52%	62%	46%	59%
Temperature	Negative	48%	63%	56%	60%	68%
	Positive	52%	37%	44%	40%	32%
Maximum thawing depth	Negative	50%	16%	29%	49%	53%
	Positive	50%	84%	71%	51%	47%

Generally, temperature had a negative contribution to alpine vegetation over the five stages except for the 1st stage, which is shown in row 2 of Figure 8. The area with negative contribution rate accounted for 56 percent to 68 percent of the study area and was mainly in the west, center to southeast (Table 2). At the 1st stage, about 48 percent of the area had positive contributions to the annual variation rate of NDVI.

The effects of maximum thawing depth on the annual variation rate of NDVI also varied with time and space. From the 2nd to 4th stages, the area with positive contributions is more than half (Table 2). The areas with higher NDVI, such as in the east, the maximum thawing depth had a negative contribution to the annual variation rate of NDVI. While the maximum thawing depth had a positive role in lower-NDVI vegetation in the west, which are shown in the averaged contribution (Figure 8). This typical spatial distribution pattern can also be found at the 3rd and the 5th stages, as shown in row 3 of Figure 7.

The zone-averaged contribution rates in the three subzones over the five stages are shown in Figure 9. The contributions of precipitation and maximum thawing depth are positive in in the three subzones at the 3rd and 5th stages (Figure 9) where the slope of NDVI is positive (Figure 10). During these two stages, in the region of changing zone and seasonal frozen ground zone, the maximum thawing depth had a positive contribution to the annual variation rate of NDVI. However, with the climate warming, the positive contribution of thaw depth became smaller (Figure 9b and c). The contribution of maximum thawing depth was even changed from positive to negative in the permafrost zone from the 3rd stage to the 5th stage (Figure 9a). Within the two stages, in the regions of changing zone and seasonal frozen ground zone, temperature had a negative contribution to the variation rate of NDVI (Figure 9b and c). As the climate warming, the negative contribution of temperature increased. In the permafrost zone, the contribution of temperature on vegetation growth is similar to that of maximum thawing depth (Figure 9a).

Figure 9. The zone-averaged contribution of the three driving factors at the five stages. (a)–(c) Contributions in the permafrost zone, changing zone, and seasonal frozen ground zone, respectively.

Figure 10. The zone-averaged NDVI trend in the three subzones at the five stages.

The NDVI trend can help us understand vegetation cover change over time. Figure 10 shows that NDVI in the three subzones had similar trends at each stage. In the 3rd and 5th stages, NDVI increased strongly in all three subzones. It decreased slightly in the 2nd and 4th stages and stayed about the same in stage 1.

Discussion

Permafrost map accuracy valuation

In fact, it is difficult to obtain enough field observation data to outline the distribution of permafrost in the Tibetan Plateau. The physically based modeling for permafrost is also difficult because of the difficulty in acquiring accurate model parameters for the whole region. On the basis of the available data, we selected the Frost Number model to carry out large-area permafrost simulations in the Tibetan Plateau. The main data used in the model is ground surface temperature. In addition to the effects of ground surface temperature, permafrost changes are also affected by vegetation cover, soil and bedrock thermal conductivity, moisture content, and other factors. When these physical parameters are considered, the accuracy of permafrost simulation can be improved (Brown 1966; Guo, Wang, and Wang 2017).Compared with previous studies, the simulated areas of permafrost may be different at different stages, but the general trend of permafrost degradation is in consistent with previous studies (Li and Cheng 1999; Nan et al. 2012; Guo, Wang, and Wang 2017). The degradation of permafrost is characterized by the reduced area, as well as the increased surface burial depth (increasing the active layer thickness). We compared our simulation results with the 1:3,000,000 permafrost map (Li and Cheng 1996). The similarity of 90.5 percent strongly suggests that this statistical-empirical model can be successfully used to simulate the distribution of permafrost. However, there is still 9.5 percent area where there is difference between the simulation result and the survey map of permafrost distribution in the Tibetan Plateau. These different regions are sporadically distributed around the simulated permafrost, and mainly concentrated in the east region, especially near the source region of the Yellow River (Figure 3). The reason for this inconsistency is likely the inaccuracy of the survey map itself, because most of the permafrost boundary lines were empirically drawn on the topographic map by hand. Although the sparse borehole observations and the easy-access terrain survey can obtain accurate data, these survey data accumulated from nearly two decades till 1996, so the timeliness maybe another factor that cause the two maps different. The last reason may be the scale, the scale of the survey map is 1:3,000,000, while the modeling resolution is 1,000 m. Despite these differences, the simulated map is generally in good agreement with the survey map.

Frozen ground classification

The thawed permafrost surface would result in lowering of the groundwater table above the permafrost surface. In an extreme case, when the permafrost disappeared, it would affect recharging of the groundwater over the permafrost surface into deep aquifers. These changes will eventually affect the alpine vegetation that depends on the groundwater over the permafrost (Cheng and Wu 2007; Cheng and Jin 2013). Therefore, we separate permafrost degradation zones from the frozen ground map. Here, we defined the changing zone as the region transforming from permafrost into seasonal frozen ground over the past 31 years. In this way, the frozen ground was divided into seasonal frozen ground zone, changing zone, and permafrost zone. Seasonal frozen ground zone, changing zone, and permafrost zone account for 46.1 percent, 23.3 percent, and 30.6 percent of the study area, respectively (Figure 4). Over the past 31 years, permafrost areas have been decreasing and the large area of continuous permafrost has been degraded into pieces. The changing zone is mainly distributed between the permafrost zone and the seasonal frozen ground zone, especially around the river basins and large lakes.

Factors correlation to NDVI variation

To some extent, permafrost and its overlying alpine vegetation are interdependent. Changes of alpine vegetation can reflect variations of permafrost indirectly. In some typical areas, vegetation change is considered as an indicator of the existence and variation of permafrost (Frauenfelder et al. 1998; Gruber and Hoelzle 2001). To distinguish the impacts of climate and environment on vegetation, we calculated the correlations of precipitation, temperature, and maximum thawing depth with NDVI in these three typical regions. We also divided the past 31 years into five typical periods according to the changing trend of NDVI. These data provided a new basis to determine the relevance of these factors on alpine vegetation and their changes in space and time. During the past 31 years, NDVI has undergone four mutations and consequently had five different stages. Correlations with NDVI, ranked in turn, is precipitation, maximum thawing depth, and temperature.

Precipitation was the main factor driving NDVI variation and had a positive relationship with NDVI in the three subzones (Figure 7a–c). This is consistent with many previous studies (Song, Zhou, and Ouyang 2005; Shen et al. 2011). However, we used the average data to draw this conclusion. In fact, different regions and times would come to different conclusions. We must be careful when we analyze the relationship between precipitation and NDVI. Therefore, we also calculated the contributions of precipitation in different regions and different stages to NDVI changes. Results show that the contributions varied with time and space. We will discuss this in more detail in the next section.

The maximum thawing depth is the second factor and has a positive relationship with NDVI, except for in changing zone when NDVI is <0.3 and in seasonal frozen ground zone when NDVI values are between 0.44 and 0.7 (Figure 7g–i). The affected area where the maximum thawing depth has a positive relationship with NDVI exceed a half in the Tibetan Plateau, and gradually decreases from the second stage (Figure 8, row 3). In general, based on the averaged contribution rate in the three subzones, the maximum thawing depth had a positive relationship with NDVI during the vegetation increasing periods (the 3rd and 5th stages). In the end, the effect of active layer deepening is similar to the effects of temperature rise. When the active layer is very shallow, certain plant species or vegetation types may expand when the active layer deepens. Further deepening of the active layer leads to drying of the soil at rooting depth, thus preventing water and nutrient uptake. The thawed layer can provide available liquid water, nutrients, or minerals for alpine vegetation, which cannot be utilized by vegetation when frozen. Only in the thawed layer can vegetation grow and develop. From this point of view, the maximum thawing depth had a positive relationship with vegetation (Cuo et al. 2015; Qin et al. 2016). Permafrost provides water source for arid vegetation by controlling the water table. When the large area of continuous permafrost is degraded into pieces—that is, the maximum thawing depth is increasing—the thawing of permafrost results in lowering of the water table (Cheng and Wu 2007). The alpine vegetation that depends on the water table would degenerate. In this case, the maximum thawing depth had a negative impact on the vegetation growth (Jin et al. 2009). This negative effect, such as in the changing zone with NDVI <0.3, can be interpreted as the dependence of drought-tolerant vegetation on groundwater. Permafrost will degenerate as global warming, and vegetation will continue to be endangered by groundwater scarcity, even being replaced by drought-tolerant vegetation.

In cold areas, accumulated temperature directly affects vegetation coverage (Shen et al. 2011; Sun et al. 2013). In areas with higher accumulated temperature, the vegetation coverage is generally higher, and vice versa. In areas with high vegetation cover (NDVI >0.5), accumulated temperature can meet the needs of local vegetation. Our study shows that a further rise in temperature does not lead to a further increase of NDVI. The influence of rising temperatures on vegetation composition and cover depends on the temperature optima for specific species, plant communities, and biomes and on how close the temperature is to the limits of these optima. In areas with lower vegetation cover, such as the northwest of the plateau or NDVI <0.5, mainly steppe, it is cold and accumulated temperature is still an important demand for the growth of vegetation. In areas with higher vegetation cover, such as the southeast of the plateau or NDVI >0.7, mainly forest or shrub, accumulated temperature is not a necessary demand for the growth of vegetation. It can be inferred that the effects of temperature on alpine vegetation may depend on the fractional vegetation cover and species. If the fractional vegetation cover is lower, the rising temperature may promote the vegetation growth; otherwise, the rising temperature may against the vegetation growth.

Driving factors of NDVI change

We calculated the contribution rates of the three main factors (precipitation, temperature, and maximum thawing depth) pixel by pixel over the five stages. In general, in the 1st, 2nd, and 4th stages, the NDVI reduced slightly (Figure 10). During these stages, the rising temperature was not conducive to the growth of alpine vegetation. Temperature may result in vegetation degradation, and this impact was accelerated with time (Figure 9). The maximum thawing depth played a minor role in preventing vegetation degradation (Figure 9), and the increasing maximum thawing depth is not conducive to the vegetation growth in these stages. In these three stages of vegetation degradation, the effect of precipitation on vegetation varies in different regions. The contribution of precipitation would change and even have negative effects on vegetation in different regions at different stages. When the precipitation amount is sufficient, the alpine vegetation change may not be related to the precipitation. In arid areas, the negative correlation may be due to lagging of vegetation change relative to the precipitation. The impact of precipitation on the vegetation in cold and arid regions is different spatially and temporally.

During both the vegetation increasing stages (the 3rd and 5th stages) and the vegetation degrading stages (the 1st, 2nd, and 4th stages), temperature is not conducive to vegetation growth (Figure 9a–c). Temperature hinders the growth of vegetation or accelerates vegetation degradation. The enhanced negative effect over stages may be due to climate warming. The effect of temperature rise on alpine vegetation has two sides. The increase of temperature can meet the accumulated-temperature demand of the alpine vegetation in cold area. On the contrary, the high temperature also can lead to the increase of evaporation, result in the loss of water resources, and ultimately limit the vegetation growth in the arid region. Over the past few decades, there have been differences in the temperature rise in different areas of the Tibetan Plateau (Guo and Wang 2012; Gao et al. 2015; Kuang and Jiao 2016; You, Min, and Kang 2016). This understanding can explain why the effects of temperature on vegetation varied at different stages and different places. The finding is consistent with previous studies (e.g., Liang et al. 2012).

There are still some driving factors (such as wind, solar radiation) on NDVI dynamics have not been considered in this study. The effects of these driving forces on NDVI dynamics may be the topic for further studies. Different human activities are also important factors that influencing NDVI changes, especially in the National Nature Reserve of Three Rivers Source.

In this study, we used the GIMMS-NDVI data, with 8 km spatial and 15-day temporal resolution (Holben 1986; Pinzon and Tucker 2014). It is suitable for large-scale research, but the coarse resolution cannot reflect the local changes. In the Tibetan Plateau, vegetation growth season is very short, and only data with higher temporal resolution can accurately reflect the characteristics of rapid seasonal vegetation changes in a year. NDVI also cannot reveal the vegetation replacement that has occurred in the Tibetan Plateau over the past few decades (Klanderud and Birks 2003; Song, Zhou, and Ouyang 2005; Ma et al. 2015; Seddon et al. 2016). NDVI does not characterize sparse vegetation well (Elmore et al. 2000). Future studies on vegetation with small NDVI should adopt other indexes (e.g., vegetation coverage) to correctly reflect the real variation of sparse vegetation.

We only analyzed the annual average data (precipitation, temperature, maximum thawing depth, and NDVI data), and the results may contain some smoothed effects resulted from average. However, based on multi-annual averages, precipitation has a good positive correlation with vegetation in most cases, though its negative effects can also be found in some stages and some spaces. If we increase the temporal resolution from 1 year to semi-month (consistent with the GIMMS resolution), it may be possible to find more precise changes in the impact of the factors on the alpine vegetation.

Conclusions

In this study, we simulated the frozen ground distributions by employing a modified Frost Number model and divided the frozen ground into seasonal frozen ground zone, changing zone, and permafrost zone in the Tibetan Plateau. The area of permafrost demonstrated a fluctuation but presented a generally decreasing trend from 1982 to 2012. The decreased areas is 31.76 × 10⁴ km², and the degradation rate is about 1.02 × 10⁴ km²/yr. The greatest rate of degradation occurred around 2000 due to the accelerated climate warming in the Tibetan Plateau.

We quantitatively analyzed the relations between precipitation, temperature, the maximum thawing depth and NDVI in the subzones, and contributions of the three factors to the annual variation rate of NDVI in different stages and different zones from 1982 to 2012. During the entire study period, precipitation has the largest average contribution to the variation rate of NDVI, followed by the maximum thawing depth and temperature.

Acknowledgments

The authors thank the two reviewers and the editor for the insightful and constructive comments.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Strategic Priority Research Program of Chinese Academy of Sciences [XDA20100103] the National Natural Science Foundation of China (NSFC 91747204, 41330634, and 41072191), the Fundamental Research Funds for the Central Universities, and the Guangdong Provincial Key Laboratory of Soil and Groundwater Pollution Control (No. 2017B030301012).