Forest Cover and Sustainable Development in the Lumbini Province, Nepal: Past, Present and Future

: The analysis of forest cover change at different scales is an increasingly important research topic in environmental studies. Forest Landscape Restoration (FLR) is an integrated approach to manage and restore forests across various landscapes and environments. Such restoration helps to meet the targets of Sustainable Development Goal (SDG)–15, as outlined in the UN Environment’s sixth Global Outlook, which includes the sustainable management of forests, the control of desertiﬁcation, reducing degradation, biodiversity loss, and the conservation of mountain ecosystems. Here, we have used time series Landsat images from 1996 to 2016 to see how land use, and in particular forest cover, have changed between 1996 and 2016 in the Lumbini Province of Nepal. In addition, we simulated projections of land cover (LC) and forest cover change for the years 2026 and 2036 using a hybrid cellular automata Markov chain (CA–Markov) model. We found that the overall forest area increased by 199 km 2 (2.1%), from a 9491 km 2 (49.3%) area in 1996 to 9691 km 2 (50.3%) area in 2016. Our modeling suggests that forest area will increase by 81 km 2 (9691 to 9772 km 2 ) in 2026 and by 195 km 2 (9772 km 2 to 9966 km 2 ) in 2036. They are policy, planning, management factors and further strategies to aid forest regeneration. Clear legal frameworks and coherent policies are required to support sustainable forest management programs. This research may support the targets of the Sustainable Development Goals (SDG), the land degradation neutral world (LDN), and the UN decade 2021–2031 for ecosystem restoration.


Introduction
Forest ecosystems are vital sources of food, medicine, and fuel for human beings. Achieving sustainable forest management is an important step in the transition towards sustainable development as well as for providing economic, social, and environmental outputs [1]. Global forest coverage of the Earth's terrestrial surface declined from 31.6% to 30.6% between 1990 and 2015 [2], and every year, 10 million hectares (ha) of forest are lost [3]. To maintain a sustainable environment, regeneration and reforestation programs have become a top priority globally [4,5]. The Global 2030 Agenda for Sustainable Development Goals (SDG 15) was accepted by all countries [6]. The 15th SDG aims to "Protect, restore and promote the use of terrestrial ecosystems, environmental conservation, sustainably manage forest resource, combat desertification and halt and reverse land degradation and biodiversity loss" and has prioritized the monitoring of forests and the sustainable management of forest resources. It further endorses the guarantee of resilient agricultural practice, sustainable food production, and imaginative use of natural resources [7]. Many countries have collectively committed to restore 150 million ha of degraded land by 2020 to make a positive contribution in controlling forest degradation and protecting vulnerable landscapes [34]. Watershed management plans also contribute to landscape conservation and include slope management through terracing, trenching, and re/afforestation [35]. There are 118 ecosystem services have been recorded in Nepal [36]. The National Parks and Wildlife Conservation Act (1993), the fourth amendment, provisioned the investment of 30-50% of the total annual revenue of protected areas in local communities for biodiversity conservation and livelihood improvements [37]. Recently, local-level administrations have introduced urban plantation programs in urban areas, green road programs for major urban areas, and plantation programs in barren lands in rural areas of this region. All of these provide support for forest restoration, and we should expect forested area to continue to increase in the future. Private forest programs (PF) provide another option to increase the forested areas in Nepal [38] and are already in place in Lumbini Province. Against this backdrop, we attempted to extract land cover (LC) changes for the Lumbini Province of Nepal at decadal intervals from 1996 using remote sensing technology and to project future LC for the years 2026 and 2036 using an LC change simulation model.
Several studies have examined the spatiotemporal land use land cover (LULC) change of various parts of Nepal, such as the Sagarmatha National Park region during the period of 1992-2011 [39], in Kathmandu during the period of 1990-2010 [40], in the Tanahun district during the period of 1976-2015 [41], in Koshi river basin during the period of 1992-2010 [42], in the Bagmati river basin [43,44], and at the trans-boundary of the Gandaki river basin [45]. In this study, we have chosen Lumbini Province to complete the LC change analysis for several reasons. First, the forest in this region is the dominant LC in the Province, particularly in the northern belt, while the southern part of the study area comprises the rapidly urbanizing Tarai region, which is experiencing remarkable LULC changes. Since landscape changes and anthropogenic factors affect habitat quality and distribution [46,47], the monitoring of forest cover is crucial. Further, these socioeconomic and environmental factors have characterized the heterogeneous composition and overtime dynamics of forest cover. Third, to our knowledge, no other study has examined historical forest cover changes or has made predictions for future changes in the Lumbini Province of Nepal, and our research attempts to fill this gap. Understanding the dynamics of LULC and making reliable projections are imperative in resource governance and for improving land use [48] since forest resources and their management are essential for human well-being and for the maintenance of ecosystem services [49]. As such, the obtained research outputs will be the base for establishing the legal framework and for formulating coherent policies in support of the global targets of the Sustainable Development Goals (SDGs), the land degradation-neutral world (LDN), and the UN decade 2021-2031 for ecosystem restoration.
LC changes are a major driver of global environmental change and occur due to deforestation, reforestation, urban expansion, and the intensification of cultivated land [50,51]. LC change models can be used to predict the location and frequency of LC change [52]; to aid land use planning and conservation [53], urbanization [54,55], Lake area evaluation [56]; and to monitor environmental changes [57]. A number of simulation models are widely used, such as DINAMICA [53,58], SLEUTH [59], SERGoM [60], CLUE [61], GEOMOD [62], LUCAS [63] ANN-CA [64], and the CA-Markov change model [65][66][67][68]. Here, we evaluate past landscape change and predict future changes in forests and built-up areas. We use the Cellular Automata (CA)-Markov chain (MC) (CA-Markov) model, as these models appear to be excellent at predicting future LC changes and transitions [69] due to its strong hybrid functions [70].

