Time Series Analysis of Land Cover Change in Dry Mountains: Insights from the Tajik Pamirs

Greening and browning trends in vegetation have been observed in many regions of the world in recent decades. However, few studies focused on dry mountains. Here, we analyze trends of land cover change in the Western Pamirs, Tajikistan. We aim to gain a deeper understanding of these changes and thus improve remote sensing studies in dry mountainous areas. The study area is characterized by a complex set of attributes, making it a prime example for this purpose. We used generalized additive mixed models for the trend estimation of a 32-year Landsat time series (1988–2020) of the modified soil adjusted vegetation index, vegetation data, and environmental and socio-demographic data. With this approach, we were able to cope with the typical challenges that occur in the remote sensing analysis of dry and mountainous areas, including background noise and irregular data. We found that greening and browning trends coexist and that they vary according to the land cover class, topography, and geographical distribution. Greening was detected predominantly in agricultural and forestry areas, indicating direct anthropogenic drivers of change. At other sites, greening corresponds well with increasing temperature. Browning was frequently linked to disastrous events, which are promoted by increasing temperatures.

Globally, anthropogenic factors, such as fertilization and irrigation, are the drivers with the highest impact on greening and by far outweigh the effect of climate change [4,13,29]. Yet, on a regional level, there are strong differences. Drivers of greening caused by climate change play a major role in the high latitudes, in high mountains, and in drylands [4,11,[25][26][27][30][31][32][33]. A prolonged growing season caused by increased temperatures is responsible for most greening in the boreal zone and the arctic tundra [10,13,[15][16][17][18][19][20]29,34], whereas increased precipitation is the major driver for greening in arid and semi-arid areas [6,11].
numbers. Thereby, the presented study provides a novel example for a time series analysis of irregular data in a dry and heterogeneous environment.

Study Area
The study area covers 1280 km 2 of the western part of the Rushan Range in the Western Pamirs of Tajikistan ( Figure 2). Administratively, it belongs to the Rushan and Shugnan districts of the Gorno-Badakhshan Autonomous Oblast (i.e., Province; GBAO). It is delimited by the river Gunt in the south, the river Panj in the west, the river Bartang in the north, and the 72°E meridian in the east. The highest peaks reach more than 5000 m a.s.l., while the minimum elevation is found in the northwestern corner at 1980 m a.s.l. The harsh climate is characterized by severe, long winters with minima below −35 °C, followed by short summers. The annual mean temperatures range between −7.1 °C on the Fedchenko Glacier and 8.7 °C at the meteorological station in the region's center in Khorog [51,[92][93][94]. Breckle [94] reports an annual precipitation in the Western Pamirs of between 90 and 400 mm, with a distinct winter maximum. This region is also part of the upper watershed of the Amu-Darya, rendering it crucial for the water supply for millions of people living downstream [95]. The Western Pamirs are considered to be a biodiversity hotspot [48,49,51,96,97], and large parts of it are (though insufficiently) protected by the Tajik National Park [98][99][100]. Around 2000 vascular plant species, including 160 endemics and a high number of threatened animal species, such as the snow leopard (Panthera uncia) and the Marco Polo sheep (Ovis ammon polii), are native to this area [49,101]. Furthermore, numerous wild progenitors of cultured fruit trees are found in the region, constituting a valuable genetic resource [50,93,96]. The Western Pamirs are also the habitat of rare juniper forest patches and Tugai forests (azonal forests in river valleys), which play important roles in the water regulation of rivers. Moreover, a large number of endemic

Study Area
The study area covers 1280 km 2 of the western part of the Rushan Range in the Western Pamirs of Tajikistan ( Figure 2). Administratively, it belongs to the Rushan and Shugnan districts of the Gorno-Badakhshan Autonomous Oblast (i.e., Province; GBAO). It is delimited by the river Gunt in the south, the river Panj in the west, the river Bartang in the north, and the 72 • E meridian in the east. The highest peaks reach more than 5000 m a.s.l., while the minimum elevation is found in the northwestern corner at 1980 m a.s.l. The harsh climate is characterized by severe, long winters with minima below −35 • C, followed by short summers. The annual mean temperatures range between −7.1 • C on the Fedchenko Glacier and 8.7 • C at the meteorological station in the region's center in Khorog [51,[92][93][94]. Breckle [94] reports an annual precipitation in the Western Pamirs of between 90 and 400 mm, with a distinct winter maximum. This region is also part of the upper watershed of the Amu-Darya, rendering it crucial for the water supply for millions of people living downstream [95]. The Western Pamirs are considered to be a biodiversity hotspot [48,49,51,96,97], and large parts of it are (though insufficiently) protected by the Tajik National Park [98][99][100]. Around 2000 vascular plant species, including 160 endemics and a high number of threatened animal species, such as the snow leopard (Panthera uncia) and the Marco Polo sheep (Ovis ammon polii), are native to this area [49,101]. Furthermore, numerous wild progenitors of cultured fruit trees are found in the region, constituting a valuable genetic resource [50,93,96]. The Western Pamirs are also the habitat of rare juniper forest patches and Tugai forests (azonal forests in river valleys), which play important roles in the water regulation of rivers. Moreover, a large number of endemic dwarf shrub species dominate the upper vegetation belts. Further information on the vegetation in the study area is provided in Section 2.2.5.
We collected Landsat images of the study area and built a raster time series of the modified soil adjusted vegetation index (MSAVI; 2.2.1). Additionally, we gathered time series of temperature and precipitation (2.2.2), population (2.2.3), and livestock data (2.2.4). Furthermore, we used the vegetation and land cover data of a historic Soviet vegetation map and ground-truthing records from fieldwork (2.2.5), as well as topographic variables derived from the shuttle radar topography mission (SRTM) digital elevation model (DEM; 2.2.6).

Landsat Data
We used the Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) on demand interface (https://espa.cr.usgs.gov/, accessed on 6 July 2020) provided by the United States Geological Survey (USGS) to obtain preprocessed satellite scenes and derived the MSAVI [71] from the Level-2 Surface Reflectance [103] of the sensors Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI with a spatial resolution of 30 m. The calculation formula of the MSAVI is as follows: where ρ is the surface reflectance of the near-infrared band and ρ is the surface reflectance of the red band.
Furthermore, a cloud and cloud shadow mask was calculated according to the CFMask algorithm [104][105][106]. Altogether, we obtained 1267 scenes of the study area between 1988 and 2020, which we combined in a 32-year raster image time series.

Data
We collected Landsat images of the study area and built a raster time series of the modified soil adjusted vegetation index (MSAVI; 2.2.1). Additionally, we gathered time series of temperature and precipitation (2.2.2), population (2.2.3), and livestock data (2.2.4). Furthermore, we used the vegetation and land cover data of a historic Soviet vegetation map and ground-truthing records from fieldwork (2.2.5), as well as topographic variables derived from the shuttle radar topography mission (SRTM) digital elevation model (DEM; 2.2.6).

Landsat Data
We used the Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) on demand interface (https://espa.cr.usgs.gov/, accessed on 6 July 2020) provided by the United States Geological Survey (USGS) to obtain pre-processed satellite scenes and derived the MSAVI [71] from the Level-2 Surface Reflectance [103] of the sensors Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI with a spatial resolution of 30 m. The calculation formula of the MSAVI is as follows: where ρ N IR is the surface reflectance of the near-infrared band and ρ red is the surface reflectance of the red band. Furthermore, a cloud and cloud shadow mask was calculated according to the CFMask algorithm [104][105][106]. Altogether, we obtained 1267 scenes of the study area between 1988 and 2020, which we combined in a 32-year raster image time series.

Climate Data
We used the 2-m air temperature and precipitation data from three different data sources. Daily climate data were obtained from the meteorological station in Khorog, located in the Panj valley at 2084 m a.s.l. The data were downloaded from the Climate Data Online (CDO) archive of the National Centers for Environmental Information [107]. Furthermore, we utilized reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis v5-Land (ERA5-Land) dataset with a spatial resolution of 0.1 • [108], and the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2) dataset with a spatial resolution of 0.5 • × 0.625 • [109], as these products proved to be the most suitable for vegetation analysis in mountain regions with limited data availability [110].

