Detecting abrupt change in land cover in the eastern Hindu Kush region using Landsat time series (1988–2020)

Land cover change in the semi-arid environment of the eastern Hindu Kush region is driven by anthropogenic activities and environmental change impacts. Natural hazards, such as floods presumably influenced by climatic change, cause abrupt change of land cover. So far, little research has been conducted to investigate the spatiotemporal aspects of this abrupt change in the valleys. In order to explore the abrupt change in land cover and floods as its possible drivers in the eastern Hindu Kush, a semi-arid mountain region characterized by complex terrain, vegetation variation, and precipitation seasonality, we analyzed long-term Landsat image time series from 1988 to 2020 using Breaks For Additive Seasonal and Trend (BFAST). Overall, BFAST effectively detected abrupt change by using Landsat-derived Modified Soil Adjusted Vegetation Index (MSAVI). The results of our study indicate that approximately 95% of the study area experienced at least one abrupt change during 1988–2020. The years 1991, 1995, 1998, 2007, and 2016 were detected as the peak years, with the peaks occurring in different seasons. The annual trend of abrupt change is decreasing for the study area. The seasonality of abrupt change at the catchment level shows an increasing trend in the spring season for the southern catchments of Panjkora and Swat. The spatial distribution patterns show that abrupt change is primarily concentrated in the floodplains indicating that flooding is the primary driver of the land cover change in the region. We also demonstrated the accurate detection of past flood events (2015) based on the two case examples of Ayun, Rumbur, and Kalash valleys. The detection of the flood events was verified by fieldwork and historical high-resolution Google Earth imagery. Finally, our study provides an example of applying Landsat time series in a dry mountain region to detect abrupt changes in land cover and analyze impact of natural hazards such as floods.

