Spatiotemporal change detection of land use land cover (LULC) in Fashiakhali wildlife sanctuary (FKWS) impact area, Bangladesh, employing multispectral images and GIS

Land cover change has posed significant concerns to biodiversity and climate change in Bangladesh and globally. Despite the country’s designation of forest regions as protected areas to conserve their valuable resources, deforestation and forest conversion remained unabated. Fashiakhali Wildlife Sanctuary (FKWS), a protected area in the Chittagong Hill Tracts, and its surrounding forested impact area have experienced considerable changes over the years, yet are deficient in extensive assessment. This study evaluated the land use land cover (LULC) changes in the FKWS impact area over almost 3 decades (1994–2021) using multispectral remotely sensed data. The Landsat images of 1994, 2001, 2010, and 2021 were classified using a maximum likelihood algorithm and analyzed for change detection. The comparative potential of vegetation indices, including Normalized Difference Vegetation Index (NDVI) and Soil Adjusted Vegetation Index (SAVI), in forest cover assessment, and the relationship between Land Surface Temperature (LST) and NDVI was also assessed. A significant forest cover loss of around 1117.17 ha (16%) was recorded in the FKWS impact area between 1994 and 2021, with the hugest proportion of 867.78 ha (12.24%) deforested in the first period (1994–2001). Agricultural land also declined by 593.73 ha (8.37%) within the entire period, despite its initial increase of 392.04 ha (5.53%) between 2001 and 2010, being the primary driver of earlier deforestation. However, in the recent decade (2010–2021), settlement expansion of 963.90 ha (13.59%) due to massive human migration in the area contributed to the most remarkable overall land cover change of 1731.51 ha (24.42%). Furthermore, NDVI provided a better and more accurate forest cover assessment than SAVI and was recommended to aid in the quick evaluation and monitoring of the future impacts of agriculture, settlement, and other sorts of land use on the forest cover. In tandem with the widely acknowledged issue of increased temperature due to climate change, an absolute negative correlation was found between the NDVI and LST, confirming the negative impact of climate change on forest loss in the FKWS impact area.