Population Data
Since the demographic development of a given region is a critical factor for both its use of resources and its potential to originate labor force for cultivation, we included the time series of population numbers in GBAO in our analysis. We combined two sources published by the State Agency for Statistics under the President of the Republic of Tajikistan (TajStat, [111,112]) in a 28-year time series (1990-2018) for the entire GBAO and in a 25-year time series  for the individual districts.

Livestock Data
For the discussion of grazing impact on the development of land cover and vegetation, we used numbers of large livestock (i.e., mainly cattle) collected by TajStat [113][114][115], which differentiate between the entire GBAO and each district. A consistency check of the data revealed that the time series could not be completely combined, as the numbers in the oldest dataset (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) are three times lower than the numbers in the two other datasets (1992-2011). Therefore, we decided to aggregate only the latest two of the datasets and analyze the resulting two time series separately from each other.

Vegetation and Land Cover Data
To evaluate the results of the MSAVI trend analyses, we used a historical vegetation map produced by the Soviet botanist Okmir Agakhanjanz in 1958-1961 (see [97] for details), which distinguishes between 12 different vegetation and land cover classes in the study area (Table 1). Additionally, to check the information from the maps and the results of the trend analyses for plausibility, we recorded 76 ground-truth field plots all over the study area.

Topographic Data
We evaluated relationships between the trends and topography using six topographic variables derived from the SRTM DEM with a resolution of 1 arc-second (30 m) [102]. Elevation was directly extracted from the DEM and used as the basis to calculate the slope, aspect, ruggedness, and distance to the valley bottom using the function terrain of the R package raster [116]. The aspect was transformed into north and east exposedness by calculating the cosine and sine of aspect, respectively [117].

