Forest Cover Changes and Trajectories in a Typical Middle Mountain Watershed of Western Nepal

There have been drastic changes in resource use practices and land-use patterns in the middle mountains of Nepal as a result of human transformation processes of the environment. This study aimed at assessing land-use and land-cover changes, especially those related to forest cover changes, in Phewa Lake watershed—a typical middle mountain watershed of western Nepal—using multi-temporal Landsat images from 1995, 2005 and 2017. Landsat images of each year were classified individually using object-based image classification into four land-use and land-cover types: agriculture and built-up, forest, waterbodies and other. Post-classification comparison was employed to quantify the extent and rate of changes, which was further extended to quantify the level of persistence, gains, losses, and swaps of forests. Furthermore, temporal trajectories of land-cover associated with forest cover changes were established, and their spatial pattern analyzed. The results show that, between 1995 and 2017, forest cover increased by 6.8% with a corresponding decrease in the extent of all other land-cover types. Dynamic transitions and internal trading among forest and agriculture and built-up category were observed, revealing more complex patterns than the commonly assumed linear and irreversible forest cover transformations in the mountains of Nepal. Our approach to assess major signals of forest cover transitions and change trajectories will help link patterns to the process of change including deforestation and forest regeneration. This would, in turn, form the basis for formulating practical conservation and management strategies for Phewa Lake watershed and other mountain watersheds of Nepal.