Introduction
The alarming rate of deforestation and its repercussions have raised concerns on a global scale (Caravaggio 2020;Hite and Seitz 2021). The illicit harvesting of natural forests, as well as the expansion of agricultural and development projects into previously forested areas, is the leading cause of forest loss (Chakravarty et al. 2012;Islam and Sato 2012;Hishe et al. 2021;Oluwajuwon et al. 2021). According to Lambin et al. (2001), the primary causes and reasons driving global and regional land cover changes are tropical deforestation, grassland alteration, agricultural intensification, and 1 3 urbanization. In particular, deforestation and forest degradation are more prevalent in developing countries due to less per capita land and poor adaptability (Iftekhar and Hoque 2005;Fagan et al. 2020).
Bangladesh is a small developing nation with a vast population size that was around 147 million as of mid-2007 and is predicted to be 222 million by 2050 (Streatfield et al. 2008;Farid et al. 2011), which has led to an increase in the threat to the country's forests (Mukul 2016;Hasnat et al. 2018;Singh et al. 2020). Accelerated population explosion has long drawn attention to the country's existing land use patterns (Biswas and Choudhury 2007;Reddy et al. 2016). Bangladesh's forest acreage has declined since the 1870s and now only covers approximately 2.33 million ha or less than 16% of the country's total land (Mukul et al. 2016;Nesha et al. 2021). Bangladesh's forests have been rapidly depleting for years, and the situation worsens daily. They might disappear shortly due to ongoing deforestation if sustainable forest management strategies are not implemented (Salam et al. 1999;Sunderland et al. 2011). The primary causes of forest degradation in Bangladesh are encroachment, illegal logging, increased fuelwood demand, conversion of forest land to settlements, agricultural expansion, and unregulated industrialization (Iftekhar and Hoque 2005;Rahman 2015; Reza and Hasan 2019;Mahmood et al. 2021). The government of Bangladesh has designated several forest regions as protected areas at various times to conserve valuable forest resources (Chowdhury and Koike 2010;Rahman and Islam 2021). Bangladesh has declared 49 protected areas, which span several broad categories, including national parks, wildlife sanctuaries, marine protected areas, vulture safe zones, safari parks, special biodiversity conservation areas, aviary parks, eco-parks, and botanical gardens (Rahman et al. 2016a(Rahman et al. , 2017. Protected areas are vital because they enhance carbon sequestration, maintain stable weather patterns, provide natural habitats for flora and wildlife, and prevent species extinction (Stolton and Dudley 2015;Xu et al. 2017).
The Fashiakhali Wildlife Sanctuary (FKWS), located within Fashiakhali Forest Reserve, is a protected area in the Chittagong Hill Tracts, Bangladesh and was declared a wildlife sanctuary in 2007 (Das et al. 2018). The Fashiakhali Forest Reserve has an undulating topography with numerous hills of varying heights that are covered with bushy vegetation, low valleys, waterbodies, and marshes in addition to natural and plantation forests (Uddin et al. 2011). Most of the forest covers in the FKWS are tropical wet, evergreen/ semi-evergreen, and deciduous forests, which provide habitats for high biodiversity (Billah et al. 2021). However, over time, people of the 30 villages that encircle the sanctuary have harmed this important ecosystem through the conversion of forested areas and rendered it more vulnerable to climate change (BFD 2015). In general, hill forests in Bangladesh, which comprise 43% of the nation's total forest acreage, have been heavily degraded and destroyed due to shifting cultivation, excessive illicit logging, encroachment, settlement, and rapid urbanization (Rahman et al. 2012;Ahammad and Stacey 2016;Hossain et al. 2020). The FKWS, like other forest ecosystems in Bangladesh, has lost tree populations, resulting in the demise of significant flora and fauna species (Das et al. 2018). Even if protected area principles, like those of wildlife sanctuaries, serve as the cornerstones of all regional biodiversity conservation initiatives, effective management and the proper execution of laws are mostly lacking in protected areas (Masum and Hasan 2020;Ullah et al. 2022a). As a result, the biodiversity in these crucial conservation areas continues to be impaired by the aspects of increasing land cover change.
The land cover describes the physical and environmental characteristics of the land surface area, including the presence of water, crops, forests, constructions, etc. (Turner 1994), while the human intention related to these characteristics is referred to as land use (Nendel et al. 2018). The land use data are required for evaluating land cover changes in a geographic site. The quantitative analysis of the land use patterns is required to obtain a greater understanding of land cover changes and aid decision-makers in setting program goals and adopting suitable strategies while remaining consistent with other disciplines of sustainability (Verburg and Chen 2000;Diouf and Lambin 2001;Lambin et al. 2001). Assessment of land cover change is a fundamental technique for evaluating changes in vegetation at different spatial and temporal dimensions. Empirical studies have shown that anthropogenic deforestation can impact the ecosystem in several ways, including destroying species' habitats, causing desertification and soil erosion, upsetting the water cycle, and increasing environmental risks due to forest fragmentation (Islam and Weil 2000;Romijn et al. 2015;Reza and Hasan 2019). The potential of forest biomass to store carbon is further impacted by changes in forest cover, which also affects the local climate by altering the diurnal temperature variation and raises the risks of global climate change (Newell and Stavins 2000;Sangermano et al. 2012). The primary objectives of research on global forest and environmental change are to catalog, monitor, and model the effects of changes on the forest environment, associated ecosystems, and forest properties at various scales. This procedure, known as change detection, relates to changes in land use land cover (LULC) (Reid et al. 2000;Abd El-Kawy et al. 2011). Changes in land cover, both naturally occurring and caused by human activities, have an impact on global and regional climate because of their interactions with terrestrial ecosystems, biodiversity, and landscape ecology (Houghton 1994;Reid et al. 2000). So, it follows that for environmental planning, conservation activities, comprehending the effects on the terrestrial ecosystem, and attaining sustainable development, monitoring and modeling of LULC is imperative (Redowan et al. 2014;Rawat and Kumar 2015). Furthermore, effective forest management requires evaluation and monitoring of forest land cover changes. Poorly vegetated or bare ground areas express concern about the forest (Park et al. 2017;Oluwajuwon et al. 2021). Meanwhile, in Bangladesh, the forestry sector places less emphasis on monitoring and determining the extent of land cover than on the ecological inventory of conventional forest plots. LULC assessment can be challenging for a vast forest ecosystem like FKWS, where satellite images allow for forecasting and more accurate results with less expenditure.
Remote sensing is a fundamental approach for studying spatial and temporal changes in LULC (Pastor-Guzman et al. 2018;Islam et al. 2018;Mamnun and Hossen 2020). Using medium to high-resolution satellite imageries in remote sensing and GIS applications, dynamic changes in the surface of the planet can periodically and easily be detected (Mallupattu and Reddy 2013;Rawat and Kumar 2015;Islam et al. 2021). Landsat satellite images, which have a spatial resolution of 30 ⋅ 30 m, can be a valuable economic data source to collect information from a particular area and identify changes in LULC (Gounaridis et al. 2018). Detecting and monitoring changes in forested land using widely available remote sensing data and methods also contribute to implementing climate change mitigation initiatives like Reducing Emissions from Deforestation and Forest Degradation (REDD+) and Clean Development Mechanism (CDM) (Sangermano et al. 2012;Potapov et al. 2014). The use of vegetation indices, such as the Normalized Difference Vegetation Index (NDVI) and Soil Adjusted Vegetation Index (SAVI), has frequently been reported for mapping forest cover (Nath and Acharjee 2013;Nath 2014;Bera and Prakash 2018;Islam et al. 2021;Oluwajuwon et al. 2021). Studies have also demonstrated how land surface temperature (LST) affects vegetation indices, such as NDVI (Alam et al. 2022;Hussain et al. 2022).
Moreover, until now, several studies have been conducted to detect changes in LULC of different Bangladesh's hill forests (Redowan et al. 2014;Islam et al. 2018;Mamnun and Hossen 2020;Masum and Hasan 2020;Hasnat 2021). Specifically, in the case of the wildlife sanctuary in Bangladesh, Chunati Wildlife Sanctuary was analyzed for land cover changes using satellite images, and it was discovered that the sanctuary had severely lost high-density forest cover (Islam et al. 2016). Another study found a strong correlation between forest fragmentation and the conversion of land to non-forest uses while using geo-informatics to track land use changes for Chunati Wildlife Sanctuary in Bangladesh (Rahman et al. 2016b). However, despite Billah et al. (2021)'s findings on LULC changes in Fashiakhali Forest Reserve as a whole and importantly its impact on human-elephant conflict, extensive scientific research on long-term LULC changes in the FKWS and its critically important surrounding impact area is lacking. In addition to the FKWS, the assessment of land use patterns in the impact region is crucial because it primarily describes how socio-economic and environmental structures function at the landscape level while making the necessary sustainability concessions to conserve the core ecosystem. Therefore, this study aims to narrow down this research gap and assess the LULC changes in the FKWS impact area. The objectives of this study are to assess the spatiotemporal changes in LULC in the FKWS and surrounding impact area over a 3-decade period, analyze the extent of the forest cover using NDVI and SAVI indices, and demonstrate the correlation between NDVI and LST indices within the FKWS.