Time Series Analysis Using Generalized Additive Mixed Models
Generalized additive mixed models (GAMMs) [86,118] were applied as an innovative approach that has rarely been used in remote sensing studies, using the function 'gamm' of the R package 'mgcv' [119]. Depending on the type of time series used, our response variables were MSAVI, temperature, precipitation, population, and livestock numbers, respectively. Their expected values were estimated by the explanatory variable 'year'. We additionally included the explanatory variables 'day of the year' for the MSAVI time series and the variable 'month' for the temperature and precipitation time series. In doing so, we were able to separate the seasonal from the inter-annual component and hence could uncover gradual as well as abrupt land cover changes. In contrast to seasonal changes, which are mostly dependent on annual temperature and rainfall interrelations altering plant phenology, gradual changes can be caused by inter-annual climate variability or land degradation and land management. Abrupt changes are usually caused by disturbances and extreme events. Due to the different periods and spatial scales at which the different response variables operate, the trend comparisons had to be limited to qualitative discussions.
GAMMs are extensions of generalized additive models (GAMs; [120]) and, as such, are able to estimate non-linear trends. They are a preferred method for ecological studies since they allow for complex shapes of response and are well-suited to cope with irregular spacing in a time series. Similar to other regression models, in a GAMM, a univariate response variable is related to one or multiple explanatory variables. As in generalized linear models (GLMs; [121]), the mean of the response variable in GAMMs is related to the explanatory variables by a link function. However, as opposed to parametric models, GAMMs make use of smoothing functions of the explanatory variables. The resulting smoothing curves can take on many forms and allow for non-linear relationships between the response and explanatory variables [66]. The general equation of a GAM is: where g is a specified link function, (µ) is the expected mean response, and ∑f(X) is the smooth functions of the covariates (X). In GAMMs, a random factor is added to this general GAM equation. For a detailed description, readers are referred to Wood [86].
In this study, we used a Gaussian distribution and two types of smoothers-i.e., (1) thin plate regression splines for the development of the data values over the years and (2) cyclic cubic splines to account for the seasonality in the MSAVI, temperature, and precipitation data. The value of the smoothness parameter had to balance between a model that was too simple to sufficiently represent the data and an exceedingly complex model that overfitted the data. With this in mind, the degree of smoothness was optimized by treating it as a random effect using the maximum likelihood method [65,86,122]. As the values in our time series are mostly irregularly distributed in time, we introduced a continuous-time first-order autoregressive process (CAR(1)) to account for autocorrelation in the data that Remote Sens. 2021, 13, 3951 8 of 25 may bias the trend results [65,118]. To compensate for differences caused by the sensors, which may otherwise distort the trend results, we used 'sensor type' (a factor variable with the three levels Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI) as a random effect (random intercept and random slope) variable in the GAMM models for the analyses of the MSAVI time series [86].
We applied the GAMM function to each time series, and this approach allowed us to address the time, amount, and direction of change within the trend component. For the MSAVI time series, we present significance and amount of change as a map ( Figure 3). In doing so, the detected changes were localized on the vegetation map and subsequently further evaluated. In particular, for the 1200-pixel time series, we conducted further statistical analyses to obtain deeper insights into the development of the land cover. These 1200 comprise 100 time series for each of the 12 land cover types. Seventy-six of these locations were visited in the field, giving us a better understanding of the on-site situation. The remaining 1124 locations were stratified randomly and selected using the function 'stratified' of the R package 'fifer' [123].  This visual evaluation was supported by a quantitative analysis, which revealed that trend types of the land cover classes, to some extent, differ fundamentally from the situation in the entire study area and from each other ( Figure 4). In four classes (mountain Tugai, juniper vegetation, mountain meadows, and cultivated land), the number of time series with significant greening exceeds the number of time series with no trend. The share of time series with significant browning is particularly high in mountain Tugai and cultivated land. Time series with no significant trend dominate, especially in mountain steppes and cryophytic and subnival vegetation. Nival areas also show a large proportion of browning and no significant change. However, the very low number of valid observations in these time series suggests spurious trends. Therefore, we decided to refrain from further interpretation.
In the next step, we deepened the analysis based on the characteristics of 1100 test time series, with 100 from each land cover class (except nival areas). First, we examined whether the single years showed significant greening, significant browning, or no significant change. We ascertained that 55% have no significant change in each of the 32 analyzed years. In 36.7%, MSAVI increased significantly in all years-i.e., we detected a continuous linear trend over the entire study period. A continuous browning trend was revealed for less than 1% of the time series. Non-linear shapes-i.e., a mix of greening, Furthermore, we identified points of significant change in the time series. In contrast to linear models, where each point in the time series always has the same slope of the trend, in non-linear models, such as GAMMs, each time point may have a different slope. To solve this difficulty, we estimated the first derivative of the fitted splines by calculating finite differences for numerous time points in the time series [65]. Then, we evaluated if the first derivatives at these time points comply with the null hypothesis of no change. This was achieved by calculating a 95% confidence interval for the derivative estimates. Points in the time series where zero is outside the confidence interval of the first derivative indicate times of significant change. For a detailed description of this method, readers are referred to Simpson [65].