Study Area
Lumbini Province is located in Western Nepal and is geographically located between 81. 10 (CBS, 2014). We have chosen Lumbini Province as the subject for the LC change analysis, as the Province has integrated national parks, hunting reserves, and conservation areas (Banke National Park, Bardiya National Park, the Dhor Patan hunting reserve, and the Krishnasar conservation areas) (Figure 1). It also includes Jagdishpur, a Ramsar-listed lake, and Lumbini, the birthplace of Lord Gautam Buddha-a place listed as a UNESCO World Heritage site in 1997. Additionally, Devdaha lake, the Anoma river, the Jitgadh fort, Tansen, and Ridi Ruru are also located within the Province.

Study Area
Lumbini Province is located in Western Nepal and is geographically located between 81.10-84.072 E to 27 323-28.833 N, covering about 19,256 km 2 . Elevation ranges from ± 90 m to 5600 m above sea level (masl). It includes a total of 12 administrative districts (six in Tarai and six in Hill region), covering 109 local administrative units. The total population of this Province was 3.45 million in 1991 and 4.96 million in 2011 (CBS, 2014). We have chosen Lumbini Province as the subject for the LC change analysis, as the Province has integrated national parks, hunting reserves, and conservation areas (Banke National Park, Bardiya National Park, the Dhor Patan hunting reserve, and the Krishnasar conservation areas) ( Figure 1). It also includes Jagdishpur, a Ramsar-listed lake, and Lumbini, the birthplace of Lord Gautam Buddha-a place listed as a UNESCO World Heritage site in 1997. Additionally, Devdaha lake, the Anoma river, the Jitgadh fort, Tansen, and Ridi Ruru are also located within the Province.

Data
In the study, we collected Landsat satellite Level 1 data (TM and OLI) for the years 1996, 2006, and 2016 from the USGS website (earthexplorer.usgs.gov) ( Table 1). All of the images were verified for image processing. Radiometric correction of remote sensing data mainly involves the correction of digital image errors [71,72], where image enhancement and registration are conducted for image correction [72]. The L1T (level one terrain corrected) satellite images were converted from DN (digital number) to radiance. The digital number was converted, and the FLASH model and was applied for image processing using ENVI software. ENVI software is organized with radiance calibration,