Study area
Fashiakhali Wildlife Sanctuary (FKWS) is situated in Cox's Bazar District of Chittagong Division, located in the South-Eastern part of Bangladesh. The FKWS is situated within Fashiakhali Forest Reserve (3068.7 ha) and lies between 21°45' to 21°40' N and 92°4' to 92°8' E (Billah et al. 2021;BFD 2015). The FKWS was declared a wildlife sanctuary according to the Bangladesh Wildlife Amendment Act in 2007 and as a part of the Fashiakhali Forest Reserve, it covers a core protected area of 1302.52 ha, which is further divided into blocks: Dulahazra (287.50 ha), Ringbong (613.00 ha), and Fashiakhali (402.02 ha) (BFD 2015). However, the study area purposively includes the core FKWS and its surrounding impact area (Fig. 1). Hence, the considered FKWS impact area covers a total area of around 7093 ha with the inclusion of a buffer zone (1366 ha) and an impact zone (4384 ha) to the main core FKWS area. The FKWS impact area (core, buffer, and impact zones) was delineated by the management plan of the Bangladesh Forest Department in order to encourage co-management with the active involvement of locals to enhance their livelihoods, lessen forest degradation, and reduce the susceptibility to climate change. The impact area was identified to consider the 30 villages that encircle the designated FKWS territory and whose inhabitants have a history of harming ecosystems and making them more susceptible to climate change (BFD 2015). A moist tropical maritime climate prevails in the study area, where the average humidity is 79%, the mean annual rainfall is 741 mm, the mean annual temperature is 27 °C, and elevation ranges from 2 to 60 m with mostly flat terrain and scattered small hills of similar geology (BBS 2012). The FKWS region has a variety of different types of land cover. Forests, including natural stands and plantations, agricultural land, water bodies, aquaculture, and settlements are some examples of these land uses. Reportedly, the land cover has altered dramatically in each zone of the FKWS impact area as a result of anthropogenic activities over the past few decades, including migration, illegal felling and harvesting, land conversion to agricultural uses, and invasion (BFD 2015). The human settlement in the sanctuary area started with 112 people under only two villages, which now includes a total population of about 50,000 under 30 villages around the core FKWS area. Agriculture is the primary occupation in the area, with 60% of the population engaging in it, 15% in fishing, 20% as day laborers, and the remaining 30% in other jobs. The forest in the study area falls under the tropical rainforest category, and the dominant tree species are Garjan (Dipterocarpus turbinatus) and Dhakijam (Syzygium grande). Among the plantation species, Teak, Eucalyptus, and Acacia species are the most prevalent. This sanctuary is crucial for maintaining the habitat of a wide range of local endangered and threatened wildlife species, especially Asian Elephants (Elephas maximus). People frequently enter and extract forest resources without proper supervision, posing impacts on the natural regeneration of important trees along with impacts on the lives of other