Vegetation Greening and Browning Trends as Indicated by MSAVI Time Series
In this study, we refer to a significant positive (i.e., increasing) trend of MSAVI values over the study period (1988-2020, p < 0.05) as greening, whereas a significant negative (i.e., decreasing) trend (1988-2020, p < 0.05) is called browning. In total, we analyzed more Remote Sens. 2021, 13, 3951 9 of 25 than 1.4 million pixel time series of an area of 1280.2 km 2 . For 16.5% of these time series, there were too few observations to perform statistical analyses. In the remaining time series, we detected greening for 385.9 km 2 (36.1%). In contrast, only 17.2 km 2 (1.6%) showed browning. The majority of the time series (62.3%, 666.5 km 2 ) did not exhibit a significant trend (p > 0.05) of MSAVI increase or decrease (Figure 4).   Then, we compared the amount of change-i.e., the difference in MSAVI of the component smooth functions for the variable 'year' between 1988 and 2020 and between the different types of times series based on their median values: for all non-NA time series, the median is 0.01, for time series with a significant greening trend it is 0.02, and for time series with a significant browning trend it is −0.03 (Table 3). To put these rather small values into perspective, we need to consider that the median MSAVI calculated over all 1100 test plots for the whole time series is only 0.06. The comparison of the absolute values did not unveil patterns of change between the land cover classes. In such cases, a common approach is to compare the classes in terms of some measure of location [124,125]. We A visual evaluation of the trends suggests that the time series with no significant changes are distributed all over the study area and, in particular, on slopes, ridges, and high plateaus. Greening could be detected predominantly in the western part of the study area-i.e., in Panj and its adjacent valleys, along the valleys of Bartang in the north, its tributary valleys Geisev in the north-central, and Ravmed in the northeast-and in Gunt Valley in the south of the study area. Time series with significant browning are partly clustered in valley bottoms, with the biggest area located in Red Valley in the northwest of the study area. The parts of the study area with too few observations for meaningful statistical analyses predominantly belong to the nival areas with a permanent snow or ice cover across the entire year ( Figure 3).
This visual evaluation was supported by a quantitative analysis, which revealed that trend types of the land cover classes, to some extent, differ fundamentally from the situation in the entire study area and from each other (Figure 4). In four classes (mountain Tugai, juniper vegetation, mountain meadows, and cultivated land), the number of time series with significant greening exceeds the number of time series with no trend. The share of time series with significant browning is particularly high in mountain Tugai and cultivated land. Time series with no significant trend dominate, especially in mountain steppes and cryophytic and subnival vegetation. Nival areas also show a large proportion of browning and no significant change. However, the very low number of valid observations in these time series suggests spurious trends. Therefore, we decided to refrain from further interpretation.
In the next step, we deepened the analysis based on the characteristics of 1100 test time series, with 100 from each land cover class (except nival areas). First, we examined whether the single years showed significant greening, significant browning, or no significant change. We ascertained that 55% have no significant change in each of the 32 analyzed years. In 36.7%, MSAVI increased significantly in all years-i.e., we detected a continuous linear trend over the entire study period. A continuous browning trend was revealed for less than 1% of the time series. Non-linear shapes-i.e., a mix of greening, browning and insignificant change over the years-were found for 7.7% of the time series. However, the individual land cover classes partly exhibit considerable differences from this result ( Table 2). Nonlinear time series occur predominantly in mountain Tugai (19%) and cultivated land (16%). Another remarkable land cover class is the mountain steppes. They have the highest proportion of no change in all the years (76%), followed by cryophytic and subnival vegetation (63%). Constant greening over the entire study period was especially detected in juniper vegetation (46%) and Rosaceae scrubs (46%), followed by barren land (43%). Then, we compared the amount of change-i.e., the difference in MSAVI of the component smooth functions for the variable 'year' between 1988 and 2020 and between the different types of times series based on their median values: for all non-NA time series, the median is 0.01, for time series with a significant greening trend it is 0.02, and for time series with a significant browning trend it is −0.03 (Table 3). To put these rather small values into perspective, we need to consider that the median MSAVI calculated over all 1100 test plots for the whole time series is only 0.06. The comparison of the absolute values did not unveil patterns of change between the land cover classes. In such cases, a common approach is to compare the classes in terms of some measure of location [124,125]. We calculated the quartiles that are standard statistical location parameters to compare the differences between land cover in the four different categories in terms of the amount of change, as described below.
First, we calculated the quartiles of all change values in the test plots with a significant positive change (n = 486). We call values above the 75% quantile a 'very strong change', values above the median but smaller than the 75% quantile 'strong change', values above the 25% quantile to the median 'intermediate change', and values below the 25% quantile 'weak change'. Then, we compared the land cover classes based on their share of change values in the four categories ( Figure 5). In particular, cultivated land has many plots with a very strong positive change, followed by mountain deserts. Mountain meadows is the land cover class with the highest number of plots with strong positive change, followed by floodplain meadows. The category intermediate change is led by juniper vegetation, followed by Rosaceae scrubs. Plots with weak positive change are predominantly found in juniper vegetation, followed by mountain Tugai and barren land. many plots with a very strong positive change, followed by mountain deserts. Mountain meadows is the land cover class with the highest number of plots with strong positive change, followed by floodplain meadows. The category intermediate change is led by juniper vegetation, followed by Rosaceae scrubs. Plots with weak positive change are predominantly found in juniper vegetation, followed by mountain Tugai and barren land.   Analogous to the plots with a significant positive change, we calculated the quartiles of the modulus of the change values in the test plots with a significant negative change (n = 26). Again, we compared the land cover classes based on their share of change values in the four categories ( Figure 6). There are no plots with a significant negative change among Rosaceae scrubs, cushion plant vegetation, and cryophytic and subnival vegetation. The highest number of negative changes was found in mountain Tugai, though most plots fall into the category of weak negative change. Very strong and strong negative changes predominantly occur in mountain meadows, followed by cultivated land. = 26). Again, we compared the land cover classes based on their share of change values in the four categories ( Figure 6). There are no plots with a significant negative change among Rosaceae scrubs, cushion plant vegetation, and cryophytic and subnival vegetation. The highest number of negative changes was found in mountain Tugai, though most plots fall into the category of weak negative change. Very strong and strong negative changes predominantly occur in mountain meadows, followed by cultivated land. With regard to topography, we found significant differences (p < 0.05, Mann-Whitney U-test) between the trend types for all analyzed variables (Table 4). Regarding elevation, there is a significant difference between all three trend types, with the highest values discovered for time series with a significant browning trend, followed by time series with no significant trend and time series with a significant greening trend. The slope values also differ significantly between the three trend types, with the steepest slope for time series with significant greening, followed by time series with no significant trend and time series with significant browning. Furthermore, the time series of the three trend types differ in terms of their north exposedness. The highest values could be detected for time series with a significant browning trend, followed by time series with a significant greening trend and time series with no significant trend. Moreover, the time series of the three trend types show a significantly different east exposedness. The highest values appear in time series with no significant trend, followed by time series with a significant browning trend and time series with a significant greening trend. Considering the ruggedness of the terrain, we also found significant differences between the time series of the three trend types, with the highest values for greening time series, followed by no trend time series and browning time series. Greening time series are the trend type closest to valley bottoms, followed by browning time series and no trend time series. With regard to topography, we found significant differences (p < 0.05, Mann-Whitney U-test) between the trend types for all analyzed variables (Table 4). Regarding elevation, there is a significant difference between all three trend types, with the highest values discovered for time series with a significant browning trend, followed by time series with no significant trend and time series with a significant greening trend. The slope values also differ significantly between the three trend types, with the steepest slope for time series with significant greening, followed by time series with no significant trend and time series with significant browning. Furthermore, the time series of the three trend types differ in terms of their north exposedness. The highest values could be detected for time series with a significant browning trend, followed by time series with a significant greening trend and time series with no significant trend. Moreover, the time series of the three trend types show a significantly different east exposedness. The highest values appear in time series with no significant trend, followed by time series with a significant browning trend and time series with a significant greening trend. Considering the ruggedness of the terrain, we also found significant differences between the time series of the three trend types, with the highest values for greening time series, followed by no trend time series and browning time series. Greening time series are the trend type closest to valley bottoms, followed by browning time series and no trend time series.
Finally, with regard to the geographic position as indicated by Universal Transverse Mercator (UTM) easting and northing, we also detected significant differences between all three trend types. For UTM easting, the highest values occur in browning time series, followed by no trend and greening time series. Regarding UTM northing, the values of greening time series are significantly higher than the values for no trend and browning time series, neither of which differ from each other.