Data
In the study, we collected Landsat satellite Level 1 data (TM and OLI) for the years 1996, 2006, and 2016 from the USGS website (earthexplorer.usgs.gov) ( Table 1). All of the images were verified for image processing. Radiometric correction of remote sensing data mainly involves the correction of digital image errors [71,72], where image enhancement and registration are conducted for image correction [72]. The L1T (level one terrain corrected) satellite images were converted from DN (digital number) to radiance. The digital number was converted, and the FLASH model and was applied for image processing using ENVI software. ENVI software is organized with radiance calibration, geometric, and atmospheric correction for satellite images. A root mean square (RMS) error of the geometric rectification images that was less than 0.5 (<15 m) pixels was accepted. A 30-m shuttle radar topography mission (SRTM) digital elevation model (DEM) was used for the image registration for the years 1996, 2006 and 2016. The overlay function was conducted to fill the no data gaps in the ETM+ scenes where the TM data were unavailable (Table 1). There are multiple methods that can be used to fill in the gaps of the SLC-off Landsat images for years after 2003. The local histogram matching method (LHMM) was applied by Storey et al. in 2005 [73]. Similarly, the geostatistical interpolation method is another option that can be used to fill in the gaps of missing data [74]. A further option, a deep convolution neural network (CNN) model, can also be used to recover missing information resulting from the SL-off problem [75]. To fill in the gaps of the SLC-off images, auxiliary images were used to recover the missing data [76]. In our study, a small part of the covered area comprised SLC-off images for 2006, so we collected the topographical data developed by the Survey Department of Nepal for the scale of 1:25,000 [77]. Similarly, the LC data for the year 2000 [78] and further verified satellite images (Landsat TM) got 2008, November Nine (11-09) were used to verify the LC information of the missing data. Additionally, we collected Google images of the study area for the year 2006, and this provided the best information for the missing data for the year 2006. After all of the available auxiliary images were verified, the classified satellite images for 2006 were overlaid with the corrected missing data, and the final LULC map of 2006 was prepared.
Topographical data developed by the Survey Department of Nepal, 1996, and Google Earth images were used as the reference data for the image classification [79] of the entire study area. For the extraction of the LC data, we used the modified LC classification scheme recommended in Anderson et al. [80] (Table 2) and explored nine major LC classes: agriculture, forest, shrub, grassland, sand, barren land, water bodies, ice and snow cover, and other areas (including settlement road networks, industrial areas, infrastructure, and other planned areas). To extract the nine LC classes, all of the images for different tiles were stacked, subset, and analyzed using ENVI v5.3 software. We applied a non-parametric supervised support vector machine (SVM) to extract the LC for a specific area. A number of parametric and non-parametric algorithms [81] are frequently used for LC classification, such as minimum distance (MD), maximum likelihood (ML), support vector machines (SVM), artificial neural networks (ANN) and decision trees, and ML classifiers such as MD, BC, ANN, and fuzzy classification (FC) [82]. We applied supervised-learning SVM non-parametric and nonlinear approaches for the highly accurate extraction of LC change data [79,83,84].
Generally, this approach is organized into four major kernels functions, such as polynominal, linear, radial, and sigmode. In this study, the radial basic function (RBF) kernel [85,86] was applied with a 100 penalty parameter to be addressed in the ENVI software, and each kernel equation is listed in the following equations: where, x i , y i are training vectors, and g, d, and r are the user-controlled parameters of kernel function.