Satellite images and field data
Landsat images of respective sensor and satellite types for 1994, 2001, 2010, and 2021 were obtained from USGS (United States Geological Survey) Earth Explorer (path/ row: 136/45). Since there was a paucity of cloud-free images accessible for the study area in USGS Earth Explorer's open sources, it was not feasible to obtain an equal year interval of Landsat images-a more detailed description of these satellite images is listed in Table 1. To avert the effect of cloud imaging in the land cover classification process, the satellite images used in this study were obtained during the winter season of the considered years: November-February. This choice aligns with previous remote sensing studies in Bangladesh, where most Landsat images from the winter period in the country exhibited negligible or no clouds (Islam et al. 2018Chowdhury et al. 2020;Billah et al. 2021). Since the study area is primarily an old-growth forest environment, seasonal change, such as sun sensor geometric variation, was considered low or insignificant in this study. A ground truthing was done-GPS (Global Positioning

Fig. 1 Study area-Fashiakhali Wildlife Sanctuary (FKWS) and its impact area
System) data for the predetermined specific land use/land cover classes (described in Table 2) were recorded and used as the basis for training samples for the land cover classification of the FKWS impact area. Additionally, Google Earth Pro's high-resolution real-time satellite images and a false-color composite image of the study area were used in addition to ground-based training samples from 2021 for a more precise land cover categorization.

Image pre-processing and classification
In this study, all the remote sensing analyses (i.e., image pre-processing, classification, change detection, vegetation indices, surface temperature mapping, and accuracy assessment) were performed in QGIS version 3.24.2 (Fig. 2)-a broadly used, freely available software (Congedo 2016). Before land cover classification, DN (digital numbers) values for all the considered raster bands of Landsat images were converted to TOA (top of atmosphere reflectance), which is the ratio of reflected and incident energy on a surface. When a satellite image is converted to reflectance, it can be mosaicked and compared to other satellite images from different sensors (e.g., Landsat-8 and Landsat-5), subsequently enhancing classification outcomes. Moreover, the atmospheric correction of reflectance values is critical during the pre-processing step because atmospheric factors, such as absorbance or dispersion, impact the electromagnetic energy captured by a satellite. The Semi-automatic Classification Plugin (SCP) in QGIS was used to perform a simple atmospheric correction following the DOS1 (Dark Object Subtraction 1) method. Given that the forest cover was the primary goal, it was essential to create a false-color composite of the study area to render the forest vegetation in red, which allowed to cross-check training samples (Fig. 3) and a prior visual understanding of the current conditions of forest cover in the study area. Finally, upon clipping the FKWS impact area from the satellite image grid, supervised classification following the maximum likelihood algorithm was carried out using the SCP plugin in QGIS (Congedo 2016;Karlsson et al. 2016). This involved creating ROIs (regions of interests) from training samples while considering the spectral variability or signature of specific land cover classes ( Table 2). The maximum likelihood algorithm is used to determine the likelihood of each pixel belonging to each of the predetermined LULC classes, and the pixel is then assigned to a class based on the highest probability. The area covered for each class per study year was determined by the pixel number.