Trends of Temperature and Precipitation
The visual inspection of the ERA5 and MERRA annual mean temperature values shows increasing trends for both time series (Figure 7 left). This observation is supported by the results of the GAMMs, where the component smooth functions of the variable 'year' revealed a significant (p < 0.01) increasing trend from 1988 to 2020 for all three analyzed temperature time series (CDO Khorog, ERA5, MERRA). The CDO time series, however, has many missing values; therefore, its result should be interpreted with caution, and meaningful annual means could not be calculated and plotted. The biggest temperature increase of 1.45 • C over the study period was indicated by the MERRA time series, followed by the ERA5 time series with 0.95 • C and the CDO time series with 0.8 • C.  In contrast, we could not find consistent trends in the precipitation data (Figure 7  right). The results of the GAMMs for the ERA5 and MERRA time series also support the assumption of no significant change over the study period (p > 0.05). Only the GAMM for the CDO precipitation data of the meteorological station in Khorog exhibits a significant increase in precipitation after 2007. However, this result is influenced by many missing values before this date.

Trends of Population Development
The population in the study area and its surroundings increased during the study period ( Figure 8). This is visible in the time series for the entire GBAO from 1990 to 2018, including the city of Khorog, and for the rural population (2000-2018). The numbers steeply increased until 2001, then slightly decreased from 2003 to 2007, and finally further increased after 2008. The smooth functions in the GAMMs confirm a significant (p < 0.01) positive trend. The population numbers  in the districts of Rushan, where the major part of the study area is located, and Shugnan/Roshtqala, the district of the southern part of the study area, show a positive trend as well. Rushan's population significantly increased until 2001, then remained stable, and lastly slightly increased after 2006. In Shugnan/Roshtqala, the numbers increased strongly until 2000, then remained stable on a high level, declined in 2010, and again increased in recent years.