Introduction
Sustainable development of forest resources in Nepal has acquired added importance with the deepening anxiety at the disruption of the ecological balance in the Himalayan mountains and its deleterious effects on the environment and the economy of the adjacent Terai plains [1].The last National Forest Inventory (NFI) carried out in the late 1990s reported that, of the five physiographic regions, the middle mountains witnessed the highest annual rate of deforestation (2.3%) between 1978 and 1994 [2].Deforestation and forest conversion in the Nepalese mountains is not a recent phenomenon but have a long history, caused mainly by the joint attack of government land-use policy and subsistence agriculture [3].In last decades, however, major driving factors of deforestation cannot be attributed merely to the land-use policy and agriculture land expansion.Various socioeconomic, governance and institutional factors interact to cause deforestation in the mountains of Nepal (for more detail on drivers of deforestation in the mountains of Nepal, see the Ministry of Forests and Soil Conservation [4]).Past efforts to control deforestation and degradation have failed mainly because many driving factors of deforestation operate outside the forest sector [5].Although many of these driving factors persist, several successful efforts in recent years have led to a reversal of the ongoing trends with some remarkable success, notably the one achieved through community-based forest management.Previous studies have pointed to the success of Nepal's community management in restoring forest cover in many locations of the middle mountains [6][7][8][9].In other areas, forests were found being depleted and degraded due to deforestation, agricultural expansion, excessive grazing and expansion of rural road networks [1,10].
Previous studies have highlighted the need for research on location-specific land-use and land-cover (LULC) dynamics for decision-making processes related to the management of forest resources in the mountains of Nepal [3,7,[10][11][12].Remote sensing based methods make it possible to study LULC dynamics in less time, at low-cost and with better accuracy over large geographic areas.Its integration with the geographic information system (GIS) provides a suitable platform for data analysis, update and retrieval in a regular way [13][14][15].Over the past few decades, there have been great improvements in available remote sensing data, and a parallel advance in GIS that has allowed processing a large quantity of remote sensing data [13].Remote sensing data have been widely used in generating and updating LULC maps for several decades now, and spatial-temporal change detection has become one of their major applications [16].Particularly in areas of rugged topography and poor accessibility, remote sensing combined with GIS serve as a valuable tool for monitoring spatial-temporal LULC dynamics and their impacts [17].
The conventional method of assessing LULC change involves a transition matrix obtained from the comparison of bi-temporal maps.It is a common practice to derive the net change by LULC category from the transition matrix [18][19][20].However, net change can be substantially less than the total change because net change fails to account for simultaneous loss and gain of a category, which is a phenomenon known as swap [21,22].Therefore, assessing the traditional transition matrix beyond the size of each transition can reveal information on possible processes that determine the patterns in a landscape.For instance, swap changes reveal the amount swapped for each land-cover category, and persistence offers additional information concerning the vulnerability of the land-cover category to transition to other categories [22,23].In recent years, major landscape changes observed in the mountains of Nepal are (i) decrease of forest cover in the valley floors and foot-slopes due to deforestation and conversion processes; and (ii) increase of forest cover in the upper slopes due to participatory conservation and marginal agriculture land abandonment [24].This process of forest redistribution involves substantial ecological changes, which cannot be explained by a simple analysis of net changes [25].Therefore, in addition to swap changes and persistence, assessing dominant signals of landscape changes require reporting temporal trajectories of land-cover changes that would help link spatial patterns to change processes [18].The trajectories, defined as successions of land-cover types for a given area over more than two observation years, take widely different forms and allows a better projection of regions with a high probability of land-cover change, e.g., deforestation, than based on observations from only one previous period [18].
Watersheds in the middle mountains of Nepal have undergone notable LULC changes over the past few decades [10,11,17,[26][27][28].They have witnessed substantial population and economic growth, which has led to exacerbation of point and non-point source pollution, deforestation, hydrological alterations, sedimentation and direct habitat destruction.Located in western Nepal, Phewa Lake watershed was a notable example of degradation in the middle mountain region, where forest lands were overgrazed and degraded and stripped of trees for fodder and firewood [28].Since the late 1980s, there has been a shift in focus of conservation towards participatory and community-based conservation and management [9].This has led to mixed success in terms of forest coverage with forest recovery in some parts, while continued deforestation and forest degradation in other parts of the watershed [29].The results of previous studies suggest that the spatial distribution of forest cover change is not homogeneous throughout the landscape and that the spatial patterns of change are conditioned by complex interactions between environmental and socioeconomic factors.However, many previous studies in the middle mountains have ignored complex sequences of forest cover changes as they simply measured the conversions from one land-use/land-cover category to another between only two time points.To that end, our study used satellite image processing of Landsat images from 1995, 2005 and 2017, and GIS techniques to analyze LULC changes, especially those related to forest cover changes, in Phewa Lake watershed of western Nepal.The specific objectives of the study were: (i) to detect and quantify spatial-temporal trends of land-use/cover changes in general and forest cover, in particular between 1995 and 2017, with special focus on analysis of the level of persistence and swaps of forest; and (ii) to analyze land-cover change trajectories associated with forest cover changes (e.g., deforestation) in the study area landscape for the given time frame.The trends of land-use and forest cover changes in Phewa Lake watershed can be representative of many middle mountain watersheds of Nepal that harbor substantial population and economic growth, and an increasing pace of urban expansion.More importantly, a study of spatial-temporal changes of land-use and forest cover is imperative for conservation and management of Phewa Lake watershed because the watershed provides habitat for several species of the locally endangered flora and fauna [30], provides livelihood options to tens of thousands of households in both upstream and downstream areas [27], and, together with Pokhara city, serves as the backbone of the regional economy of the western Nepal.

Study Area
Located within 28 • 11 39" N to 28 • 17 25" N latitude and 83 • 47 51" E to 83 • 59 17" E longitude in western Nepal, Phewa Lake watershed (Figure 1) forms a unique geographical entity and represents the typical characteristics of a mountain watershed of Nepal.Administratively, the watershed extends over the whole or parts of six Village Development Committees and seven wards of the southwestern part of Pokhara Sub-Metropolitan City.It covers an area of 120 km 2 with the geometrical east-west length of 18.32 km and north-south width of 9.53 km.Complex and rugged ridges and spurs and valley bottoms stand out as dominant landforms of the watershed [31].Topography is steep with slope averaging 40%.The outlet of the watershed is situated at 763 m above sea level, whereas Panchase ridge summit (2471 m above the sea level) marks its highest point.
The watershed is composed of five sub-watersheds [32] that are drained by 19 streams and small brooks into Phewa Lake, the second largest and one of most prominent lakes of Nepal [28].The watershed climate is based on annual monsoon events, which bring more than two-thirds of the yearly rainfall between June and September [33].Average yearly precipitation ranges from 3811 mm on the valley floor to about 4700 mm on the ridges.Primarily due to elevation differences, there is substantial spatial variation in the temperature range [10].The annual mean temperature on the valley floor is 21 • C, with monthly means ranging from 13 • C in January to 26 • C in July.The temperature decreases gradually from the valley floor to the ridge.The temperature in the ridge average 16 • C, with monthly means ranging from 9 • C in January to 20 • C in July [27].Extreme natural events contribute to the natural degradation of the steep terrain including landslides and flash flooding [33].There are three distinct micro agro-ecological zones in the watershed [27].The narrow river valleys have a sub-tropical climate where lands are being utilized mainly for irrigated rice and maize cultivation.Hill slopes, with a sub-tropical climate in the lower elevations and temperate climate in the higher elevations, are sandwiched between valley floors and the ridges.Rain-fed rice is a dominant crop in the lower elevations, whereas mixed cultivation of maize, millet, vegetables, wheat, beans, potatoes and fruit trees is a common practice in the higher elevations.Ridges, characterized by a temperate climate, are utilized mainly for millet and maize cultivation [27].
A study by Sharma et al. [30] identified vegetation (49%) as the dominant land cover in the watershed followed by agriculture (41%), water bodies and swampland (5%), built-up areas (3%) and sand (1%).The vegetation is dominated by broad-leaved forests constituting 98% of the total forest area in the watershed.Major vegetation types found include hill Sal (Shorea robusta) forest, Katus (Castanopsis indica) and Utis (Alnus nepalensis) forest, Khasru (Quercus species) forest and Gurans (Rhododendron species) forest.Participatory watershed conservation efforts in the 1980s led to the devolution of forest management with over 60% of the forested land handed over to the communities as community forest [28].As of 2016, seventy-five community forestry user groups (CFUGs) manage over 2700 hectares of forest [9].Decades of participatory watershed conservation and community-based forest management has led to a restoration of degraded forests in many parts of the watershed [28].Agriculture and built-up land occupy most of the flat and gently sloped areas, while forest cover dominates the upper slopes and ridges [32].

Data Used
For LULC mapping and spatial-temporal change analysis, time-series of Level-1 Precision and Terrain (L1TP) corrected Landsat data over a period of 22 years (1995-2017) were obtained from the EarthExplorer database (https://earthexplorer.usgs.gov/) of the United States Geological Survey (USGS).The dataset includes two Landsat 5 Thematic Mapper (TM) scenes (185 km × 185 km) of November 1995 and October 2005, and a Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) scene (185 km × 185 km) of October 2017 (Table 1).The study area is entirely contained within Landsat path 142, row 140.The data from the Landsat TM sensor have seven spectral bands with a spatial resolution of 30 m. Data from the Landsat 8 have nine spectral bands with a spatial resolution of 30 m in bands 1-7 and 9, while its eighth band (panchromatic) has a resolution of 15 m.All visible and infrared bands (except the thermal infrared) were included in the analysis.Acquisition dates were chosen between October and November, a transition period between the rainy season and the dry season during which classes of interest can be distinguished with minor confusion.After assessing visually for contamination by in-scene components such as clouds and haze, maximum contaminant free images were selected.Autumn (October-November) is the appropriate season for acquiring cloud and haze free satellite images for Nepal [34].Landsat satellite images are widely used to map and monitor land-use and land-cover changes [35].
In addition to Landsat data, various other spatial data sources were used to facilitate land-use and land-cover classification and change detection.Four topographic maps (sheet number: 2883-12c 12d, 16a, 16b) of 1:25,000 scale with contour intervals of 20 m published in 1995 by the Survey Department of Nepal, and digital topographic data produced by the same agency were utilized.Topographic sheets were scanned and saved in Tiff format.The digital topographic data contained information on land-use, administrative boundaries, roads, and hydrography.

Satellite Data Pre-Processing
Preprocessing of satellite images prior to actual change detection is essential and has unique goals of establishing a more direct linkage between the data and biophysical phenomena by removing data acquisition errors and image noise [36] in order to produce a corrected image that is as close as possible to the radiometric and geometric characteristics of the original scene.Of the various requirements of preprocessing for change detection, radiometric, atmospheric and geometric corrections are the most important [37].Radiometric correction routines applied in this study were at sensor radiance and Top of Atmosphere (ToA) reflectance conversion using Band Math function in ENVI 5.3 (Exelis Visual Information Solutions, Boulder, CO, USA) following Chander et al. [38].Values of coefficients required to perform the radiometric correction including band specific rescaling factors, Earth-Sun distance, and solar zenith angle were obtained from the corresponding Landsat header files.Exoatmospheric solar irradiance values were obtained from Chander et al. [38].Dark Object Subtraction (DOS) [39] was then employed to convert TOA reflectance to surface reflectance using Dark Subtract with a band minimum function in ENVI 5.3.As Landsat data used in this study were precision registered and orthorectified to UTM coordinates (Zone: 45N, Datum: WGS 84) within <12 m radial root mean square error by USGS through a systematic process that involves ground control points and a DEM [40], no further geometric processing were required for this study.To ensure consistency, all Landsat data were re-sampled to a common nominal spatial grid of 30 m using the nearest neighbor resampling algorithm to maintain the spectral integrity of the data [41].

Creation of Subset (Study Area) Map
Traditionally, watershed boundaries are drawn manually onto a topographic map using topographic features on the map to locate the watershed divide [42].With the availability of DEM and GIS tools, watershed properties (such as area, slope, flow length, etc.) can be extracted using automated procedures.Phewa Lake watershed boundary for this study was delineated and extracted using source Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) and several functions available in the Spatial Analyst toolbox of ArcMap 10.2 (Earth System Research Institute, Redlands, CA, USA) following Bajjali [43].The individual preprocessed Landsat bands for each year were imported and layer stacked into a single file in ENVI 5.3.Finally, the layer stack was clipped using an area of interest (AOI) file derived from the vector dataset of the watershed boundary and exported as ERDAS Imagine (.img) file format.

Image Transformation
Spectral indices are among the most common data transformations for mapping land-use and land-cover.The principal advantages of spectral indices over single-band radiometric responses are their ability to considerably reduce the data volume for processing and analysis, and their inherent capability to provide information not available in any single band [36].Moreover, spectral indices are potentially less sensitive than single band values to artifacts due to differences in light conditions, exposed soil or terrain slope and orientation.One widely used spectral index is the Normalized Difference Vegetation Index (NDVI) [44] which has been used to monitor vegetation conditions on regional and global scale [13].It is used to distinguish healthy vegetation from others or from non-vegetated areas using red and near-infrared reflectance values and this was integrated in object-based image analysis to discriminate between vegetated (forest) and non-vegetated areas.Mathematical formula for calculating NDVI is where ρ NIR is surface reflectance value for the near-infrared band, and ρ RED is surface reflectance value for red (visible) band.Theoretically, NDVI value ranges from −1 to +1.Areas of water, barren rock, sand, or snow usually show very low NDVI values, while sparse vegetation such as shrubs, and grasslands show moderate NDVI values (approximately 0.2-0.5).High NDVI values (approximately 0.7-0.9)correspond to dense vegetation.For improved classification, other spectral indices, Ratio Vegetation Index (RVI) [45], Soil Adjusted Vegetation Index (SAVI) [46], and Normalized Difference Water Index (NDWI) [47], were calculated as: where ρ GREEN is surface reflectance value for the green band and L is the soil brightness correction factor and its value range between 0 for very high vegetation cover and 1 for low vegetation cover [46].L value of 0.5 which represents intermediate vegetation cover is used in most research studies.

Classification Scheme and Class Description
A modified Anderson Level I scheme [48] was used in this study.It utilizes a hierarchical LULC classification with classes at Level I able to be mapped from medium resolution satellite imagery such as Landsat 5 TM, whereas the extraction of information for level II, III and IV require the use of high-, medium-and low-altitude aerial photographs, respectively [49].Accordingly, four information classes were derived: agriculture and built-up, forest, waterbodies and other (Table 2).The class other consisted of grassland, wetland and barren and degraded land that individually covered a small part of the study area.The choice of the classification scheme was guided by the objective of the research, the easiness of identifying each class on the Landsat images, and expected degree of accuracy during image classification.

LULC Type Description
Agriculture and built-up Includes planted and bare agriculture land, terraced or valley; orchards, residential and commercial buildings, roads, construction sites, temple, and school/college Forest All forest types found in the study area; includes broadleaved mixed hardwood forest at lower elevations and coniferous forest at higher elevations [17] Waterbodies Includes all water bodies (lake, pond, stream, brook, and reservoir) Other Includes grasslands, wetlands, and barren lands (including sand and degraded land)

Object-Based Image Classification
eCognition Developer 9.0 (Trimble Inc., Sunnyvale, CA, USA) was used to perform object-based image classification.Object-based image classification incorporates two main steps: segmentation, which defines the image objects, and the classification itself [50].Image segmentation is the principal step that splits an image into separated regions or objects, also called "Image Object Primitives", based on parameters specified [51].A sophisticated segmentation algorithm known as "Multi-resolution Segmentation" which is based on the Fractal Net Evolution Approach (FNEA) [52] was utilized for image segmentation in eCognition.The algorithm minimizes the average heterogeneity for a given number of objects and maximizes their homogeneity based on defined parameters.eCognition provides three key parameters, namely shape (Sh), compactness (Cm), and scale (Sc), to segment objects or pixels having similar spectral and spatial properties [50].These parameters are generally defined through trial and error.As objects of interest exist within a range of different sizes and shapes, no single spatial resolution is sufficient to capture all of their spatial characteristics.Conveniently, eCognition allows segmentation at multiple levels with different scale parameters leading to a formation of a hierarchical network of objects [53].The larger number of scale factor (e.g., 100) generates large homogeneous objects (smaller scale), whereas the smaller number of scale factor (e.g., 10) generates smaller and less homogeneous objects (larger scale).The segmentation process stops when the most modest growth exceeds the threshold defined by the scale parameter.
Initially, all the pre-processed bands including transformed bands for each reference year were uploaded in eCognition Developer 9.0.After testing for various scale levels and parameter values, scale factors from 10 to 50 were considered to be appropriate for the study.A hierarchical scheme of three levels with scale factor values of 10, 25 and 50 was finally chosen to facilitate the object-based classification depending on the spatial characteristics of LULC classes.For instance, level 1 was to handle large size landscape features such as lakes, whereas level 3 was for small size features such as ponds.The shape parameter was set to 0.1 while the compactness parameter was set to 0.5 for all three levels.Near-infrared (NIR) and NDVI bands were assigned a weight of two while all the remaining bands were assigned a weight of one.Level 1 segmentation was carried out at the pixel level, while Level 2 and Level 3 segmentations were carried out at image object levels, the output of which is shown in Figure 2. Classification of the image segments was the sequel to multi-resolution segmentation.There are two approaches in eCognition to assign classes to segmented objects: membership function classifier, and the nearest neighbor (NN) classifier [50].The membership function classifier uses a user's expert knowledge to define rules and constraints in the membership function to control the classification procedure.On the other hand, the nearest neighbor classifier uses a defined feature space, for example, using original bands or customized bands and a set of samples that represent different classes, to assign class values to segmented objects.Application of the nearest neighbor method is advantageous when using spectrally similar classes that are not well separated using just one feature or a few features [50].Both approaches were used in this study to classify the image segments.The membership function classifier was employed for classification at Level 1 (scale factor 50) and Level 2 (scale factor 25), while NN classifier was employed at Level 3 (scale factor 10). Rules were developed for each TM and OLI scene to overcome the inherent variability that arises in reflectance values across the scenes because of the different times of acquisition.Training samples for NN classifier were collected using the topographic maps (Landsat 5 TM 1995), and high-resolution Google Earth images (Landsat 5 TM 2005 and Landsat 8 OLI 2017).At the end, sub-classes generated from the hierarchical classification were assigned to their respective LULC classes defined in the class hierarchy using Assign Class function available in eCognition.

Classification Accuracy Assessment
The accuracy assessment in this study was performed through a standard procedure described by Congalton and Green [54].Google Earth images from December 2005 and October 2017 were used as the reference for the classified maps of 2005 and 2017 respectively, and a 1:25,000 scale topographic map as the reference for the classified map of 1995.Google Earth-based high-resolution images and topographic maps of appropriate scale are widely adopted for accuracy assessment of classified land-use and land-cover products [55][56][57].A disproportionate stratified random sampling technique was employed to generate reference pixels in ENVI 5.3.The reference pixels were then exported to ArcMap 10.2 as ESRI Shapefile (.shp) where each pixel was converted to point and XY coordinates assigned for subsequent comparison against the reference data.A total of 249 reference points with at least 50 points per class (the number of points increasing with the size of land-cover classes) were generated independently for the three classified images.It has been suggested that a minimum of 50 reference points for each land-use/land-cover class be collected for a successful classification accuracy assessment [58].
The reference points were then used to produce an error matrix.Since the sampling intensities are different in each row of the error matrix, proportion correct are likely to be biased when computed directly from the entries in the matrix [59].Therefore, sample error matrix of each year was converted to population error matrix following Pontius Jr. and Millones [59] to compute unbiased summary statistics of classification accuracy including omission intensity and commission intensity for individual LULC types, and overall error.Producer's and user's accuracy for individual LULC types were computed from the population error matrix following the procedure of Congalton and Green [54], which were then used to derive omission intensity (100-producer's accuracy) and commission intensity (100-user's accuracy) statistics.Overall error was calculated as the sum of the off-diagonal values of the population error matrix multiplied by 100.If commission intensity is higher than the omission intensity for a particular LULC type, then the classified map overestimates the quantity of that LULC type.Conversely, if omission intensity is higher than the commission intensity for a particular LULC type, then the classified map underestimates the quantity of that LULC type [60].

Post-Classification Comparison and Change Analysis
A pixel-based change detection called post-classification comparison (PCC), where images from different dates are classified independently and then compared pixel-by-pixel using map overlays [61], was selected to determine changes in land-use/land-cover during periods 1995-2005, 2005-2017 and 1995-2017.PCC is the most commonly used change detection method for comparing data from different sources and dates [62].Classified maps of all three time points were compared pixel-by-pixel using Thematic Change Workflow module available in ENVI 5.3.The output table from the comparison was processed in Microsoft Excel to produce a change matrix.From the matrix, a summary table of overall changes per class was created and values were presented in terms of hectares and percentages.
Patterns of forest cover change were summarized in terms of persistence, gain, loss, swap and net change from the change matrix following Pontius Jr. et al. [21].Persistence measures the area of a LULC category that did not change spatially during a period.On the other hand, gain measures the amount of area gained by a particular category during a specified period, and loss measures amount of area lost the category during the same period.Total change for a category is the sum of gain and loss for the category.On the contrary, total absolute net change for a category is the absolute value of the category's gain minus its loss.Lastly, swap refers to the loss amount by a land-use/land-cover class at a location that is accompanied by the same amount of gain at another location and is computed by subtracting net change from the total change [21].Annual change is derived by dividing the periodic change by the number of years between the first and the second time point of the period.Intensity analysis at the interval and the category level was carried out following Aldwaik and Pontius Jr. [63,64] to determine whether the forest category experienced change with the same annual intensity (fast vs. slow) and the same dynamics (active vs. dormant) at different time intervals.Active means that the intensity of a category's annual loss or annual gain is greater than the intensity of annual overall change during a particular time interval.Dormant means that the intensity of a category's annual loss or annual gain is less than the intensity of annual overall change during a particular time interval [63,64].In addition, we assessed gain-to-persistence ratio (Gp), loss-to-persistence ratio (Lp), and net change-to-persistence ratio (Np) indices proposed by [23] to determine the vulnerability of the forest cover change.The annual rate of forest cover change was calculated based on the compound-interest-rate formula proposed by Puyravaud [65]: where P is the percentage of forest cover change per year and A 1 and A 2 are the amount of forest cover at time t 1 and t 2 , respectively.

Trajectories of Forest Cover Change
To better understand general trends of land-cover during two consecutive periods, an analysis of land-cover trajectories concerning forest cover changes was carried out.The concept and methodology of land cover change trajectory have been developed already [18,66,67].Here, trajectories of forest cover change refer to the successive transitions between forest and non-forest cover categories over the observation years.Following Mertens and Lambin [18], all possible land-cover trajectories (T = n x , where n is the number of land cover classes and x is the number of observation years) were reduced to eight main sequences based upon the forest and non-forest categories (Table 3).The description of the trajectories in this study is based on the analysis of classified maps with only two binary aggregate classes (forest and non-forest) from the three observation years, but a particular transition (for instance, "forest" to "non-forest") could have occurred anytime between the observation years.For estimation of the trajectories, the classified maps of 1995, 2005 and 2017 with the aggregate forest and non-forest classes were first vectorized and imported into ArcMap 10.2 as ESRI Shapefile (.shp) where they were successively intersected starting with the map of 1995.A trajectory map was created to visualize spatial trends of forest cover change.

Land-Use and Land-Cover Classification and Accuracy
Multi-temporal Landsat images of 1995, 2005 and 2017 were independently classified and classification maps generated (Figure 3). Figure 3a-c depicts the spatial distribution of various land-use and land-cover types in Phewa Lake watershed for years 1990, 2005 and 2017, respectively.The results of LULC distribution for the three time points show that forest and agriculture and built-up were the dominant land-use/land-cover type in the watershed (Table 4).Forest cover was more prevalent on the southern part of the watershed and the upper slopes elsewhere, whereas agriculture and built-up was more prevalent on valley floors, hill terraces, and foot-slopes.Phewa Lake, located in the eastern part of the watershed, was responsible for most of the water cover in the area.Spatial distribution of type other (barren land, wetland, and grassland) shows that it is distributed mainly in the vicinity and along the waterbodies.Overall, over a duration of 22 years, forest cover increased while the area of all other land-use/land-cover types decreased.Omission intensity and commission intensity for each LULC type and the overall error were computed from the error matrix (Table S1) and summarized for each year in Table 5.The classification maps have overall classification errors of 7%, 10%, and 8%, respectively, for 1995, 2005 and 2017.The omission intensity for individual LULC type ranged between 2% and 47% for the classified map of 1995, between 7% and 81% for 2005 and between 6% and 64% for 2017.On the other hand, commission intensity ranged between 2% and 10% for 1995, 5% and 15% for 2005 and 4% and 11% for 2017.All classified maps consistently underestimated waterbodies and other.In particular, the omission intensity of LULC type other was high for all years (i.e., low producer's accuracy), implying a higher probability (proportionate to the errors) that ground reference points for this category were classified incorrectly.

LULC Change
The area statistics on various LULC types of Phewa Lake watershed are summarized for each year (Table 4).Forest, which is the primary vegetation type, occupied 48.7% of the watershed area in 1995, slightly decreasing to 48.6% in 2005 before increasing to 52.1% in 2017.On the other hand, the area under agriculture and built-up increased to 46% of the watershed in 2005 from 44.8% in 1995 before decreasing back to 42% in 2017.The area under waterbodies progressively reduced, from 4.09% (489 ha) in 1995 to 3.97% (474 ha) in 2005 to a further 3.81% (455 ha) in 2017.Other type comprising the area under wetland, barren and degraded land and grassland is the least prevalent LULC type in Phewa Lake watershed and experienced changes in two phases: it decreased from 2.32% (278 ha) in 1995 to 1.30% (154 ha) in 2005 and then increased back to 2.12% (253 ha) in 2017.
Forest cover witnessed a net annual change of 0.31% between 1995 and 2017 (Table 6).During the same period, agriculture and built-up use changed by −0.29% annually.The area under waterbodies and other also experienced a net annual change of −0.31% and −0.41%, respectively between 1995 and 2017.Percentage difference between maps of 1995 and 2005, 2005 and 2017, and 1995 and 2017 were observed at 16.24%, 17.60%, and 18.84%, respectively.Overall errors (Table 5) derived from the error matrix of 1995, 2005 and 2017 were smaller than the observed map differences between the time points.

Forest Cover Change
Patterns of forest cover were summarized in terms of net change, gross gains and losses with a special focus on persistence and swaps (Table 7).An investigation of changes in forest cover of the watershed revealed that of the total 5824 ha of forest 1995, approximately 86% persisted during 1995-2005, and 87% persisted during 2005-2017.Between 2005 and 2017, 5107 ha (87.84%) forest cover persisted.The remaining forest cover was converted to non-forest uses.Figure 4 shows the gains, losses, and persistence of forest cover in the study landscape between 1995 and 2017.Between two periods 1995-2005 and 2005-2017, the latter period observed a more substantial gross gain in forest cover of 1110 ha, while the former witnessed a gross gain of 813 ha.The gross gains of forest category are accompanied by gross losses in another category, so the total gross gain is equivalent to the total gross loss in the landscape.In our study, these gains, for the most part, came at the cost of agriculture and built-up use.On the contrary, large amount of forest cover was also lost to remaining LULC types, primarily to agriculture and built-up during both periods.The period 1995-2005 witnessed a greater forest loss compared to the 2005-2017 (Table 7).Consequently, forest cover decreased by 9 ha in the first period (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) and increased by 403.29 ha in the second period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017).Results of gain and loss intensity show that forest cover gained actively during both time intervals (Table 7).In terms of losses, forest cover lost actively during the first period from 1995 to 2005, and then dormantly during the second period between 2005 and 2017.Regarding the extent, annual forest gain was highest during the second period, while the annual forest loss was highest during the first period.Concerning annual change, 164 ha of forest cover changed annually during the first period, whereas, in the second period, annual forest cover change lowered to 151 ha.
Given the dynamics, simply observing and reporting the value of net changes can lead to a biased conclusion about total landscape changes.In our study landscape, this is evident from a small net change value of forest cover between 1995 and 2005, but over 175 folds larger value of swap changes.Forest cover for both periods showed a Gp and Lp values less than one indicating that both gains and losses during the periods are lower than the amount that persisted.Neither the gains nor the losses in forest cover exceeded 25% of the area that persisted.The net change to persistence (Np) was negative, approximately 0.002, for 1995-2005, and positive, approximately 0.08, for 2005-2017 implying a net decrease of 0.2% of the persistence in the first period while the net increase in the second period was 8% of the persistence.Table 7 also reveals that the annual rate of forest cover change between 1995 and 2005 was −0.02%.However, between 2005 and 2017, the trend of forest cover change in Phewa Lake watershed reversed.During this period, the annual rate of forest cover change was 0.56%.

Trajectories of Forest Cover Change
Forest cover change trajectories were computed using binary recoded LULC classes: forest and non-forest.Eight trajectory classes were observed for 1995, 2005 and 2017 (Figure 5).The change trajectory analysis reveals that, over a 22-year period with three observation years, 4622 ha of the initial forest cover remained unchanged (Table 8).The area that remained unchanged for non-forest class was similar.Together, the stable trajectories, i.e., stable forest and stable non-forest accounted for 77.5% of the total landscape area.The remaining 22.5% of the landscape changed at least once during the observation period.Between 2005 and 2017, 773 ha (6.4%) of existing forest was removed (sequence 2 and sequence 3).On the other hand, 486 ha (4.06%) of old and permanent forest regrowth was recorded between 1995 and 2007 (sequence 6).In addition, between 2005 and 2017, 1110 ha (9.3%) of forest was added to the landscape.

Discussion
Land cover samples were taken from calibrated and atmospherically corrected images and then plotted against their wavelength (µm).Six spectral bands were plotted for Landsat 5 TM 1995 and 2005, and seven spectral bands were plotted for Landsat 8 OLI 2017 (Figure S1).An upward trend of the blue band (≈0.48 µm) as observed in Figure S1a-c could be considered a prominent atmospheric effect, likely caused by atmospheric aerosol scattering.Atmospheric correction applied compensated for the aerosol scattering to produce spectra that closely represent the surface reflectance.The derived reflectance curves of agriculture and built-up, forest, water and other LULC types were, therefore, closer to their typical curves (Figure S1d-f).Forests showed characteristic absorption features in blue and red band, a slight peak in the green band and a sharp edge leading to higher near-infrared reflection.Compared to forests, agriculture and built-up showed higher reflectance in the visible region, but a lower reflectance in the near-infrared region.Waterbodies showed a good absorption in the short wave infrared region, distinguishing itself from other LULC types.Classification accuracies indicated that the object-based image classification employed in this study performed well for all LULC types with the exception of other type.The omission intensity of the LULC type other, which includes wetland and grass and barren land, was consistently high (i.e., low producer's accuracy) for all year.This could be a result of spectral mixing with agriculture and built-up land.In a similar study, Myint et al. [68] reported that spectral confusion between barren land and agriculture and urban lands resulted in a lower categorical accuracy of the former land-use type.
More than 18% of the watershed landscape experienced changes during the last 22 years.LULC distribution showed that forest and agriculture and built-up were the dominant LULC type in the watershed.The general trend observed between 1995 and 2017 was an increase of forest cover while a decrease in the extent of remaining LULC types that included agriculture and built-up, water and other.Over a period of 22 years, 15 ha of land under agriculture and built-up use had been converted annually to other LULC types, mainly due to the abandonment of marginal agriculture lands [9,17].Recently, an increasing trend towards de-intensification and abandonment of agricultural land has been reported in many parts of the Nepalese middle mountains [69][70][71] due to increased outmigration and a shortage of labor force working in the agriculture sector.Waterbodies (dominantly Phewa Lake) decreased by 6.87% from 1995 to 2017.This decline is a clear result of shoreline encroachment for agriculture and commercial purposes [31].Annually, the area under waterbodies decreased by two hectares between 1995 and 2017.
Total change in forest cover observed between 1995 and 2017 was 1923 ha, of which approximately 60% were gains from other LULC types.The spatial pattern of forest cover changes revealed that until the early 2000s, forest loss occurred primarily in the western part of the watershed where logging for fuelwood and heavy cattle grazing is widespread [28].After the restoration of peace in the mid-2000s, the eastern part witnessed substantial forest loss due to rapid urbanization, increased population pressure and tourism infrastructure development.Restoration of peace also brought increased forest cover, however, in the western part of the watershed as a result of forest growth in marginal agricultural land abandoned by rural farmers who left their homestead in search for better livelihood opportunities in nearby cities and Gulf countries.Of the total forest cover change, over three-fourths (about 80%) could be attributed to swap changes.Swap changes between 1995 and 2005, and 2005 and 2017 accounted for 99.43% and 77.80% of the total forest cover changes in the study landscape, respectively.Higher values of forest swap indicate features of ecological damage and ecological restoration at the same time but different locations.In the river valleys and foot-slopes suitable for paddy cultivation and terrace agriculture, local people have destroyed the forest for agriculture, grazing, and fuelwood.In upper slopes, however, local people have abandoned marginal agricultural lands which saw a succession of shrubs and other inferior woody vegetation.Regmi et al. [29] reported that the areal extent of bush/shrub and open forest in the study landscape has increased dramatically in recent years at the cost of dense forests.Consequently, while forest cover as a whole is increasing, high swap change values may indicate reduced ecological and ecosystem benefits from the forests [72].Our findings highlight that both net and swap changes are essential variables to consider to understand total landscape changes, which is in agreement with Pontius Jr. et al. [21], Manandhar et al. [22], and Barimoh et al. [23].
The annual rate of forest cover change measured for the first period (1995-2005) was −0.02%.However, considering forest clearing only and ignoring the forest regrowth that took place during 1995-2005, the annual forest cover change was computed at −1.4%-a value higher than the change rate presented in Table 7.The latter rate provides a better indication of potential impacts of deforestation on ecology and biodiversity, whereas the former rate is more relevant to establish a carbon budget or to monitor the availability of wood resources [18].Nevertheless, these values are still small compared to the countrywide deforestation rate of 1.7% or an even higher middle mountain deforestation rate of 2.3% [2].Overall, between 1995 and 2017, forest cover increased at an annual rate of 0.3% suggesting forest re-growth in previously deforested areas.These results corroborate the with findings of Gautam et al. [7], Gautam et al. [8], and Fleming and Fleming [28] that the forest cover in the middle hills of Nepal has increased and that deforestation has slowed or even reversed during the past few decades due, in large part, to the people's involvement in forest management and forest restoration.While these findings are encouraging, the deforestation rate varies over time and is strongly controlled by exogenous forces of demography and socio-economy that also change over time [18].Focusing exclusively on the rate of net forest change would overshadow many critical land-cover changes and trade-offs happening on the landscape.
The concept and methodology of change trajectory analysis were extended to the study of forest cover changes.The spatial pattern of permanent forest clearing showed that most of these clearings were made near agricultural and settlement lands.Local people clear trees, primarily for fuelwood and for construction materials to build a new house and livestock and poultry shed.In fact, for a country where 75% of the total energy demand is fulfilled by fuelwood [73], high dependency on forest resources for domestic use has been cited as one of the primary drivers of deforestation and forest degradation [5,74].In addition, rural road networks expansion and infrastructure projects (e.g., dams) have led to the permanent removal of forest cover in the hills of Nepal [75].Phewa Lake watershed witnessed a considerable increase in rural earthen roads in the past four decades, from about 23 km in 1979 to 310 km in 2013 [33], at the cost of hill forests spread across the watershed.Concurrently, recent forest clearings occurred primarily in the hills adjoining Phewa Lake in the eastern and central part of the watershed.Sarangkot, a prominent hill located in the northeast part of the watershed, in particular, experienced a large amount of recent forest clearings owing to the extensive tourism infrastructure development after peace restoration, as well as migrant influx from rural villages in the watershed [76].The drivers of recent forest clearings in all other parts of the watershed include agriculture expansion and extension of rural earthen roads.Meanwhile, the spatial pattern of regrowth in the watershed shows their occurrence in the eastern and central part due, in large part, to current and past efforts towards reducing soil erosion and downstream sedimentation of Phewa Lake through participatory watershed and forest conservation [28].Pandey [77] reported that conservation of forests in upland areas of middle mountains of Nepal reduces soil erosion, sedimentation, and flooding in the downstream.Marginal agricultural land abandonment accompanied by afforestation, reforestation, and conservation of existing forests by community forest user groups, who control over 60% of the total forest land of the watershed [28], have led to the restoration of forest cover in the landscape [17,78].Similar trends were observed in other mountainous watersheds of Nepal [7,79,80].Overall, the spatial pattern of forest cover change in the study area followed a process of diffusion with expansion in some regions, while contraction in the other, and that changes have mainly occurred in transition zones between forested areas and non-forested areas as illustrated in Figures 4 and 5. Our findings concur with Mertens and Lambin [18] that taking into account the forest cover change trajectories allow for more reliable projections of areas at risk of deforestation compared to considering only previous observations.Pontius Jr. and Lippitt [81] noted that classification errors might cause maps to show more agreement than exists on the ground.To conclude that the differences in maps of two time points essentially denote land cover changes without adequately considering potential sources of errors is misleading [82].Fuller [82] demonstrated that to map changes less than or equal to 20% with 90% reliability, classification errors should be around 1%.In this study, we found the overall errors in the classified maps of 1995, 2005 and 2017 to be smaller than the map differences between the time points.However, a lower value of classification errors compared to the observed map differences cannot be a precondition to making reliable deductions about changes [82].Therefore, the map difference results reported in the study should be treated with caution.Due to the propagation of classification errors, the amount of change on the ground might be different from that obtained from post-classification overlay procedures.
Although studies on the spatial-temporal patterns and relationships between forest cover and local land-use practices will help to find areas where changes have occurred and predict future changes, they are not sufficient for explaining local driving factors of change so that informed, science-based decisions could be made.Only a detailed understanding of the local driving factors can help to reveal the processes of forest cover changes, including deforestation and forest regeneration in the Himalayan Mountains at different spatial scales.Future studies should, therefore, focus on why and how the changes in land-use and forest cover have occurred within the socioeconomic, policy and institutional and historical context of this landscape, and how those are interconnected with the broader changes occurring at the regional and global level.

Conclusions
Information about changing patterns of land-use and land-cover through time is essential, not only for planning and management, but also for a better understanding of the relationship between landscape dynamics and ecosystem responses.Our approach of combining multi-temporal remote sensing data and GIS techniques has quantified and characterized the spatial-temporal pattern of LULC changes, especially those related to forest cover changes, in a mountainous topography.The results showed that more than 18% of the watershed landscape experienced LULC changes between 1995 and 2017.Overall, Forest cover increased between 1995 and 2017, whereas land under agriculture and built-up, waterbodies and other use decreased.The annual rate of forest cover change was 0.3%.The spatial pattern of change revealed that the forest cover experienced gains and losses with different annual intensity and dynamics between 1995-2005 and 2005-2017.While forest cover loss was active during the first interval (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) and dormant during the second interval (2005-2017), the annual gain was active during both time intervals.The net change and total change in forest cover between 1995 and 2017 were 394 ha and 1923 ha, respectively.Swap change accounted for 79.51% of the total forest cover change during the same period.Temporal trajectories revealed that forest cover change in the study area followed a process of diffusion and that changes were mainly occurring in transition zones between forested areas and non-forested areas.Quantitative evidence of net forest cover changes and rates presented here agrees with the findings of some earlier studies carried out in the middle mountains of Nepal.However, high forest swap changes in the study area indicate that, while forest cover may be expanding at present, its continuation in the future is equivocal, especially with increasing human population and urban pressure on the existing forest cover.This is substantiated by non-linear forest cover changes uncovered by trajectory analysis against the popular belief of linear trends of forest cover changes in the Himalayan Mountains.Substantial areas subject to a forest conversion in one period were subject to another change in the following years.

Figure 5 .
Figure 5.The trajectory classes of forest cover changes in Phewa Lake watershed for three time points 1995, 2005 and 2017.
: Error matrix of classified vs reference data for: (a) 1995; (b) 2005; and (c) 2017, Figure S1: An example of spectral profile curves of agriculture and built-up, forest, waterbodies and other LULC types derived from: (a-c) the original digital number (DN) image; and (d-f) the atmospherically corrected and rescaled reflectance image of 1995, 2005 and 2017, respectively.Author Contributions: P.B. and Y.W. conceived and designed the paper; P.B. performed data collection and field verification; P.B. and N.N.U.analyzed data and wrote the paper; and Y.W. modified the paper.

Table 1 .
Details of Landsat data and the corresponding solar geometries.

Table 2 .
Land-use and land-cover classification scheme.

Table 3 .
Forest Cover Change Trajectories between 1995 and 2017.

Table 4 .
Summary of Landsat classification area statistics for 1995, 2005 and 2017 in hectares (ha) and proportion of the total landscape (%).

Table 5 .
Summary of Landsat classification errors for classified maps of 1995, 2005 and 2017.

Table 6 .
Change statistics (in ha and %) of LULC in Phewa Lake Watershed from 1995 to 2017.

Table 7 .
Summary of forest cover change between 1995 and 2017.
Values presented in hectares (bold) and percentage of area in Year 1 of each period (italics).For annual gain and loss, superscript a means active, and superscript d means dormant.* Percentage of swaps are expressed relative to the total change.

Table 8 .
Estimated area (ha) and proportion (%) of forest cover change trajectories in Phewa Lake watershed from 1995 to 2017.
F stands for "forest," and N stands for "non-forest."