Change detection and accuracy assessment
The study employed pixel-by-pixel cross-tabulation analysis (Jensen et al. 1987) to map and quantify pixel changes from one land use/land cover category to another within three different time frames : 1994-2001, 2001-2010, and 2010-2021. The post-processing tools of the SCP plugin in QGIS were used to perform the change detection analysis. Furthermore, the classification accuracy of the images was assessed to evaluate the validity of the information obtained from the data using stratified random sampling (Chowdhury et al. 2020;Islam et al. 2021). To reflect a sizable amount of validation data for each year, 149 randomly selected validation points from across all categories were purposefully sampled. The accuracy was tested using an error matrix for each year of investigation. Also, in a bid to measure the level of accuracy, two class-based accuracy assessment approaches (i.e., producer accuracy: PA and user accuracy: UA) and the total accuracy of the classified images were determined. While UA measures the ratio of the number of correctly classified pixels in each class by its total number of classified pixels, PA considers the total number of reference pixels in the said class as the denominator. Another common image classification accuracy measure is Kappa coefficient (K). Typically, its values range from 0 to 1, where the higher the value, the higher the agreement and accuracy (Billah et al. 2021;Islam et al. 2021;Hasnat 2021). This K statistic was also computed in this study using the following formula: where TP is total pixels, and TCP is total corrected pixels.

Vegetation indices
Several vegetation indices can be employed to detect and analyze the existence and extent of vegetation and forest cover. Among the commonest ones are NDVI and SAVI (Huete 2012;Vani and Mandla 2017;Huang et al. 2021;Islam et al. 2021). These indices are measures of vegetation and soil surface reflectance. In this research, the viability of these two widely used vegetation indices was therefore investigated to further determine the extent of forest cover and greenness in the ecosystem and its adjourning impact area. To do so, the NDVI and SAVI values of all pixels were extracted, respective to each year of assessment, using the raster calculator tool in QGIS based on the formulas below (Eqs. 2 and 3). NDVI formula was sourced from Huang et al. (2021), while the SAVI was from Huete (1988). In these equations, NIR is the DN value from the near infrared band; R is the DN value from the red band; L is the soil brightness correction factor. L values range from 0 for very high vegetation cover to 1 for very low vegetation cover, and the most typically used value is 0.5, which is for intermediate vegetation cover. Both NDVI and SAVI values range from + 1 to -1, with areas of no or sparse vegetation characteristically recording low values while highly vegetated or forest areas are attributed to higher values.

Land surface temperature (LST)
LST is an essential parameter that measures the emission of thermal radiance (temperature) from the land surface or, in vegetated areas, the canopy surface (Alam et al. 2022). The higher the surface temperature value for an area, the less probability of it having sufficient vegetation cover. To calculate the LST index for the FKWS impact area over the course of the study period, the TM thermal band 6 (10.4-12.5 μm), ETM + thermal band 6 (10.4-12.5 μm), and OLI thermal band 10 (10.6-11.19 μm) were used. The first stage involved converting the thermal band digital number (DN) values for each Landsat image, according to the respective sensor types, to sensor spectral radiance using the following formula (Eq. 4) (Amiri et al. 2009).
The following equation (Eq. 5) is then used to convert spectral radiance to temperature in Kelvin (Alam et al. 2022), and then to Celsius by subtracting 273.15. The calibration constants K 1 and K 2 were acquired from the Landsat data user's manual (Table 3).

Relationship between NDVI and LST
The linear relationship between the NDVI and LST was computed using R (R Core Team, 2013). The relationship between the two indices was investigated by randomly choosing 500 NDVI and LST point data within the FKWS boundary, where NDVI was regarded as the dependent variable and LST as the independent variable.

LULC patterns of FKWS impact area (1994-2021)
The LULC classification of the study area for the years 1994, 2001, 2010, and 2021 are shown in Fig. 4, while the area statistics of the LULC patterns are summarized for the respective years in Table 4. In 1994, forest accounted for Therefore, according to the LULC analysis of the 27-year Landsat data of the FKWS impact area, the forest remained the dominant land use category for all the assessment periods, but with a substantial decreasing trend. The record of a considerable tract of remnant forest cover in the protected area could be attributed to the co-management activities and Bangladesh Forest Department's initiatives to avert deforestation in the Cox's Bazar forest range, wherein our study area falls (Ullah et al. 2022a). Nonetheless, the recent influx of Rohingyas in Cox's Bazar is drastically changing the status of land usage and has resulted in decreased forest cover (Hassan et al. 2018;Hasan et al. 2021). Although agriculture took the next largest share of the reserve and its impact area up till 2010, being a principal livelihood source in Bangladesh (Hossain 2013;Misbahuzzaman and Smith-Hall 2015), its cover increase and degrading impact on the forest is reducing, with settlements now taking over. According to our findings, 1117.17 ha of forest land was depleted, agricultural and grassland had declined by 593.73 ha, while settlement had increased by 1731.51 ha between 1994 and 2021. Interestingly, only the settlement was responsible for destroying 404.64 ha of forest land in the FKWS impact area (Table 5). These findings partly align with Billah et al. (2021), who by 2015 found the primary cause of deforestation and land conversion in the study area to be the increasing demand for settlement as well as rural livelihoods. Nevertheless, these reports are not so different since we identified agricultural expansion for livelihoods to have predominantly driven forest loss by 2010 before being displaced by settlement expansion.