Trends of Large Livestock Numbers
The numbers of large livestock in the study period first sharply decreased, at least until 1996, but exhibited a strong increasing development after 2000 (Figure 9). The time Figure 8. Population time series: the entire province of GBAO (total and rural population, left) and the districts of Shugnan/Roshtqala and Rushan (right; * indicates significant trends, p < 0.01; data sources: [111,112]).

Trends of Large Livestock Numbers
The numbers of large livestock in the study period first sharply decreased, at least until 1996, but exhibited a strong increasing development after 2000 (Figure 9). The time series for the entire GBAO from 1981 to 1996 shows stable numbers of around 42,000 large livestock heads until the end of the 1980s, followed by a significant (p < 0.01) decrease until 1996 to 19,440 animals. The GBAO time series from 1992 to 2011 exhibits a significant (p < 0.01), near linearly increasing trend from 2000 to 2011. While the numbers stay below 80,000 in the 1990s, 105,011 animals were counted in 2011. The decline in the livestock numbers at the beginning of the 1990s was also detected in the time series of the districts of Shugnan/Roshtqala and Rushan, where the numbers decreased from 11,814 to 3136 and from 4113 to 1610, respectively. The increase from 6904 to 8677 animals in Rushan after 2000 appears rather flat in the visualization but is nevertheless significant (p < 0.01). For Shugan/Roshtqala, we detected a significant steep increase from 2000 to 2003, followed by a further moderate increase.

Discussion
To our knowledge, this is the first study that successfully uses GAMMs to analyze vegetation and land cover change with a raster time series in complex dry mountain terrain. By combining this approach with a relatively long Landsat time series and a particularly suitable vegetation index, we addressed several common challenges of applied remote sensing studies in dryland mountain regions, such as irregular data, high background influence, and non-linearity.

Detection of Greening and Browning Trends
The analysis of the 32-year pixel-based MSAVI time series revealed significant land cover changes in about 40% of the study area from 1988 to 2020. Concerning the direction of change, the area exhibiting a greening trend by far exceeds the area with a browning trend. This finding coincides well with observations on a global scale, where trends of increasing green vegetation cover were detected in many regions of the world in recent decades [3,4]. About 60% of the analyzed area showed no significant land cover change.