a geo-ecologically important mountain region as it provides ecosystem services to 240 million people living in this area (Wester et al. 2019). HKH region is prone to several natural hazards (Rusk et al. 2022). It is highly susceptible to geological hazards such as earthquakes, landslides, rockfalls, and erosion, mainly because it is a relatively young geological formation (Hewitt 1992;Kamp et al. 2004;Nibanupudi and Shaw 2015). In addition to geological threats, hydrometeorological hazards pose a significant risk to the people living in the HKH region. Floods constitute a considerable danger throughout the entire area (Elalem and Pal 2015;Tsering et al. 2021). Riverine and flash floods are common recurring disasters Nibanupudi and Shaw 2015;Tsering et al. 2021). In particular, on the slopes of the highly elevated upper catchment areas of rivers, sudden and intense rainfall causes flash floods (Attaur-Rahman and Khan 2011;Lu et al. 2022). Such events are also supplemented by glacial melting and heavy snowfall, resulting in flash floods in the downstream valleys. The sudden bursting of glacial lakes also causes outburst floods in the valleys (Allen et al. 2016;Ashraf et al. 2021;Sati 2022;Veh et al. 2018). The presence of many glacial lakes in the mountainous region increases the flood risk for the settlements in the floodplains (Ashraf et al. 2014;Ashraf et al. 2021;Rehman et al. 2014). Moreover, riverine flood events also destroy farmlands, irrigation infrastructure, housing, roads, vegetation, and livestock (Ashraf et al. 2021;Khan 2011, 2013;Hewitt and Mehta 2012;Khalid et al. 2021;Nüsser 2001). The impacts of natural hazards further exacerbate societal problems, including the out-migration from mountain areas of HKH (Gioli et al. 2014;Maharjan et al. 2021;Rusk et al. 2022).
Several studies have investigated the impacts of flood events in the HKH region (Ashraf et al. 2012;Ashraf et al. 2014;Sati 2022;Shan et al. 2021;Veh et al. 2018). Damage assessments were carried out using remote sensing data and techniques. For instance, the impacts of the 2010 floods in Punjab and 2011 floods in Sindh provinces of Pakistan were analyzed (Gaurav et al. 2011;Haq et al. 2012). In addition, research has been conducted on spatiotemporal changes in forest cover in the HKH region (Haq et al. 2018;Ullah et al. 2016;Zeb 2019). However, there is a lack of research on the long-term spatiotemporal analysis of abrupt change in land cover in the valleys. Abrupt changes in land cover are generally caused by anthropogenic activity (e.g., deforestation and urbanization etc.) and natural processes such as floods and landslides (Chen et al. 2014;Xu et al. 2020). These abrupt changes can be driven by many factors, including past events' impact (Chen et al. 2014;Fang et al. 2018;Watts and Laffan 2014). The impact and timing can vary in different catchments for different seasons. Therefore, both seasonal and catchment level analyses of abrupt changes are crucial to understand land cover change dynamics in the region.
Remote sensing methods provide an ample opportunity to detect land cover changes in remote, complex mountain regions and on a large scale (Anderson et al. 2020;Geng et al. 2019;Mishra and Chaudhuri 2015) and have been widely used in land cover change detection (Brown et al. 2006;Brown et al. 2012;Running et al. 1995;Tian et al. 2016). The Landsat program offers satellite imagery with a long historical record and is very suitable to detect past changes in the land cover (Campbell 2011;Young et al. 2017;Zhu 2017Zhu , 2019. In particular, time series analysis of vegetation indices derived from Landsat sensors has been extensively applied to analyze changes in land cover and their causes (Franklin et al. 2015;Panday and Ghimire 2012;Zhu 2017Zhu , 2019. Recent literature shows that Landsat time series have been used in studies on glacier dynamics, detecting glacial lakes, their expansion, and associated floods in different mountain regions such as HKH, Central Asia, and European Alps (Huang et al. 2018;Ma et al. 2021;Veh et al. 2018;Yan et al. 2018;Zheng et al. 2019). In addition, changes in a forest in wetland areas of China were investigated, and flood impact on Swiss floodplains was also analyzed (Milani et al. 2020;Wu et al. 2020). However, we found that no studies are using Landsat time series to analyze abrupt changes in a dry mountain region, such as the eastern Hindu Kush. Our study aims to fill this research gap by using the Landsat time series  to detect abrupt change in land cover in dry mountain regions, like the eastern Hindu Kush.
The main objective of this study was to investigate the abrupt changes in land cover in the eastern Hindu Kush and then compare their occurrence and impacts for the different catchments. This research attempts to answer the following specific aims: (1) to detect and map abrupt changes, including peak years, in the trend component of the Landsat time series; (2) to analyze detected abrupt changes at the catchment level; (3) to extract past flood events in the valleys and map their timing and magnitude of change in the land cover. As such, our study represents a prime example for analyzing land cover change using raster time series analyses in dry mountain area that could also be applied to similar regions.

Study Area
The eastern Hindu Kush region is surrounded by the Karakoram and Himalayan ranges in the east and southeast and the Pamir mountains in the north and northwest. Our study area is situated in the north of Khyber Pakhtunkhwa province and western part of Gilgit-Baltistan in Pakistan, and eastern parts of Nuristan and Badakhshan provinces of Afghanistan, between 35°20'24' ' and 36°43'48'' N, and 71°19'12'' and 72°59'24'' E. It has an extent of approximately 22,000 km 2 (Fig. 1). The north of the study area is part of the Amu Darya basin with Wakhan and Warduj catchments. The central part drains into the Chitral catchment, while the Basghal catchment lies on the south-western side. The south of the study area has Panjkora, Swat, and Kandia catchments, while the Gilgit catchment is in the east.
Chitral River (also called Kunar) is the main river in the eastern Hindu Kush. In addition to Chitral, other rivers in the study area are Panjkora, Swat, Gilgit, Kandia, Wakhan, Warduj, and Bashgal. All these rivers are tributaries of the Upper Indus River Basin except the Warduj and Wakhan rivers which are part of the Amu Darya in Afghanistan and Central Asia. Tirich Mir (7,708 m), the tallest peak in the Hindu Kush range, is located among other mountains and glaciers in the study area. Most of the settlements are along the riverbanks of the Chitral River and other tributaries. The majority of the population is rural, while Chitral town (35°51'17.19" N, 71°47'21.74" E), with about 49,000 inhabitants, is the largest urban settlement located in the southwest of our study area (Fig. 1). Chitral weather station recorded an annual mean temperature of 16.1°C and total annual precipitation of 405 mm from 1970 to 2019 (PMD 2020). The north and northwest part of the study area receives less rain with a peak in winter and spring, while the south and south-eastern side receive more monsoon rains in summer (Nüsser and Dickoré 2002). Forests are located in the southern part, whereas central and the northern Chitral is largely treeless (Nüsser and Dickoré 2002;Zeb 2019).

Landsat data
We choose the Modified Soil Adjusted Vegetation Index (MSAVI) due to its strength to minimize soil background effects caused by the arid and semi-arid conditions, sparse vegetation, and the existence of bare soils in several valleys and surrounding mountains in our study area (Qi et al. 1994). The MSAVI is derived from the reflectance in the nearinfrared (NIR) and red bands and is denoted by the following equation: We acquired MSAVI data derived from 811 images of the sensors Landsat 4-5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI) for Path 151 and Row 35 (World Reference System-2) for the years 1988 to 2020 (Table 1). We selected a full-season archive of images from 1988 to 2020 without applying any seasonal or time-specific filters. The MSAVI data is processed to Level-2 Surface Reflectance data, which means it has undergone radiometric, geometric, and atmospheric corrections by USGS (2015). This data comes with a spatial resolution of 30 m and a temporal resolution of 16 days. It was accessed from the United States Geological Survey's (USGS) Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) platform (https://espa.cr.usgs.gov/).

Digital elevation model
To derive the river network and catchments in our study area, we used the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) Version 3 with a spatial resolution of 30 m. It was downloaded from USGS Earth Explorer (https://earthexplorer.usgs.gov/).

Validation data
To validate the results of flood detection in this study, we gathered data for selected events through fieldwork in the Lower Chitral district in August 2020. This entailed a review of information on these events in disaster reports of local organizations and government (PDMA 2015;Shah 2015Shah , 2016. The fieldwork also included visiting Ayun (35°43'29.38" N, 71°46'33.82 E") which was hit by a riverine flood in 2015, and Bhumboret (35°41'13.88'' N, 71°39' 55.13'' E) affected by a flash flood in the same year. We took photographs of the areas impacted by 2015 floods and gathered information from local people through unstructured interviews on the timing and location of these flood events. Furthermore, high-resolution historical imagery was accessed from Google Earth and geo-referenced.

Data pre-processing
Data pre-processing steps are illustrated in the flowchart (Fig. 2). The Landsat-derived MSAVI images were cropped for the study area, and clouds and cloud shadows were removed using the quality assessment band provided by USGS. Sinks in the DEM were filled, and the river network was generated using Strahler order. Since natural hazards (land erosion, debris flows, etc.) triggered by flood events are concentrated in the valleys along the streams, high slopes, glaciers, and areas remote from streams were masked out using buffers of 250 -1000 m width along the river network.

Elbow criterion test
Another important consideration before time series analysis was the minimum number of valid observations per pixel. The history of sparse observations by the Landsat program in the late 1980s and 1990s and the persistence of clouds and cloud shadows affected the number of valid observations. Pixels with few valid observations may lead to the detection of false breakpoints. To determine an appropriate threshold of minimum valid observations, we carried out the Elbow criterion test by setting all pixel values to NA if the number of valid observations is less than this threshold (Ketchen and Shook 1996). We applied this test on seven tiles (out of 100) in our study area with a range of thresholds between 0 and 500. We noticed that valid observations up to a threshold of 200 are stable across the tiles but reduce significantly if the threshold is increased (Fig. 3). We assume that a minimum of 200 valid observations in the time series is sufficient to detect valid breakpoints for a pixel. All pixels with less than 200 observations are turned to NA (non-applicable) and hence excluded from the breakpoint detection.

Abrupt change detection with BFAST
We used BFAST (Breaks for Additive Season and Trend) for time series analysis of Landsat data from 1988 to 2020 in our study region in the eastern Hindu Kush. BFAST detects abrupt changes on top of seasonal variations and possible trends (Verbesselt et al. 2010;Watts and Laffan 2014). It is robust against noise and is not influenced by the changes in the seasonal component of the time series (Verbesselt et al. 2010). BFAST has been used to detect abrupt and long-term changes in areas such as deforestation, vegetation, land use, land cover changes, and lake monitoring (Geng et al. 2019;Gholamnia et al. 2019;Lambert et al. 2013;Wu et al. 2020). It has also been applied to detect vegetation changes for selected pixels affected by flood events (Watts and Laffan 2014). However, BFAST has not been used to detect and map abrupt land cover changes in arid and semiarid environments like the eastern Hindu Kush. This high mountain region has a complex terrain, varying vegetation cover, and strong seasonality of the precipitation (Ahmad et al. 2018;Nüsser and Dickoré 2002). In addition, the remote sensing imagery of the area contains noise due to persistent cloud cover, especially in the winters. In such conditions, BFAST can be a suitable method because of the aforementioned robustness to noise and strong seasonal amplitudes.
BFAST decomposes a time series into three components: trend, seasonal, and the remainder (or noise) and does not require any historical period or the definition of a trajectory (i.e., a presumption on the expected trend of the data) (Verbesselt et al. 2010). The following equation represents the general model:

Fig. 3
Elbow criterion test used to determine minimum valid observation for a pixel. The plot shows a reduction of valid observations (y-axis) for each tile, with an increasing threshold (x-axis).
where is the observed data at time , is the trend component, is the seasonal component, is the noise or remainder component and is the total number of observations in the time series. The remainder component contains the remaining variation in the data beyond the seasonal and trend components (Verbesselt et al. 2010).
BFAST uses an additive decomposition model which iteratively fits a piecewise linear function with breakpoints * , … , * and defined * = 0, therefore the trend component can be expressed as: where = 1, … , and * < ≤ * . The magnitude and direction of the abrupt change at a breakpoint is derived using the intercept ( ) and slope ( ) of consecutive linear models by calculating the difference between at * and * , as expressed below: BFAST also uses the piecewise linear seasonal model on a seasonal dummy variable to fit the seasonal component (Verbesselt et al. 2010). The timing of breakpoints detected in the seasonal component can differ from those detected in the trend component (Verbesselt et al. 2010). If the seasonal breakpoints are # , … , # , and # = 0 , then for # < ≤ # , the seasonal component is expressed by following equation: where and denote the period of seasonality and effect of season , respectively. Since the sum of the seasonal components for all successive times is zero for # < ≤ # , can be presented as: where , = 1 when is in season and 0 for otherwise. In case, is in season 0, then , − , = − 1. For other seasons, when is in season ≠ 0 then , − , = 1. Moreover, , is a dummy seasonal variable with 0 and 1 values to account for seasons in a regression model (Verbesselt et al. 2010).
BFAST iteratively fits the seasonal-trend decomposition procedure to detect breakpoints through the following four steps; (i) the ordinary least squares residuals-based moving sum (OLS-MOSUM) test detects the breakpoints in the trend component, the number and position of breakpoints are computed from the seasonally adjusted data; (ii) the trend coefficients are then calculated through a robust regression that is based on M-estimation; (iii) in case OLS-MOSUM test detects breakpoints in the seasonal component, their number and position are estimated from trend adjusted data; and (iv) the seasonal coefficients are calculated using a robust regressionbased on M-estimation (Verbesselt et al. 2010). For a more detailed explanation of the BFAST method, readers are referred to the original paper (Verbesselt et al. 2010).
This study focused on the time series' trend component to explore the abrupt changes. We implemented BFAST package version 1.5.7 (Verbesselt et al. 2014) in the R statistical software version 3.6.1 (R Core Team 2020). The irregular time series were extracted, converted to daily and then to monthly time series. We used a seasonal dummy model on a median monthly time series of Landsat data (Verbesselt et al. 2010). This method derives the seasonal component to detect changes in the trend component (Watts and Laffan 2014). An important parameter to select in implementing BFAST is the ℎ parameter, which determines the minimum segment size, leading to the number of potential trend segments and the number of breakpoints in the trend component (Watts and Laffan 2014). The higher values of ℎ lead to the omission of several breakpoints, and low values may lead to the detection of false breakpoints. Therefore, we tested different h parameter values on various locations with a known history of flood events to determine the optimal h value. Finally, we assumed h=1/9 (i.e., 44 months as minimum segment size) as the appropriate value for this parameter used in this study. The time series results were the number of breakpoints for each pixel, the timing and the magnitude of these breakpoints.

Post-processing and validation of floods
Since we focus on floods that disrupt the vegetation in the valleys resulting in negative abrupt change, we extracted breaks with negative change from BFAST results (Fig. 2). A mosaic was generated from the tiles, and seasonal analysis was carried out at the catchment level. The detected abrupt change was further analyzed to extract specific disaster events. We demonstrated this approach successfully by exploring known flood events in the village of Ayun (35°43'29.38" N, 71°46' 33.82" E) and three adjacent valleys of Ayun, Rumbur and Kalash (35°39'36' '-35°49'1.2'' N, and 71°37'58.8''-71°49'55.2'' E) in Lower Chitral. The timing and magnitude of the flood events in these sites were verified using spatiotemporal information collected during the fieldwork. The results were also confirmed through visual interpretation using high-resolution Google Earth imagery before and after the event and documented by photographs taken during the fieldwork.

Number of breakpoints
The number of breakpoints detected in the time series' trend component varied between zero and eight for each pixel, with most pixels exhibiting between two and four breakpoints (Fig. 4). Approximately 4.7% (273,139) of the pixels in the valleys did not experience any change, while 95.3% (5,524,637) showed at least one abrupt change from 1991 to 2017.
The spatial distribution of the breakpoints detected in the time series' trend component is shown in Fig. 5. The statistical results show that 95% of the pixels in the Chitral (3,316,657), Wakhan (107,127) and Swat (498,098) catchments experienced at least one breakpoint. About 91% of the pixels in the Panjkora (410,378) catchment experienced abrupt change, which is the lowest among all the catchments. On the contrary, the Bashgal catchment had 98% pixels (323,555) with breakpoints, the highest among all catchments. About 96% of the pixels in Warduj (279,702) and 97% of the pixels in Kandia (222,447) and Gilgit (366,673) catchments experienced breakpoints.

Annual abrupt changes
Annual negative abrupt change in MSAVI for breakpoints is shown in Fig. 6. Peak years that experienced significant abrupt change are 1991, 1995, 1998, 2007, and 2016. In 1991 alone, approximately 1,160,613 pixels (20%) exhibited abrupt change making it the year with the highest number of   breakpoints. The second significant year was 1998 (1,096,131; 19% of the pixels) followed by 2007 (716,168; 12% of the pixels) as the third year with the most noticeable change. Eleven years had pixels with an abrupt change between 0.3 and 0.6 million each. The abrupt change is concentrated in different seasons for different years (Fig. 6). For example, abrupt change in 1991 was strongly concentrated in the autumn while, e.g., in 1998, spring and summer showed the highest number of breakpoints.
Overall, the annual trend of abrupt change is on the decline, but each season has different peak years (Appendix 1). For instance, the trend in the autumn season is in fact on the decline; 1991 is the year with the highest peak in this season. The winter season recorded peaks in 1991, 1998 and 2007, also with an overall downward trend. The spring season showed the top three peaks in 1995, 1998 and 2007 respectively with a downward trend, whereas for summer, 1991, 1995 and 1998 remained peak years with an overall decreasing trend.
The spatial distribution of peak years varies in the different catchments, as shown for the examples of 1991 and 2007. In 1991 (Fig. 7), the abrupt change scattered across eight river catchments in the study area. Approximately 25% of the pixels were affected in Swat (121,205 pixels), Kandia (59,399 pixels), Gilgit (93,330 pixels), Panjkora (112,575 pixels), and Bashgal (82,848 pixels) river catchments. The Chitral river catchment area showed about 19% of the pixels (653,539 pixels) that changed in 1991. Warduj (28,620 pixels) and Wakhan (9,097 pixels) river catchments in eastern Afghanistan were the least disturbed catchments (equal or less than 10% of the pixels). The magnitude of the change for breakpoints was between 0 and -0.5, suggesting a loss of vegetation cover. The negative magnitude of the breakpoint shows the strength of the loss of vegetation.
The timing and magnitude of abrupt change in 2007 are shown in Fig. 8. Most river catchments experienced abrupt change during the winter and spring seasons. Kandia (47,955 pixels) and Swat (113,391 pixels) catchments were the most disturbed (20% of the pixels), while Wakhan (8,833 pixels) and Warduj (24,484 pixels) were least changed (less than 10% of the pixels). Approximately 9% of the pixels changed in the Chitral (316,501) catchment, while 18% of the pixels exhibited abrupt change in the Panjkora (79,732) and Bashgal (58,746) catchments in 2007. The magnitude of abrupt change varied between 0 and -0.5, with a large portion of pixels falling in the range of 0 and -0.1. Seasonal analysis of abrupt change suggested the peak months varied across different catchments during 1991 -2017 (Fig. 9). November and December were the peak months in Bashgal, Panjkora, Kandia, Chitral, and Gilgit catchments during the 1991-2000 period. June was the peak month in Swat, Wakhan, and Warduj catchments in this decade. Warduj, Gilgit, and Wakhan catchments recorded peak in December during 2001-2010. Bashgal and Chitral maximum breakpoints were in April, while Panjkora, Kandia, and Swat catchments experienced peaks in February. During the 2011-2017 period, Warduj and Wakhan river catchments experienced extremes in February, while the highest breakpoints in Bashgal were in March. In the south and southeast of the study area, three adjacent catchments, Kandia, Swat, and Panjkora, showed their peaks in April. The Gilgit and Chitral river catchments' peaks were in May during the last decade. Overall, the number of breakpoints has declined significantly over the decades, especially from October to December.

Analysis of abrupt change at catchment level
Analyzed at the catchment level, the trends of abrupt change for summer, autumn and winter are decreasing for all the catchments (Appendixes 2-9). The trends in the spring season are decreasing for Warduj, Wakhan, Chitral and Gilgit catchments. On the contrary, Kandia and Bashgal have no trend while Panjkora and Swat catchments showed increasing trends in the spring season.

Detection of flood events
In the Chitral river catchment, a pronounced flood event that was detected in the three adjacent Ayun, Rumbur and Kalash valleys (35°39'36''-35°49'1.2'' N, and 71°37'58.8''-71°49'55.2'' E) was used for verification. The timing and magnitude of the flood impact extracted from our results are shown in Fig. 10. The event's timing suggests that the flood occurred in July 2015 and triggered breakpoints in all three valleys. The magnitude of abrupt change was between 0 and -0.4. This flood event damaged housing, crops, bridges, and infrastructure along the streams in the valleys.
This result could be confirmed by the interviews and the visual inspection of Google Earth imagery. Fig. 11(a), dated 27/08/2010, shows the state of vegetation, river and other infrastructure before the flood event in a trench of Kalash valley. Fig. 11(b), dated 24/09/2019, which is after the July 2015 floods, clearly shows the widening of the riverbanks and damage to vegetation and infrastructure by flood along the riverbanks. The devastation caused by this flooding in Batrik village of the Kalash valley was documented by field photographs (Fig. 12).
Furthermore, we investigated the impact of floods in the village of Ayun (35° 43' 29.38" N, 71° 46' 33.82" E). It was affected by floods in the Chitral river in February and July 2015, resulting in the loss of cropland and trees. The area affected by each flooding spell is shown in Fig. 13(a), while the corresponding change in magnitude is shown in Fig. 13(b). In our results the timing and magnitude of the event were accurately detected. The negative change of MSAVI is attributed to the impact of the flood events due to the loss of vegetation cover caused by erosion. This result could also be confirmed using high-resolution historical imagery from Google Earth (Fig. 14) and field photographs (Fig. 15).     Fig. 11 (b).

Fig. 13
BFAST results produced for Ayun using MSAVI, NDVI, and NDWI. The vegetation loss due to floods was detected in February and July 2015 (a), with a reduced magnitude for MSAVI and NDVI and increased magnitude for NDWI (b). The location of the Ayun site is shown in Fig. 14.

Discussion
Breakpoints, affecting 95.3% of all pixels, are spatially distributed across all the catchments in the study area (Fig. 5). The high number of pixels with breakpoints suggests that vegetation cover has experienced extensive change. This change can be attributed to both anthropogenic and environmental factors. Besides the recurring flood events affecting the vegetation in the valleys along the eastern Hindu Kush rivers (Nibanupudi and Shaw 2015;Tsering et al. 2021), boulder fall is a widespread hazard as indicated by field observations but is unlikely to affect an entire pixel or several pixels. Similarly, landslides and avalanches are also recurring hazards in our study area. Still, they have a different spatial pattern than the abrupt changes by floods that occur simultaneously in connected areas along the valley floors. Anthropogenically induced changes include construction works, clear-cuts, and abandonment or redevelopment of land, which are expected to cause patchy and incoherent patterns.
We noticed that the upper catchments experienced more breakpoints than the lower catchments (Fig. 5). The high number of breakpoints is possibly due to the occurrence of flash floods in the upper catchments, which are usually caused by intense rainfall in addition to heavy snow and glacial melting (Ashraf et al. 2021;Nibanupudi and Shaw 2015). The pixels experiencing vegetation loss are aligned along the river courses, which suggests flooding as the cause. Smaller catchments in the upper valleys react more to locally heavy precipitation than to large-scale inundation, which is why we assume flash floods to be the cause there. This interpretation is supported by our field observations showing high sediment transport typical of hyper-concentrated flows and debris flows. Flash floods could cause more damage than riverine floods due to their high velocity, little or no early warning, and relatively higher concentration of debris than riverine floods. In the smaller valleys, fewer settlements are affected by flash floods because they are sparsely populated.
On the other hand, the lower catchments are susceptible to the riverine floods caused by water overflows in the rivers and their impact on crops, lands, infrastructure, and human lives in the floodplains. The breakpoints with a negative magnitude of change are concentrated in the valleys along the rivers rather than on the mountains' slopes which shows that flooding is, quantitatively, the most critical hazard in this area. This confirmed our presumption that the negative changes co-occur in connected areas along the rivers, while positive changes are more scattered with single pixels.
Negative abrupt change is scattered across the entire time series. 1991, 1995, 1998, 2007, and 2016 appeared as the peak years during the 1991-2017 period (Fig. 6). Breakpoints in these peak years are distributed across different seasons. The abrupt change is concentrated in autumn for 1991, in winter for 2007, and in spring for 1998 and 2016. The spatial distribution of abrupt change shows that the breakpoints were concentrated in southern catchments of Swat, Kandia, Panjkora, Gilgit and Bashgal in 1991. Kandia and Swat, the southern catchments, again were affected by the abrupt change in 2007. The region is characterized by complex climatology where various local and regional climatic conditions influence events in different seasons. The exact determination of the drivers of these seasonal peaks requires a sophisticated climatological approach beyond the scope of this study. In general, summer events are caused by summer monsoon rainfall extremes which are influenced by the Madden Julian Oscillation, Indian Ocean Dipole, and indirectly by El Niño Southern Oscillation (Gadgil et al. 2004). Winter and spring events are influenced by the Northern Atlantic Oscillation and the Siberian High (Syed et al. 2006). Spring events are possibly caused by heavy snowfall in winter and the early melting of glaciers and snow in spring (Ahmad et al. 2018;Shahid and Rahman 2021). Moreover, intense summer events are related to a strong summer monsoon (Webster et al. 2011). Interestingly, e.g., 2010 with a strong monsoon has comparably low abrupt changes in the summer season.
The trends for summer, autumn, and winter are decreasing for all the catchments, while mixed trends can be noticed in the spring season (Appendixes 2-9), with two catchments (i.e., Panjkora and Swat) showing increasing trends, and Bashgal and Kandia exhibiting no trend. This upward trend could be due to the different climatic conditions (relatively high rainfall, vegetation patterns etc.) in the southern part compared to other regions in our study area. Chitral catchment has diverse climatic conditions, with southern valleys receiving more rainfall than the valleys in the north. It may also have different seasonal trends for the southern and northern parts. Overall, the decreasing annual trend can be attributed to the decreasing annual rainfall trend in the Upper Indus Basin as argued by other hydrometeorological studies (Latif et al. 2018;Yaseen et al. 2020). The increasing trend in the spring season in some of the catchments is probably influenced by the increase in glacial and faster snow melt caused by a rise in temperature in spring (Baig et al. 2021;Shahid and Rahman 2021).
Through the spatiotemporal exploration of the BFAST results, we were able to trace the timing and spread of specific flood events, e.g., the flooding of 2015 in Ayun, Rumbur, and Kalash valleys (Fig. 10,  13). The results were interpreted using a mixed approach comprising visual interpretation and historical imagery and fieldwork. Overall, both the historical floods were reflected by major breakpoints, but there are many other isolated breakpoints for which no corresponding record was found. Natural hazards such as rockfalls and landslides are challenging to detect due to the lack of accurate information on their occurrence and smaller size. As stated earlier, these isolated events could also be triggered by anthropogenic activities such as deforestation and undercutting of the slopes for building roads and other infrastructure. Flooding can only be discerned from anthropogenic activities through its coherent patterns and proximity to the river network. This limitation can be overcome by using land cover classification to understand changes in land cover types. However, the BFAST method cannot be used on its own to explore the change in land cover classes. For future studies, we think of applying new techniques such as Prophet and dynamic time warping based time series classification approach and deep learning based time series classification that already showed promising results (Rußwurm and Körner 2020;Yan et al. 2019).
This study is the first attempt to explore abrupt change in the land cover using BFAST and Landsat image time series in the complex terrain (i.e., arid, semi-arid environment, high mountain region, diverse vegetation cover influenced by strong seasonality of precipitation) of the eastern Hindu Kush. Our findings indicate that the region and its catchments experienced significant, abrupt changes in land cover from 1988 to 2020. The methodological approach enables us to understand the spatiotemporal patterns and to extract past flood events. Moreover, our research method can furnish a flood disaster database for this area. In addition, it provides a historical spatiotemporal analysis of the past floods and their impact in various valleys. Finally, the results and methodology are meant to support policymakers and disaster management practitioners to prepare future hazard risk assessments and disaster management plans for the eastern Hindu Kush mountain communities.

Conclusions
In this study, the spatial and temporal patterns of abrupt change in land cover in the eastern Hindu Kush valleys were investigated using BFAST based on Landsat-derived MSAVI times series from 1988 to 2020. The study confirmed the applicability of the BFAST method to detect abrupt change in the complex environment of this high mountain region.
Our study showed that all three research questions were answered, which is reflected in the conclusions: 1. The study results indicate that only 5% of pixels remained stable, while 95% experienced a minimum of one breakpoint from 1988 to 2020. The occurrence of a high proportion of breakpoints suggests that vegetation has undergone extensive change in the eastern Hindu Kush region. Moreover, many breakpoints were detected in the upper catchments of streams and rivers, indicating their high susceptibility to flash floods. 2. 1991, 1995, 1998, 2007, and 2016 appeared as the peak years for abrupt change. There were a few years with heavy rainfall events which are not detected, e.g., the strong monsoon rains in 2010. At the same time, the peak years vary for each season. Peak years in different seasons are probably caused by a complex interplay of several climatic circulation patterns, which influence precipitation patterns in the region. Summer, autumn, and winter seasons showed downward trends for all catchments, while Kandia and Bashgal have no trend, and Swat and Panjkora catchments showed upward trends in the spring season.
3. The concentration of negative abrupt change in the floodplain suggests that floods are major natural processes that trigger debris flows and river erosion in our study area. These natural hazards cause damage and loss of agricultural land, housing, and other community infrastructure. Their impacts were shown by examples of the 2013 and 2015 floods in Kalash, Rumbur, and Ayun valleys of Chitral.
4. BFAST applied to Landsat-derived MSAVI time series data can be very useful to analyze past natural processes in high mountain regions such as the eastern Hindu Kush. The approach involving analysis of abrupt change at different levels, number of breakpoints at catchment level, annual and seasonal change, and the extraction of events in specific valleys enabled us to understand the land cover change in the region. 5. Future research should investigate the abrupt change in land cover by integrating BFAST results with social research (e.g., survey questionnaires and interviews) to explore the impacts of natural hazards on mountain communities. manuscript. Thanks also go to the two anonymous reviewers and editors whose comments and suggestions improved this manuscript.

Open Access
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
Funding note: Open Access funding enabled and organized by Projekt DEAL.