Land cover change detection of FKWS impact area
The change detection analysis for each land cover category is mapped in Fig. 5 and summarized in Table 5. The relative 1 3 changes in the land cover classes are graphically presented in Fig. 6. From the findings, it was observed that the land cover of the reserve and its impact area had undergone a considerable change within the assessed overall epoch 1994-2021. The forest, water, and agriculture/grassland classes predominantly experienced negative changes, although agriculture/ grassland had recorded a considerable positive change within 2001-2010. On the other hand, the land areas considered bare land and settlement recorded a positive change over the years. The most significant transitions between 1994 and 2021 were found for forest and settlement, with a total of about 16% loss and 24% gain, respectively. Of these changes, the hugest forest change was recorded within the earliest seven (7) years evaluated (i.e., 1994-2001). This period witnessed the removal of 867.78 ha of forest cover, accounting for over 12% relative to the initial cover. The substantial loss of forest vegetation commensurated the observed extensive conversion to bare land, which recorded almost a 12.5% (883.71 ha) increase. By that time, the rates of change in agricultural lands and settlement areas had not become prominent. Although the forest conversion had significantly decreased by the following epochs (i.e., 2001-2010 and 2010-2021), the already depleted land in the reserve was substantially transformed into agricultural and structural development uses. There was a conversion of approximately 392.0 ha and 521.3 ha to agriculture and settlement between 2001 and 2010. Both land uses were the leading causes of deforestation and degradation in the protected area. However, the trend became reversed for agriculture but doubled for settlement in the last evaluation period-between 2010 and 2021. That is, while almost 808.65 ha of land managed under an agricultural system was lost, settlement expansion had recorded a significant increase accounting for 14% relative to 2010, indicating the recent pressure of the rural communities and migrants to satisfy their needs for habitation. Furthermore, a cross-tabulation of pixel-by-pixel examination of Landsat images provides more discernible findings to corroborate the LULC change analysis of an area (Chowdhury et al. 2020). Therefore, the cross-tabulation of the land use changes over three distinct periods: 1994-2001, 2001-2010, and 2010-2021 are presented in Tables 6, 7 and 8. This analysis confirms that the most extensive forest cover depletion was within the first two periods. During these periods, the majority of the lost vegetated area had been immediately converted to cropland/grassland ( (Tables 6 and 7). Meanwhile, in these periods, the conversion of forest cover for settlement was relatively minimal, cumulating to 104.49%. These findings align with Billah et al. (2021), who assessed the changes in the study area without its impact area until 2015 and attributed forest loss mainly to agricultural land. The primary source of livelihood for the surrounding rural communities of FKWS is agriculture, principally, Jhum cultivation-a traditional method of shifting cultivation. People usually clear a forest patch into bare land and gradually convert it to agricultural land. This kind of widely adopted indigenous cropping practice has been likewise reported as the main driver of deforestation and land use change in other forest areas in Bangladesh, such as in Chittagong Hill Tracts (Rahman et al. 2012). Likewise, several studies across the country have ascribed the depletion of forest resources to the people's dependency on forest lands and products, encroachment for farming, illegal occupation and grabbing of forest land, and inadequate supervision by the forest department (Salam et al. 1999;Iftekhar and Hoque 2005;Biswas and Choudhury 2007;Xu et al. 2020;Billah et al. 2021).
In the recent decade, although the direct conversion of forest for agriculture in the protected area was no longer prominent, there has been continual cultivation of a large expanse of bare land. However, the rate decreased, with settlement expansion becoming the most prominent driver of forest cover loss (300.15 ha) (Table 8). During this period, there was critical pressure for human settlement due to a significantly huge migration and intrusion of about a million Rohingya refugees (Hassan et al. 2018;Hasan et al. 2021;Dampha et al. 2022;Ullah et al. 2022b). This might stimulate even more future forest degradation as a substantial proportion (798.93 ha) of the agricultural lands available for household farming was already   (Table 8), despite the dependence of the rural people on agriculture for survival (Hossain 2011). Already, almost 500 ha and 195 ha of agricultural and forest lands, respectively, had been converted to bare land with a high potential of being transformed into buildings in future years if the settlement rate remains unchecked. Nonetheless, there has been a considerable improvement in forest management, as demonstrated by the relatively reduced forest conversion rate in recent times. Also, in the last 2 decades, some agricultural and bare lands (e.g., 64.08 ha and 135.72 ha, respectively, between 2010 and 2021) had been re-established with forest stands. This development is similar to the report of Islam et al. (2021) in Nijhum dwip national park, which might be because of the involvement of community participation through the community forest management practice as an initiative by the Bangladesh Forest Department to increase the forest cover in the eastern upland of Bangladesh (Ahammad et al. 2019). However, such achievement is still minimal, considering the overall forest cover loss of the FKWS impact area and the Chittagong Hill Tracts. Therefore, while improving forest monitoring, dealing with the critical deforestation drivers is pivotal, including adequate

Accuracy of classified images
An important step after land cover classification is accuracy evaluation, which defines the validity of the resultant   map by assessing the errors for each class and generally for the whole classified image (Anderson et al. 1976;Congedo 2016;Rwanga and Ndambuki 2017). Therefore, the error matrices for the classified images are shown in Tables 9,  10, 11 and 12, while Table 13 summarizes their overall accuracy and Kappa measures. The almost perfect number of corrected pixels and the high user and producer accuracies (mostly 80-100%) for the land cover maps across all the years demonstrated the performance and validity of the land cover classification and products. This is even further established by the high overall accuracy of 87.64%, 90.46%, 93.34%, and 92.87%, with satisfactory kappa coefficients of 0. 82, 0.87, 0.91, and 0.91 for 1994, 2001, 2010, and 2020, respectively. Such a great accuracy could be largely attributed to the classification technique adopted for this study, as Mahmon et al. (2015) had earlier reported the Maximum Likelihood classifier to have produced the greatest accuracy when compared with some other classifiers like Mahalanobis Distance. Moreover, the accuracy assessment results correspond with the standard land cover classification accuracy of 85-90% proposed by Anderson et al. (1976) and are similar to some other land cover mapping studies in adjacent forests in Bangladesh (Hassan et al. 2018;Billah et al. 2021).

Assessment of forest cover based on vegetation indices (NDVI and SAVI)
Remotely sensed vegetation indices are straightforward and efficient methods for quantitative and qualitative assessments of vegetation cover, vitality, and growth dynamics (Xue and Su 2017;Pesaresi et al. 2020;Pasternak and Pawluszek-Filipiak 2022).  (Tables 5 and 7). Nevertheless, a further comparative validation of these indices' reports was necessary to establish their forest cover assessment deviations, as subsequently provided in Fig. 9. Table 14 further summarizes the estimation of forest land cover of the study area based on these indices. In this process, all pixel values of NDVI and SAVI were extracted from the classified image of 1994, and the lowest values representing forest cover predominantly within the FKWS boundary were carefully identified (i.e., 0.31 and 0.50, respectively). Therefore, for each image, pixels with NDVI ≥ 0.31 and SAVI ≥ 0.5 were assigned to the forest class and those with lower values to the non-forest category. A comparative evaluation of forest area estimates from both indices against the supervised classification output revealed reasonably similar trends but with a relative overestimation in 1994 and underestimation in 2001 by NDVI as well as overestimation in both 2001 and 2010 by SAVI (Fig. 9). Both indices largely undervalued the forest area percentage in 2001 and reported a relatively higher forest cover values for 2010. These estimation biases could be attributed to the sensitivity of the vegetation indices to effects of soil reflectance, soil appearance, atmosphere, and cloud shadows, necessitating remote sensing calibration (Huete 1988;Vani and Mandla  2017; Xue and Su 2017;Pasternak and Pawluszek-Filipiak 2022). Nonetheless, from the comparative deviations in Fig. 9, it could be inferred that NDVI identifies and estimates vegetation cover more accurately than SAVI, which is consistent with the findings of Islam et al. (2021). NDVI is therefore preferably recommended for rapid evaluation and monitoring of forest cover in Bangladesh (Redowan et al. 2014;Hassan et al. 2018;Islam et al. 2021).

Land surface temperature of FKWS impact area
Remotely sensed Land Surface Temperature (LST) has become one of the crucial indices to characterize and understand the thermal behavior of the earth cover in response to vegetation changes and associated climate processes (Trigo et al. 2008;Srivastava et al. 2018 Nevertheless, the general increase in land surface temperatures over the years in the protected area could not be considered to follow a linear trend, similar to the vegetation indices (Fig. 9). The mean LST value (24.02 °C) recorded in 2010 was lower than the average temperature (27.00 °C) found about 10 years earlier (in 2001) (Fig. 10). Such a trend is in tandem with the land use type and intensity experienced within a given period, besides other environmental (natural) factors like slope, soil type, altitude, soil moisture, and rainfall variation (Li et al. 2013;Jiang and Tian 2010;Tafesse and Suryabhagavan 2019). Unlike within 1994-2001, there was no considerable reduction in the forest cover, and even a large proportion of the bare land had become revegetated mostly by crops by 2010 (Figs. 4 and 5). This resulted in the relatively lower land surface thermal radiance emission reported in the year. Meanwhile, several studies have emphasized the comparative potential of agricultural lands (i.e. crops) to reduce land surface temperature over other (non-vegetated) LULC classes, mainly due to plants' considerable evapotranspiration capacity and their ability to scatter incident solar radiation and absorb heat, similar as the forests (Tafesse and Suryabhagavan 2019; Morsy and Aboelkhair 2021; Alam et al. 2022).

Relationship between LST and NDVI
NDVI has been widely demonstrated as a crucial, conventional indicator for assessing land surface temperature and dryness (Roberts et al. 2015;Morsy and Aboelkhair 2021;Alam et al. 2022). This tendency is contingent upon the observable relationship between the earth's near-infrared and thermal reflections. That is, there is commonly a correlation between NDVI and LST, although the direction of the relationship might be dependent on the LULC types (Jiang and Tian 2010;Zhang et al. 2016;Tafesse and Suryabhagavan 2019). Also, in FKWS, excluding the impact area, the relationship between LST and NDVI indices was investigated by randomly extracting and regressing 500 NDVI and LST data within the wildlife sanctuary's boundary. This selection was considered statistically representative of the pixels of the vegetation indices. The output of the regression analysis (i.e., scatterplot) for each of the assessed years is presented in Fig. 11. It is evident that the NDVI had an inverse or opposite relationship with LST, depicted by the negative regression coefficient between the indices, which aligned with several other studies (Jiang and Tian 2010;Morsy and Aboelkhair 2021;Alam et al. 2022). The coefficient of determination (R 2 ) between the indices varied across the years, ranging from 0.38 to 0.53; however, it could be generally asserted that the vegetation in the protected area assumes a moderate effect on the land surface temperature for all the years evaluated. Such finding therefore posits that the remaining LST variance would be explained by and is dependent on other site conditions in the area, like soil type, geological setting, elevation and climatic factors, including rainfall (Tafesse and Suryabhagavan 2019).

Conclusion
This study assessed the land use land cover changes in FKWS impact area in Bangladesh within nearly 3 decades (1994-2021) using remotely sensed multispectral imagery. The forest cover, predominantly accounting for almost 50% of the entire land area, drastically reduced over the years, especially within the first evaluation period (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001), when over 12% of the forest was lost. Up till 2010, forest conversion to bare lands primarily for agricultural activities was a prevalent driver of deforestation with almost 6% increase in agricultural land, due to the absolute dependence of the rural communities on household agriculture for livelihoods. However, in the last decade, settlement expansion increased sporadically by 14% as a result of the population expansion in the area from the massive migration of Rohingya refugees. This development not only resulted in forest depletion but also reduced the livelihood potential of the local people owing to the conversion of even agricultural lands to human settlements. Nonetheless, the trend of forest cover loss recently decreased as some proportions of agricultural lands and bare lands were re-established with forest trees. Based on our findings and observation, we recommend that the sustainability of such development and a transition to positive forest cover change in the area will largely depend on several forest management factors. Some of these include effective control of settlement expansion, continuous prevention of agricultural encroachment and land grabbing by settlers, and improved forest monitoring. Another important finding of this study was the applicability of the two vegetation indices (NDVI and SAVI) for evaluating or mapping forest cover changes. NDVI performed better and more accurately than SAVI in forest estimation and is most recommended for rapid monitoring of the protected area. Furthermore, the relationship (i.e., a negative correlation) existing between LST and NDVI was assessed in the forest area, which implies an indirect effect of the vegetation cover on the land surface temperature and vice versa in the area.
Funding Open Access funding enabled and organized by Projekt DEAL. No funds, grants, or other support was received for this study.
Data availability Dataset used during this study are available from the corresponding author on request.

Conflict of interest
The authors declare that no conflict of interests is directly or indirectly related to this work.
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 licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence 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 licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.