Distinction of Greening and Browning Trends
We detected greening in all land cover classes, though the largest proportions were found in areas at or near the valley bottom. The amount of change was highest in cultivated areas. This finding is not in line with our expectation that valley bottoms are characterized mainly by browning due to overuse [49,51,59,63]. The highly elevated land cover classes show a smaller percentage of greening; however, due to their large extent, the absolute areas are still extensive. Cryophytic and subnival vegetation have a high proportion of time series without significant change, which shows the stability of this high-elevation land cover class. This is also an unexpected finding, as it contradicts the results of other studies that expect greening, particularly in the highest altitudes [43,44,126].
Browning could be detected for less extensive areas, though it also occurred in all land cover classes and predominantly in vegetation at or near the valley bottom. These time series mostly show a non-linear development, indicating abrupt changes. A cluster of such non-linear negative development was found in Red valley in the northwest of the study area, which resulted from a debris flow in 2011 (Figure 10a).  For the variable slope, browning sites have the lowest median, which can be explained by the highest impact of disastrous events in the relatively weakly inclined valley bottoms, whereas greening was detected in both steep and flat areas. The same explanations hold true for north-exposedness, for which the median in sites with browning is the highest. Again, this is because a relatively large cluster of browning was found in a north-exposed valley. However, we assume that north-exposedness is not a general explanation for browning trends, as similar events also occurred in south-exposed valleys outside the study area-e.g., around the village of Barsem [94,127].