Simulation of LC Change
The CA-Markov hybrid approach is a suitable method for the simulation of LC change in places where cognition and characterizing landscape relationships are difficult [87]. This hybrid model has been widely used to effectively recognize and estimate landscape changes [88]. The CA-Markov model uses Markov chain matrices to identify the quantity of changes and Cellular Automata (CA) to spatially allocate these changes. The CA model addresses spatial allocation and the location of change via five parameters: (a) cell (b) neighborhood, (c) rules, (d) time, and (e) state [65,89]. In most cases, the steps of the CA-Markov model consist of [65]: (1) the classification of satellite images to generate LC maps; (2) the computing of transition area matrices; (3) the creation of transition potential images through driving parameters; (4) the estimation the model's capability to predict future changes based on evaluation indices; and (5) simulating the LC maps for future years.
The transition area matrix (TAM) was calculated using the Markov model. The TAM shows the number of pixels that are expected to shift from one LC type to another in the coming years. TAMs were calculated based on changes made in consecutive time periods (i.e., 1996-2006, 2006-2016, and 1996-2016) to show how each LC type was projected to shift. The TAM of 2006-2016 was used to simulate the LC projection to 2026, and the period of 1996-2016 was used to project the 2036 LC map.
The preparation of suitability maps is a difficult preliminary step in modeling LC change and depends on access to information and data [65]. In our study, the transition potential maps (TPMs) of LC types were prepared using multi-criteria evaluation (MCE), analytic hierarchy process (AHP) models, and fuzzy membership functions. The TPMs demonstrate the capacity of a cell to shift to a new category or to remain unchanged in each transition based on driving parameters [90]. In our study, the effective layers in the LC of the region, including slope, distance from roads, distance from water resources, distance from built-up areas, and distance from forest were chosen as the driving parameters based on the specialized knowledge and research history (e.g., [65,87,91]. The distance layers were standardized through the fuzzy membership method ( Table 3). The factors were rescaled to special ranges (0-255) based on particular functions. Table 3 shows the standardization properties of the criteria for each factor. After standardization, all of the factors were weighed using the AHP method [88]. In this method, the factors are considered and compared in a pairwise form based on their relative importance for use. After all of the possible combinations are compared between the two factors, the module, weights, and consistency are calculated, and if the value is less than 0.1, it means that the pairwise comparisons have an acceptable level of compatibility.

Land-Cover Modeling and Validation
In order to validate the CA-Markov model, the simulated LC map of the year 2016 was compared with the real map obtained from the classification of the satellite images from that year [92]. An accuracy of more than 80% indicates the model's simulation capability. The model was verified using several kappa variables: Klocation (the location accuracy of pixels in the simulation), Kstandard (a criterion of the ability of the model to achieve a complete classification), and Kno (the number of correctly classified pixels compared to the pixels that are expected to be classified correctly, without considering the quantity or the location error [65]. Therefore, researchers consider the Kno as a modified and more reliable version of Kstandard. The two indices-quantity disagreement and allocation disagreement, suggested by Pontius and Millones [93]-were also calculated and analyzed.

Accuracy Assessment
Accuracy assessment is one of the most fundamental tasks when LC data are prepared using remote sensing tools [94][95][96]. However, there is a challenge to find high-resolution data that can be used for the accuracy assessment. In this study, a total of 2844 sample points were designed for each of the LC-classified images for the years 1996, 2006, and 2016. A minimum of 300 sample points for each LC class and the user accuracy (UA), producer accuracy (PA), and overall accuracy (OA) were identified. The accuracy assessment of the classified images was prepared based on GPS points that had been collected from field observations using topographical maps developed by the Survey Department, 1998 (scale 1:25,000 and 1:50,000) [77], and high resolution Google Earth images (http://earth.google. com, accessed on 30 August 2021). The overall accuracy of the results obtained for the individual years were 85.61% (1996), 84.95% (2006), and 86.91% (2016).

LC Dynamics
Over the period of 1996-2016, remarkable LC changes were observed. The areas of other land, mainly settlement areas, increased, particularly in the second decade. Forest cover decreased slightly in the first period but increased in the second, whereas shrub land showed the opposite trends. Changes in the barren land, water body, sand, and grassland areas were minimal. Small parts of the northern area fluctuated among ice and snow

Spatial Transitions
The gains and losses of different LULC in different classes between 1996 and 2006 are presented in Table 5. Other areas(including settlements) increased by 18.09%, from 183.24 km 2 to 216.38 km 2 . This increase was mainly due to the conversion of 26.31 km 2 of cultivated land, 1.16 km 2 of forest, 3.2 km 2 shrub land, and 1.68 km 2 of sandy areas into

Spatial Transitions
The gains and losses of different LULC in different classes between 1996 and 2006 are presented in Table 5. Other areas(including settlements) increased by 18.09%, from 183.24 km 2 to 216.38 km 2 . This increase was mainly due to the conversion of 26.31 km 2 of

Spatial Transitions
The gains and losses of different LULC in different classes between 1996 and 2006 are presented in Table 5. Other areas(including settlements) increased by 18.09%, from 183.24 km 2 to 216.38 km 2 . This increase was mainly due to the conversion of 26.31 km 2 of cultivated land, 1.16 km 2 of forest, 3.2 km 2 shrub land, and 1.68 km 2 of sandy areas into other areas, most of which were settlement areas. However, cultivated land declined by 37 km 2 (from 6542 km 2 to 6504 km 2 ) during this period, which was mainly due to conversions to shrub land (32.98 km 2 ), sandy areas (17.96 km 2 ), forest (8.30 km 2 ), and water bodies (7.88 km 2 ). Similarly, despite the conversion of an 80 km 2 area from other classes (e.g., 74 km 2 from shrub land), forest cover witnessed an overall loss of 43 km 2 (9491 km 2 to 9447 km 2 ), with roughly 132 km 2 forest being converted to other classes. In contrast, shrub land increased by 90 km 2 , mainly because of the conversion of 36 km 2 of grassland areas into grassland areas. Sandy areas increased by 80 km 2 (from 476 km 2 to 556 km 2 ) mostly due to the conversion from forest and grassland areas. Barren land mainly increased due to reductions in snow and ice cover (Table 5 and Figure 4).    The major LC changes during this period include (a) overall increases in other areas in terms of settlement, barren land, forest, water body, and grass areas and (b) declines in cultivated land, shrub land, sandy areas, and ice/snow cover. Other areas with urban areas increased by 120.52 km 2 (from 216.38 km 2 to 336.9 km 2 ), mainly due to the conversion of cultivated land (97.44 km 2 ), forest (9.45 km 2 ), sand areas (5.86 km 2 ), and shrub land (4.78 km 2 ).
Large portions of cultivated land were converted into other areas, particularly urban settlement areas, forests (9.45 km 2 ), shrub land (4.78 km 2 ), barren land (1.67 km 2 ), and The major LC changes during this period include (a) overall increases in other areas in terms of settlement, barren land, forest, water body, and grass areas and (b) declines in cultivated land, shrub land, sandy areas, and ice/snow cover. Other areas with urban areas increased by 120.52 km 2 (from 216.38 km 2 to 336.9 km 2 ), mainly due to the conversion of cultivated land (97.44 km 2 ), forest (9.45 km 2 ), sand areas (5.86 km 2 ), and shrub land (4.78 km 2 ).
Similarly, 40 km 2 of water body areas were converted to sand areas, and 37 km 2 of sand areas were converted into water bodies. Based on this transition (Table 6), there were 9 km 2 increases in water body areas and an 80 km 2 sand area decrease. Around 31 km 2 of snow-and ice-covered areas in the Himalayan region was converted from barren land alone; however, barren land still increased from a 334 km 2 to a 350 km 2 area. The majority of the grassland areas that were acquired were gained from cultivated land, forest land, shrub area, barren land, and sand and snow and ice cover areas, and grassland areas increased by 181 km 2 (471 km 2 to 652 km 2 ) ( Table 6 and Figure 5). majority of the grassland areas that were acquired were gained from cultivated land, forest land, shrub area, barren land, and sand and snow and ice cover areas, and grassland areas increased by 181 km 2 (471 km 2 to 652 km 2 ) ( Table 6 and Figure 5).

Analysis of Transition Matrix
For the analysis and prediction of the LC data, we computed a transition potential matrix based on LC conditions for the period of 1996, 2006, and 2016 to identify how each LC class was projected to change in the years 2026 and 2036. During the simulation, several LCs were converted into each other while some remained almost constant over time. The transition probability matrix predicted that cultivated land will be converted into other (urban/built-up) areas and that shrub land will be converted into forest. This probability matrix showed the transition of each LC class (Appendix A Table A1).

Analysis of the Simulation Results
Major changes predicted for 2016-2036 include increases in other areas, including settlement, forest, and water bodies; declines in shrub land, barren land, and cultivated land; and fluctuations for sand and grassland. According to the actual LC statistics and simulations, urban areas are predicted to increase from 1.75% in 2016 to 2.58% by 2026 and to 3.08% by 2036 (Table 7, Figures 6 and 7). In contrast, cultivated land will decrease from 33.38% in 2016 to 32.01% by 2026 and to 31.63% by 2036. The simulations also suggest that forests will continue to increase and will increase from 50.33% in 2016 and to 50.74% and 51.76% by 2026 and 2036, respectively. However, shrub area will continue to decrease in 2026 and 2036, mainly because of its conversion to forests. Similarly, cultivated land will also continue to be converted into forests.

Discussion
Our results show that overall forest cover increased during the period of 1996-2016. Community-based forest management and other programs were responsible for the overall increase in forest cover in the study area. Similarly, the National Biodiversity Strategy and Action Plan provided a strategic roadmap for biodiversity conservation in Nepal [97]. As a result, by 2018, the country had a total of 12 National Parks, 1 wildlife reserve, 1 hunting reserve, 6 conservation areas, and 13 buffer zones, accounting for 23.39% (34,419 km 2 ) of the total land area (http://www.dnpwc.gov.np/en/, accessed on 30 August 2021). Additionally, different local and international level donor agencies such as IUCN, UNDP, WWF, and ICIMOD [98] supported environmental management and forest conservation. This increased forest cover has enhanced the conservation of endangered species and has enhanced animal biodiversity [99].
Forest area increases were mainly due to greater community level conservation practices and forest management strategies [22]. Further, urban area increases due to population growth resulted in migration from the highlands to the lowlands [55], causing abandoned cultivated land in upstream areas, which have turned into vegetation cover, leading to increases in forest cover [100,101]. Political conflict between 1996 and 2006 caused many people to migrate to more secure urban areas from agricultural areas in the mid-hills, resulting in the abandonment of agricultural lands and increased forest cover [102]. However, forest encroachment was higher in the Tarai region over the same period. Erosion, lowland flooding, urbanization, and deforestation are major causes of forest degradation in the lowlands (Tarai) of Nepal. The forest cover in Tarai was found to decrease between 2001 and 2010, and major losses occurred in the Kapilbastu district [103]. Similarly, our results also showed an overall decline in forest cover during the period from 1996 to 2006.
About 38% of people in Asia use wood fuel [2], and 70% of the total energy that is

Discussion
Our results show that overall forest cover increased during the period of 1996-2016. Community-based forest management and other programs were responsible for the overall increase in forest cover in the study area. Similarly, the National Biodiversity Strategy and Action Plan provided a strategic roadmap for biodiversity conservation in Nepal [97]. As a result, by 2018, the country had a total of 12 National Parks, 1 wildlife reserve, 1 hunting reserve, 6 conservation areas, and 13 buffer zones, accounting for 23.39% (34,419 km 2 ) of the total land area (http://www.dnpwc.gov.np/en/, accessed on 30 August 2021). Additionally, different local and international level donor agencies such as IUCN, UNDP, WWF, and ICIMOD [98] supported environmental management and forest conservation. This increased forest cover has enhanced the conservation of endangered species and has enhanced animal biodiversity [99].
Forest area increases were mainly due to greater community level conservation practices and forest management strategies [22]. Further, urban area increases due to population growth resulted in migration from the highlands to the lowlands [55], causing abandoned cultivated land in upstream areas, which have turned into vegetation cover, leading to increases in forest cover [100,101]. Political conflict between 1996 and 2006 caused many people to migrate to more secure urban areas from agricultural areas in the mid-hills, resulting in the abandonment of agricultural lands and increased forest cover [102]. However, forest encroachment was higher in the Tarai region over the same period. Erosion, lowland flooding, urbanization, and deforestation are major causes of forest degradation in the lowlands (Tarai) of Nepal. The forest cover in Tarai was found to decrease between 2001 and 2010, and major losses occurred in the Kapilbastu district [103]. Similarly, our results also showed an overall decline in forest cover during the period from 1996 to 2006.
About 38% of people in Asia use wood fuel [2], and 70% of the total energy that is produced by people in rural Nepal is from firewood [104]. However, the consumption of fuel wood for cooking has reduced by 3.3 times over the last decade and has been replaced by the use of liquefied petroleum gas [105]. Community forest programs provide strong support for energy management and carbon storage [20]. Nepal's decade of forestry (2014-2024) aims to conserve forest resources and to create urban greenery in cities and towns throughout the country. Operational plans for community-managed forest programs emphasize the control of forest fires [106], encroachment and illegal logging, forest awareness campaigns and monitoring, and silvicultural practices for the sustainable utilization of forest resources [107]. As a result, such programs should play an important role in mitigating greenhouse gases and reducing the impacts of climate change as well as preserving forest biodiversity.
Local-level plans and strategies also provide support for the future restoration of forests. The Tilottama municipality located in the Rupandehi district (Appendix B Figure A1: Map 1), for example, has introduced a plantation (about 0.3 million plants) program that will cover the period of 2017-2022. Butwal sub-metropolitan city has a similar program for the Tinau corridor and other areas, and the Sainamaina municipality has identified forest zones where they are promoting forest development on barren land. In total, 20 local governments have plantation programs in this region. These include roadside plantations, trees and green spaces in urban areas, the "one house two plantations program", riparian plantation programs, barren land plantations, and the promotion of agroforestry and private forests [108].
The Nepalese government is promoting agroforestry in the region and has put a forestry decade into place (2014-2024) with the motto of "one house one tree". Many of the suggested changes are in keeping with landscape and watershed level plans to adapt to climate change as set out in the 2010 National Adaptation Programme of Action (NAPA) and the 2011 Local Adaptation Plan of Action (LAPA). Plans at the national level emphasize forest preservation, the control of forest fires and invasive species, community-based integrated forest management, wetland and riverine forest conservation, and agroforestry. The action plan identifies local adaptations that are needed to mitigate climate change vulnerabilities while boosting resilience [109]. The trade agreement between the Nepalese government and the World Bank's Forest Carbon Partnership Facilities (FCPF) also encourage forest conservation in Nepal. Agroforestry practices (combined agriculture practices with forest and fruits) also support the maintenance of a green environment.
Community forestry has been successful in improving forest stocks, with increases in canopy cover [18], tree species diversity [110], growing stock, and biodiversity [111] being observed, all of which are due to reduced human pressure on forests [110]. Similar trends have been observed in provincial-level studies in Province One and in the Gandaki Province in Nepal [79,112]. Several studies have shown that forest cover of Nepal has increased in different regions. Tuladhar et al., 2019, explored forest cover changes in the entire catchment areas of the Bagmati river basin and found that forest cover increased by 4.1% between 1975 and 2005 in the middle section and overall for the basin. However, the forest change ratio differs in the upper and downstream sections of the catchment area [44]. Similarly, the overall forest cover area increased in the Phewa watershed during the periods of 1995-2017 [113] and 1975-2015 [114] and in the Tanahun district during the period of 1976-2015 [41]. At the national level, the annual deforestation rate decreased from 1.31% during the period of 1930-1975 to 0.01% in 2005-2014, and forest patch areas increased from 6925 km 2 in 1930 to 42,961 km 2 in 2014 [21], whereas national forest cover increased from 38% in 1978-1979 to 40.36% by 2015 [20]. A similar trend was observed in our study, where forested areas increased between 2006 and 2016.

Conclusions
Our analysis revealed that during 1996-2016, the study area witnessed declines in cultivated land, shrub land, sand areas, and ice/snow cover and increases in other LC, including urban areas, forest cover, barren land, water bodies, and grassland. Meanwhile, our predictions suggest that other areas, forested areas, and water bodies will continue to increase and that shrub land, barren land, and cultivated land will decline during the period of 2016-2036. We analyzed national forest management practices and policies and found multiple factors that were associated with historical and future trends of forest cover change. Community forest management practices, watershed management strategies, and forest carbon partnership facilities (FCPF) also encourage forest conservation at the local level. However, settlement expansion, forest encroachment, illegal logging, and forest fires are major challenges for the Tarai region (southern plains). As such, effective forest planning, particularly urban forest management, is essential for sustaining environmental equilibrium in the future. For this to occur, it is essential that forest protection and restoration plans and programs are formulated and are implemented at the national, provincial, and local levels in order to maintain the minimum amount of forested area.
We recognize a number of limitations of our study. We applied medium-resolution satellite imagery due challenges in capturing real-time high-resolution data [101]. Further, we only identified nine major LC classes. Accurate assessment of planted forests was difficult due to the low resolution of the available satellite images. Other uncertainties are related to natural factors, which include complex topography and seasonal snow-and ice-covered areas where ground truthing was not possible. We evaluated the LC changes at the provincial level, and we did not attempt to identify individual forest types or species, which we will consider in the future. Hence, we recommend the use of high-resolution satellite images for the finer evaluation of LC classes with forest species changes in the future. Existing policies should prioritize different forest species and tree-based ecosystem service-oriented management goals [115] in different geographical regions (mountain, hill, and Tarai). These will help maximize the benefits to local people and will help to maintain sustainable forest management in different ecological zones. Similarly, decision makers and scientists will need to apply suitable quantitative tools such as machine learning techniques and artificial intelligent (AI) methods for further data analysis and for the forecasting of environmental change.