Relationship of Greening and Browning Trends to Temperature and Precipitation
For mountains, in particular, rising temperatures and a prolonged growing season were discussed as major drivers for increasing vegetation cover and species richness in alpine areas. As the temperature is one of the major limits on plant growth, warming reduces this constraint [43,44,128]. For the study area, we found significantly rising temperatures. This result is in line with the findings of Haag et al. [129,130]; the review of Robinson [59]; δ18O data from ice cores [131]; and the Intergovernmental Panel on Climate Change (IPCC) report, which expects an increase of 3.7 • C in the mean annual temperature (temperature change for 2080-2099 in comparison to [1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990][1991][1992][1993][1994][1995][1996][1997][1998][1999] in the region [132]. We conceive these results as a strong indication that warming also plays a role in the greening of the study area. In most of the pixel time series that exhibit a greening trend, the increase in MSAVI is linear. This finding corresponds well to the nearly continuous temperature increase and corroborates the role of rising temperatures as a driver of greening. In arid areas, increasing precipitation is a frequent driver of greening [6,11]. However, we did not find significant trends of change in precipitation for the study area. This result is backed by the findings reviewed by Robinson [59], whereas Deji et al. [131] found a wetting trend in δ18O ice core data, and Haag et al. [129] detected increased summer precipitation.
In the highly elevated land cover classes, the increasing temperature that causes permafrost thawing is likely to contribute to the observed greening, particularly in moraine material that is classified as barren land, as well as in cryophytic and subnival vegetation. Peng et al. [133] investigated permafrost-vegetation relationships in the northern hemisphere and were able to show that the increased number of thawing days prolongs the growing season and increases the active layer thickness, which facilitates vegetation growth in permafrost regions. In the Pamirs, this is an important point, as 84.1% of the GBAO was identified as a potential permafrost area. Projections based on a scenario of a temperature increase of 4 • C until the end of the century expect a disappearance of 40% of the recent permafrost areas in the GBAO [134].
Furthermore, thawing permafrost decreases slope stability and hence promotes gravitational mass movements. Together with increased run-off induced by intensified melting of snow and glaciers [135][136][137][138][139], this may create catastrophic mass flows that devastate the valley bottoms and thus lead to an abrupt browning signal in remote sensing data [127,134,140]. Indeed such events of gravitational mass movement regularly occur in high mountains [127,140]. However, the prospective consistent temperature increase within the next few decades [132] will accelerate such events. Additionally, heavy summer rains in the study area in recent years resulting from disturbances triggered by the Indian Summer Monsoon or originating from extratropical cyclones [139,[141][142][143] exacerbate these processes. In some cases, the enhanced melting can lead to the development of glacial lakes and subsequently to glacial lake outburst floods (GLOFs). Such an event occurred in Shakhdara valley, ca. 50 km south of the study area, where a GLOF destroyed the village of Dasht in 2002, killing dozens of people and devastating vegetation [127,134].

Relationship of Greening and Browning Trends to Population and Livestock Numbers
The strong greening trend of the land cover classes at and near the valley bottom can be explained by the increased cultivation efforts of a rising population size together with improved land management practices linked to privatization since the breakdown of the Soviet Union [59]. Another argument in favor of this explanation is that most of the greening areas are located in regions where the highest population numbers are reported. Furthermore, the restoration and extension of irrigation systems are likely to play important roles. After a major extension of these systems during the Soviet period, their maintenance ceased with independence. However, many of these systems have been reestablished-e.g., as community-based irrigation governance systems (see, e.g., [144]). Similarly, projects of joint forestry management have helped to restore the Tugai scrubs and forests [145].
In mountainous drylands, floods and debris flows are often disastrous, as arable and irrigation lands, as well as associated villages, are predominantly situated on narrow river terraces and alluvial plains in valley bottoms [144]. In this area, the Tugai scrubs and forests play the dominant role in river discharge regulation and embankment stability, as well as in the protection of soils and infrastructure. Directly after the Soviet breakdown, these forests strongly degraded because of the energy crisis that forced the local population to use the Tugai for fuelwood and timber, leading to the increased vulnerability of the valley vegetation, arable land, and local communities [144][145][146][147][148]. In some areas, the state of these forests has already improved, but there is still an increased vulnerability to natural disasters in many areas due to poor forest conditions [149].
Another factor considered to be an important driver in browning is grazing impact. Our results show that the livestock numbers in the study area declined after the Soviet breakdown but strongly increased in the last two decades. This finding coincides with the results of Robinson [59]. Furthermore, current trends towards a reduction in livestock mobility and the abandonment of remote pastures are causing overgrazing in easily accessible areas and, in particular, around villages [59,62,150]. The increasing livestock numbers are primarily driven by the increasing population that we found for GBAO since 1990. In particular, the outbreak of the Tajik civil war in 1992 is seen as an important booster initiating the rising numbers after the Soviet breakdown [63]. A higher population also increases the need for other natural resources, such as fuelwood and timber. In the study area, the Tugai forest is the most important provider for these ecosystem goods, which explains the browning trend seen in some parts of this land cover class.
Indeed, we detected browning in land cover classes extensively used for grazing livestock-e.g., in meadows-though the extent is low. Mountain steppes, another heavily grazed land cover class, exhibit the largest proportion of pixel time series without significant change (for an illustration, see Figure 10b).
According to our observations, grazing caused no changes in vegetation quantity here but rather in vegetation quality. An explanation can be found in the long grazing history of the study area and the change to vegetation dominated by prickly shrubs and herbs or plant species that contain alkaloids [96,151]. Our field data support these findings. For example, we found many unpalatable species, such as Senecio spp., and we frequently detected a strong dominance of prickly herbs-in particular, of the genus Cousinia-in areas where historical vegetation studies reported a much lower quantity of this plant type [97]. We conclude that the change in the mountain steppes is not visible in the remote sensing data but is evident on species and trait levels.
The inspection of 76 ground-truthing sites in the field and comparison with the historic vegetation map support the plausibility of the findings derived from the remote sensing analyses in most cases. However, the results for scree sites with no vegetation should be interpreted with caution, as we detected unreliable steep positive trends in some rare cases (see the example given in Figure 10c).

Conclusions
In this paper, we presented an innovative and promising approach for analyzing land cover trends in dry mountain systems by combining long-term Landsat image time series with environmental and socio-demographic data using GAMMs. Thereby, we were able to tackle several challenges of remote sensing studies in mountainous drylands, such as irregular data, high background influence, non-linearity, and sensor issues. We found that greening and browning trends coexist in the study area, with greening being predominant. The characteristics of the detected trends vary according to the vegetation and land cover classes involved, as well as their topographical and geographical distribution. With the applied approach, we were able to detect remarkable changes at a small spatial scale-e.g., in narrow valley bottoms that remained undetected in previous long-term remote sensing studies. These areas are predominantly used for agriculture and forestry, which suggests direct anthropogenic drivers of change here. We also discovered a high number of pixels with a greening trend that corresponds well with increasing temperature. Trends of increasing precipitation, which would support greening, were not observed. Detected browning was mainly attributed to disastrous events that were most probably enhanced by increasing temperatures. Contrary to our expectations, browning trends did not correspond well with an increasing population and increasing livestock numbers. Grazing increased, but associated changes were not detectable with our approach. This can be explained by our field observations, which showed that the respective effects are most active on the species and trait level. Different approaches and future studies are required in order to enhance our understanding of the processes involved at this resolution. Nevertheless, this study presents trends of land cover change in the Pamir mountains with a novel methodology that could serve as a basis for further research on vegetation dynamics for management purposes-e.g., for sustainable pasture and forest management. Thereby, the approach shows high potential in mountainous and dry areas globally.
Author Contributions: K.A.V. designed the study, conducted the analysis, and wrote the manuscript; H.Z. and C.S. revised the manuscript; K.A.V., H.Z. and C.S. discussed the results and reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript.