Co-management Effects on Forest Restoration in Protected Areas of Bangladesh: A Remote Sensing and GIS-based Analysis

1. Introduction

Forest ecosystems require innovative and sustainable management approaches to maintain their social, economic, and environmental functions [1]. In the face of growing dependency and excessive extraction of forest resources, sharing roles and responsibilities among stakeholders, from forest resource users to decision-makers, is imperative for sustainable forest management [2]. In recent years, influential actors such as governments and international donors have emphasized forest management by involving resource users and other stakeholders in forest conservation efforts [3].

Making decisions and managing forests through the equitable sharing of roles, responsibilities, authority, and entitlements among government and non-government actors is conceptualized as co-management [4]. Co-managed forests aim to achieve the dual objectives of conserving forest ecosystems while generating livelihood opportunities for forest-dependent communities. To ensure social and political equity in forest management, power, roles, and responsibilities are distributed fairly [4]. While several studies see [5,6,7] have been conducted to assess the social, economic, and environmental components of sustainability, few studies see [8,9] have considered the interactions between social and environmental components in the evaluation of the success of co-management, particularly in the context of developing countries.

The efficacy of co-management in protecting forest cover and biodiversity varies across the world. For example, co-management has successfully conserved forests in Ethiopia [10], Malawi [11], Honduras [12], and Nepal [12,13] compared to the forests without co-management. However, cases such as Miombo woodland, Malawi [14], and Haui Lu Luang, Thailand [15] have shown negligible success. To understand the underlying mechanisms of the success, many studies argue that co-management succeeds when local communities are engaged from the planning of a project to its implementation (e.g., co-management in Nepal) [16]. Since project planning and implementation contexts vary across countries, it is important to evaluate the success of co-management at a national scale for conserving and improving forest ecosystems. Also, such practices can help identify best forest management practices for replication across the world.

In 1979, Bangladesh launched its first community forestry project, Betagi-Pomora, amidst significant pressure on its forests [17]. Despite a nationwide logging ban from 1970 to 1980 failing to meet conservation goals [18], the Betagi-Pomora demonstrated community engagement’s potential in restoring forests and improving local livelihoods. However, legal constraints (i.e., Bangladesh Forest Act 1927 and the Wildlife Preservation Act 1974) and limited public participation hindered widespread co-management adoption [19]. To address this gap, the Bangladesh Forest Department (BFD) launched the Nishorgo Support Project (NSP) between 2004 and 2008, funded by the US Agency for International Development (USAID), covering five protected areas, including Lawachara National Park (LNP), Satchari National Park (SNP), Rema-Kalenga Wildlife Sanctuary (RKWS), Chunati Wildlife Sanctuary (CWS), and Teknaf Wildlife Sanctuary (TWS) [19]. The NSP aimed to involve local stakeholders in managing and conserving forests and implementing socio-economic programs to reduce forest dependency [20].

The NSP was scaled up as the Integrated Protected Area Co-management (IPAC) in 12 more protected areas in 2007, which ended in 2013 [21], focusing on stakeholder engagement and capacity building [22]. Subsequently, the Climate-Resilient Ecosystems and Livelihoods (CREL) project, running from 2013 to 2018, extended co-management to 21 protected areas [23]. Currently, the World Bank-funded Sustainable Forests and Livelihoods (SUFAL) project (2018–2023) seeks to expand co-management to 38 protected areas, enhancing community involvement and income options. Additionally, the USAID-funded Bangladesh Ecosystems Activity (July 2021–June 2026) follows similar objectives [24]. Thus, we are not considering the CREL and SUFAL timeframe (2016–2024) in the present study because our main goals are on the initial five protected areas.

Since all these projects have different management goals and the project implementation sites are distributed across different social-ecological zones, it is difficult to understand the success of the country’s co-management approach in terms of forest conservation at a national scale. Most of the studies evaluating the success of co-management have adopted a participatory approach to investigate the social and economic pillars of sustainability [6,25], reporting that co-management reduces the forest dependence of local communities for livelihoods. To bridge this knowledge gap, we have taken ‘forest cover change’ as a proxy of forest conservation to evaluate the forest conservation efficacy of the co-management [26] in five protected areas of Bangladesh (i.e., LNP, SNP, RKWS, CWS, and TWS) where co-management was officially piloted from the very beginning (i.e., 2004 during NSP) using remote sensing and GIS techniques.

Assessing the impact of co-management on vegetation cover change requires methodological precision in remote sensing and GIS techniques. Most of the studies primarily used time series of Landsat 5–8 see [10,12,27] and Sentinel [28] data to generate land-use maps or produce spectral vegetation indices. Because of the discontinuation of Landsat 5 in 2013 and Landsat 7 scan-line corrector (SLC) failure in 2003, it is challenging to create a continuous Landsat time series prior to Landsat 8 launch (2013). Many studies erroneously used Landsat 5, 7, and 8 data combinedly (e.g., [28,29]) to create a continuous time series without cross-calibrating locally, resulting in the bias in vegetation cover estimation [30]. The use of land-use maps developed from individual bands or a group of single bands to determine vegetation cover change can be misleading without proper ground-truthing in the heterogeneous vegetation canopy [31]. Spectral vegetation indices are generally better than land-use maps in determining vegetation cover gain/loss [32]. Normalized difference vegetation index (NDVI) is the most commonly used vegetation index in forest cover assessments [10,15]. However, NDVI is inherently sensitive to atmospheric conditions and background soil moisture conditions [31]. The use of multiple vegetation indices can effectively overcome this problem [31,32]. Moreover, some studies used MODIS (~1km resolution) images [15] to assess vegetation cover change over time. While MODIS images have a very high temporal resolution, they fail to detect small vegetation cover changes. It is thus paramount to use a robust method to assess the efficacy of co-management in enhancing vegetation cover.

Studies exploring forest cover change due to co-management in Bangladesh also showed some of the above-mentioned methodological limitations (e.g., [6,33]). The overall outcomes of these studies are ambiguous, and the effects appear to be specific to each protected area. For example, Islam et al., [6] observed an overall improvement in forest cover due to co-management in CWS, while Islam et al., [34] reported a rapid trend of forest cover loss in the same forest and under the same management scheme. Islam et al., [6] suggested no significant forest cover change under co-management in the LNP and the TWS. Using Landsat 5 and 8 data combinedly, Ahmed et al., [35] observed a negative relationship between forest cover and co-management in LNP, SNP, and RKWS. However, no studies have explored the spatial and temporal forest cover change in all five protected areas (LNP, SNP, RKWS, CWS, and TWS) during the successive project period (2004–2013) where the pilot co-management project began. Moreover, the efficacy of the co-management has not been compared with other forest areas without a protected area co-management policy.

A comprehensive study of forest cover changes due to co-management in the five protected areas compared to non-comanaged protected areas is critical to ascertain the efficacy of co-management in the forests of Bangladesh. To address this knowledge gap, this study assessed the effectiveness of co-management in improving forest cover in protected areas of Bangladesh. Using multiple spectral vegetation indices from Landsat images between 2003 and 2015, it compares co-managed areas (SNP, LNP, RKWS, CWS, and TWS) with a control site, Rajkandi Reserve Forest (RRF).

2. Materials and Methods

2.1. Description of the Study Protected Areas

Teknaf Wildlife Sanctuary (TWS), with an area of 11,615 ha (Figure 1a), includes 115 small villages with varying degrees of forest dependency [36]. TWS is dominated by evergreen and deciduous trees with a total of 535 species, of which 178 are herbs, 110 are shrubs, 150 are trees, 87 are climbers, and 10 are epiphytes (Uddin et al., 2013). The Asian elephant (Elephas maximus) is an important species of conservation at TWS. The sanctuary area has a moist subtropical climate (temperature: 15–32°C) with frequent and heavy rainfall (130–940 mm) during the monsoon season (May to October) [37]. In 2006, a total of eight co-management committees covering 48 villages around the TWS were formed.

Rema Kalenga Wildlife Sanctuary (RKWS) has an area of 1095 ha (Figure 1b). The sanctuary accompanies 45 villages with 24,000 inhabitants with a diverse range of dependency on forest resources [38]. A total of 634 plants, 167 birds, 7 amphibians, 18 reptiles, and 37 mammal species were recorded in this forest. The sanctuary is home to several primates, such as the globally endangered Hoolock hoolock and Trachypithecus phayrei, and globally vulnerable Ursus tibetanus, Nycticebus bengalensis, and Trachypithecus pileatus [39]. The area has a moist tropical climate characterized by high rainfall (4,162 mm), and the mean annual temperature ranges between 9.6 and 34.8°C [39]. The co-management council and co-management committee of the RKWS are primarily responsible for managing the sanctuary, which includes conserving the landscape and forest composition and increasing alternative income generation opportunities to lessen the forest dependency of local people for livelihoods and other resources [20].

Satchari National Park (SNP) has an area of 243 ha (Figure 1c). There are 73 villages in and around SNP with various degrees of involvement of local communities with the park (Mukul et al., 2017). Vegetation of the SNP is evergreen, with 200 ha of “natural” forest and the rest with secondary vegetation and planted forest. The national park is critical for protecting globally threatened mammals like the Western hoolock gibbon (Hylobates hoolock), capped langur (Trachypithecus pileatus), and Phayre‘s leaf Monkey (Trachypithecus phayrei) [25]. The average annual rainfall in this area is 4,162 mm and temperature usually range from 12°C to 32°C [40]. Co-management at SNP comprises a co-management council, co-management committee, people’s forum, village conservation forums, community patrol groups, and forest conservation clubs/youth clubs [25].

Lawachara National Park (LNP) was declared a protected forest in 1996 with an area of 1250 ha (Figure 1d). There are 26 villages, mostly with Indigenous settlements located within the LNP and the villagers have varying degrees of forest dependence [41]. LNP is a faunal biodiversity hotspot with 26 mammals (including five primates), 246 birds, 4 amphibians, and 6 reptile species [42]. The average air temperature varies between 26.8°C (February) and 36.1°C (June). The annual average rainfall of this area is ~4000 mm, with maximum rainfall from June to September [43]. In 2005, the LNP formed a co-management committee involving members from the local grassroots to administrative (BFD) level to create alternative income opportunities to reduce the forest dependency of local people.

Chunati Wildlife Sanctuary (CWS) has a hilly to mountainous physiography and an area of 7764 ha [44]. The climate of CWS is typically subtropical, with an average annual rainfall of 2,493mm, and the average temperature varies from 25.7°C in January to 32.5°C in May [21]. There are 60 villages in and around the CWS with various degrees of dependency on forest resource collection and betel leaf (Piper betel) cultivation [6]. The forest comprises 890 ha bush, 84 ha garjan (Dipterocarpus spp.), 13 ha small crown high forest, 11 ha open, 1458 ha plantation, 2761 scattered area, and 9 ha water bodies. CWS is one of the most important protected areas of Bangladesh as it is an important habitat of globally threatened Asian elephants (Elephas maximus). Co-management activities were implemented under two co-management committees to reduce over-exploitation and increase forest cover [45].

Rajkandi Reserve Forest (RRF) in north-eastern Bangladesh (24°15′0″ N, 91°55′0″ E) spans 5295.55 ha with medium to steep slopes and water streams. The brown, sandy clay loam to clay loam soils supports a diverse forest ecosystem. RRF experiences warm, humid weather, with temperatures averaging 27 °C in February to 36 °C in June. High humidity prevails, ranging from 74% in March to 89% in July. Annual rainfall averages 4000 mm, peaking during the Southwest monsoon from June to September [46].

2.2. Satellite Image Collection and Processing

2.2.1. Base Map Collection

The base maps for the six protected areas (LNP, SNP, RKWS, TWS, CWS, and RRF) were collected from the Nishorgo Support Project reports. These maps include boundary, mouza, and Upazila maps [47]. Shapefiles of these protected areas were collected from the World Database of Protected Areas [48].

2.2.2. Acquisition of Long-Term Satellite Image

To analyze the effect of co-management on vegetation cover change, we used Landsat 7 ETM+ satellite data from 2003 to 2015 for the six protected areas downloaded from the United States Geological Survey (USGS) official website (https://earthexplorer.usgs.gov). All the downloaded images had a spatial resolution of 30mx30m and a cloud cover of less than 20%. A total of 1105 images were downloaded for the five study areas (406 for TWS, 207 for SNP, 164 for RKWS, 99 for LNP, 112 for CWS, and 117 for RRF) from 2003–2015. For each study site in each year, we considered images from winter, summer, and rainy seasons to accommodate seasonal changes that can induce strong differences in vegetated scenes due to variations in phenological stages (e.g., plant senescence), physiological activity (e.g., stomatal conductance), and canopy morphology [49]. Landsat 7 data used for this study were standard level 2 terrain corrected (L2T) products that are already radiometrically adjusted and atmospherically corrected [50,51]. Through the LEDAPS project, the digital number values were finally converted to surface reflectance values [52].

2.2.3. Scan Line Error (SLC) Correction

Landsat 7 enhanced thematic mapper (ETM+) sensor had a failure of the scan line corrector (SLC) on 31 May 2003 resulting in images with wedge-shaped gaps on both sides of each scene. [53] developed a technique that can fill gaps in one scene with data from the nearest (same season) Landsat 7 scene. A linear transform is applied to the “filling” image to adjust it based on the standard deviation and mean values of each band, of each scene. We implemented this algorithm using the plugin “landsat_gapfill.sav” [54] in ENVI 5.3. Two images are required for performing the processes, one image termed as an SLC-off image which needs to be gap filled, and another by which we filled the gap are marked as SLC-on images in ENVI. We took the nearest (same season) minimum cloud cover contained scene as an SLC-off image in this plugin to fill the gap of SLC error images (Figure 2).

2.2.4. Filtering the Bad-Quality Pixels

Pixels of Landsat 7 images can be affected by land surface conditions, cloud contamination, or sensor conditions, which are termed as bad-quality pixels. Filtering bad pixels are important for time series analysis, particularly with vegetation indices. If, for example, bad pixels are included in a phenology study, the results might not show the true surface characteristics of seasonal vegetation growth [55]. Cloud-contaminated pixels will lower the indices values, and measurements such as the timing of ‘green up’ or peak maturity would appear later than they occurred [56]. To filter this affected pixel from a Landsat image a quality assessment (QA) band has been provided with every Landsat product. QA bits help improve the integrity of Landsat-derived science investigations by indicating which pixels might be affected by surface conditions, cloud contamination, or sensor conditions. The pixel values in the Level-2 QA band were converted into 16-bit binary form before use and a criteria table was developed (Table S1) to filter out the bad-quality pixels.

2.3. Calculation of Spectral Vegetation Indices

Several spectral vegetation indices have been proposed for the enhancement and extraction of vegetation information from satellite images. In this study, we used the most common indices such as vegetation index (NDVI), enhanced vegetation index (EVI), soil-adjusted vegetation index (SAVI), and modified soil-adjusted vegetation index (MSAVI) [57]. Using the surface reflectance values, NDVI was calculated as equation (1):

$$NDVI={\displaystyle \frac{NIR-R}{NIR+R}}$$

(1)

where NIR and R represent surface reflectance in near-infrared and red bands.

EVI is calculated as equation (2):

$$EVI=G{\displaystyle \frac{NIR-R}{NIR+C_{1}R-C_{2}B+L}}$$

(2)

where, NIR, R, and B represent surface reflectance in near-infrared, red, and blue bands, respectively. G is a gain factor, C1, C2 are coefficients to correct aerosol scattering in red band using blue band, and L is a soil adjustment factor [58]. Here, we used G = 2.5, C1 = 6.0, C2 = 7.5, and L =1 [57].

SAVI is used to correct NDVI for the influence of soil brightness in areas where vegetative cover is low. Landsat surface reflectance-derived SAVI is calculated as equation (3), a ratio between the R and NIR values with a soil brightness correction factor (L) defined as 0.5 to accommodate most land cover types [59]:

$$SAVI={\displaystyle \frac{NIR-R}{NIR+R+L}}\times (1+L)$$

(3)

The MSAVI minimizes the effect of bare soil on the SAVI. It is calculated as equation (4), a ratio between the R and NIR values with an inductive L function applied to maximize the reduction of soil effects on the vegetation signal [57]:

$$MSAVI=2\times \ast \rho NIR+1-{\displaystyle \frac{\sqrt{{\left(2\rho NIR+1\right)}^2-8(\rho NIR-\rho RED)}}{2}}$$

(4)

Here, $\rho$ is the atmospherically corrected or partially atmospherically corrected surface reflectance, and L is the canopy background adjustment that addresses nonlinear and different NIR and R radiant transfer through the canopy. At higher vegetation (NDVI > 0.6) this L factor value was considered 0 and in medium to lower vegetation cover (NDVI < 0.6) it was 0.5 [60].

Our study analyzed the net change (%) of vegetation indices, specifically NDVI, EVI, and MSAVI, relative to the base year of 2003. This calculation was carried out through a two-step process. Firstly, we subtracted the vegetation index raster of the years spanning from 2004 to 2015 from the base year of 2003. This process resulted in a newly produced raster containing values for every pixel index. These values can be either positive or negative, where positive values indicate an increase in indices compared to the base year, and negative values indicate a decrease. Secondly, the net change (%) of the vegetation indices was calculated annually by utilizing Equation 5 and considering the area.

$$Netchange\left(\%\right)={\displaystyle \frac{positivechangepixelnumber-Negativechangepixelnumber}{Totalpixelnumber}}\times 100$$

(5)

2.4. Land Cover Classification

Three commonly used vegetation indices, including NDVI, MSAVI, and EVI, were calculated as an indicator of LULCC over time in five protected areas that were considered in pilot projects for protected area co-management in Bangladesh. In our study, we also collected and measured the SAVI indices. However, because they were highly correlated to every MSAVI value, we excluded them from our analysis.

Based on spectral vegetation indices calculated from Landsat–7 ETM+ satellite images, each protected area used in this study was grouped into one of five different types of land use: a) barren, with waterbodies, and built-up; b) agriculture; c) a mix of agriculture and forest; d) dense forest; and e) highly dense forest. The spectral index threshold values (NDVI, MSAVI, and EVI) were used to calculate this land use class. We took the standard vegetation threshold values (Table 1) from the USGS website and adjusted the range for our study areas by testing several points implied on the same year‘s Google maps.

2.5. Statistical Analysis

Satellite image bulk downloading, pre-processing, spectral indices calculations, and statistical analyses were conducted in the R statistical software (version 4.0.3) [61]. To test the significance and the magnitude (increasing or decreasing) of trends in land use and land cover change (LULCC) in the protected areas of Bangladesh since the initiation of the co-management approach, the Mann-Kendall (M-K) trend test and the Theil-Sen slope estimator were utilized in this study respectively. Mann-Kendal method, of vegetation cover change, was conducted using the “Kendall” R-package [62]; slope analysis was done using “Robslopes” R-package [63].

3. Results

Land cover classification based on NDVI, MSAVI, and EVI values indicated that in all six protected areas, the largest proportion of land was covered by dense forests (~41–71%) followed by forest-agriculture mosaic (~19–47%), agriculture (~4–8%), highly dense forests (~0.5–2.2%), and barren-building-waterbody (~0.3–0.4%). The overall pattern in land cover was consistent for the three vegetation indices but the proportion of land cover classes varied among the protected areas (Table 2 and Figure 3).

3.1. Overall Trend of Land Cover Change in the Six Study Areas over the co-Management Period

Our overall results indicated that the proportion of dense forest decreased significantly (MSAVI and EVI: p < 0.01; NDVI: p = 0.07) over the co-management period in five protected areas, with higher slopes (–3.7 to –0.96) compared to other land cover classes. RRF showed similarly large decreases in dense forest areas in non-comanaged areas, with NDVI, EVI, and MSAVI ranging between –2.88 to –2.36 (p < 0.01) for slopes. Conversely, both agriculture and forest-agriculture mosaic land cover classes increased significantly, with agriculture having lower slope values (0.01–0.09, p < 0.05) than those of forest-agriculture mosaic (1–2.53; MSAVI: p <0.01) in co-management sites. On the other hand, highly dense forest and barren-building-waterbody land cover classes remained unchanged with fluctuating (plus/minus) but extremely low slope values. For the highly dense forest cover class, the slope values ranged from –6.9x10^-2 to 1.3e^-17, and for barren-building-waterbody classes, slope values ranged from 9.8e^-16 to –2.1e^-17 (see-Table 1 and Figure 3). In the control site RRF, the Barren-building-waterbody showed minimal change with non-significant slope values (NDVI: 0.024, MSAVI: 0.05, EVI: 0.02). Agriculture exhibited a positive trend with significant slope values (NDVI: 1.35, MSAVI: 1.02, EVI: 0.99). Forest-agriculture mosaic displayed significant growth with highly significant slope values (NDVI: 2.44, MSAVI: 2.22, EVI: 2.09). Dense forests experienced a significant decline (p < 0.001) with slope values (NDVI: –2.88, MSAVI: –2.63, EVI: –2.36). Highly dense forest showed non-significant trends (NDVI: –0.47, MSAVI: –0.11, EVI: –0.13) with p-values ranging from 0.55 to 0.20, indicating stability in this land cover class within the RRF.

3.2. Protected Area-Specific Trends in Land Cover Change during the Co-Management Period

Our analysis showed highly variable land cover change among the protected areas during the co-management period. The proportion of dense forests declined with time in the SNP, TWS, RKWS, and CWS while increasing in LNP. On the other hand, the proportion of agriculture and agriculture-forest mosaic increased in the SNP, TWS, RKWS, and CWS while declining in LNP. Patterns of highly dense and barren-building-waterbody were similar across the protected areas during the study period.

Although the dense forest cover was dominant, the proportion of dense forest in SNP reduced significantly [NDVI slope: –2.08 (p < 0.01); MSAVI: –3.64 (p < 0.01); EVI: –1.78 (p = 0.02)] from ~90% to 0.5%–1% (Figure 3d, 3i, and 3n). The largest decline slope in dense forest cover was observed in RKWS [NDVI slope: –4.58 (p < 0.01); MSAVI slope: –4.27 (p < 0.05)]. Similarly, in CWS [EVI slope: –1.11 (p < 0.05)] and TWS [EVI slope: –1.11 (p < 0.05)], the proportion of dense forest cover declined significantly. The proportion of dense forests in CWS, as indicated by NDVI, decreased from ~70% in 2003 to ~34.5% in 2015. The MSAVI and EVI also showed similar decreasing trends of dense forest cover with slightly different magnitudes in TWS throughout the period. Unlike other protected areas, in LNP, there was a significant increase [MSAVI slope: 6.55 (p < 0.01)] in the dense forest from ~1% (in 2003) to ~80% (in 2015); however, non-significant trends were found in NDVI and EVI (see-Table 1).

Similarly, forest-agriculture mosaic in SNP increased significantly (Table 1) from 10–17% in 2003 to 94.5–98% in 2015 (see-Figure 3c, 3h and 3m). A similar pattern was also observed for forest-agriculture mosaic in RKWS [NDVI slope: 5.42 (p < 0.01); MSAVI slope: 3.92 (p = 0.01); EVI slope: 1.14 (p = 0.04)] throughout the co-management period (Figure 3c, 3h, 3m), with a decreasing trend in the first few years and then increased to reach ~80–90% cover. Forest-agriculture mosaic in CWS also showed a significant increasing trend [NDVI slope: 2.43 (p = 0.20); MSAVI slope: 1.34 (p = 0.01); EVI slope: 2.94 (p < 0.01)]. In LNP, however, this pattern was the opposite and showed a significant decreasing trend [NDVI slope: –0.017 (p = 0.85); MSAVI slope: –4.62 (p ≤ 0.01); EVI slope: –0.22 (p = 0.85)] throughout the study period.

The agriculture land cover in SNP increased slightly from 0% (2003) to 4.5-7.5% (2015) with a slope ranging from –0.12 to 0.02 and was only significant (Table 1) for NDVI-based classification (slope = 0.021; p = 0.01) (Figure 3b, 3g, 3i). Although the overall trend was not significant in RKWS, NDVI-based classification showed agriculture land cover was below 1% in 2003, increased to ~3% during 2005—2006, decreased to ~1%, and remained the same afterward. MSAVI and EVI-based classification also showed a similar pattern for RKWS. In TWS, the agriculture land cover trend showed a significant positive slope with values ranging from 0.57–1.87. Agriculture land cover increased significantly [NDVI slope: 0.25 (p = 0.03); MSAVI slope: 0.91 (p < 0.01); EVI slope: 0.63 (p = 0.02)] in CWS during the co-management. In 2003, the proportion of agricultural land use accounted for ~4%, ~19%, and 11.5%, while in 2015, they reached ~24%, ~25%, and ~28% based on NDVI, MSAVI, and EVI classifications, respectively. Unlike other protected areas, in LNP, there was a significant decreasing trend in agricultural land cover throughout the study period [NDVI slope: –0.0006 (p = 0.31); MSAVI slope: –0.61 (p = 0.03); EVI slope: –0.023 (p = 0.50)].

There were no significant trends in highly dense forest land cover in SNP in 2003–2015 except a slight temporary gain in 2007—2009 (~2% based on MSAVI) (Figure 3e, 3j, 3o) with trend slope values ranging from 0.00—0.023 (Table 1). In CWS, highly dense forest cover decreased significantly with very low slope values [NDVI slope: —0.04 (p ≤ 0.01); MSAVI slope: —0.02 (p ≤ 0.01); EVI slope: —0.009 (p = 0.20)]. The proportion of highly dense forest cover in RKWS was below ~1% at the beginning of co-management and remained the same throughout the study period with no significant trend. Likewise, based on both NDVI and EVI, the proportion of highly dense forests in TWS declined from 6–13% (2003) to below ~1% (2015) although the trend was not significant. EVI-based classification showed a marginally significant (p = 0.09) decline with a low slope value (–0.016). Unlike other protected areas, a significant increasing trend was observed for highly dense forest cover in LNP from 2003 (0%) to 2015 (~5.5%) based on MSAVI (slope = 0.03, p = 0.03), while NDVI-based trend (slope = 0.158, p = 0.09) was marginally significant.

Moreover, no trend was observed for the proportion of barren-building-waterbody class over time in SNP, CWS, LNP, and RKWS, which was below ~1% at the inception of the co-management approach in 2004 and remained the same throughout the study period based on all spectral indices (Table 1; Figs 3a, 3f, 3k). In contrast, in TWS, this land cover class increased significantly [NDVI slope: 1.73 (p < 0.01); MSAVI slope: 0.45 (p < 0.01)] throughout the study period. The proportion of barren-building-water bodies increased from ~1% to ~17.5% (NDVI), ~1% to ~5.5% (MSAVI), and ~1% to ~4.5% (EVI).

In RRF, our analysis of land cover dynamics revealed distinctive patterns across various classes. Barren-building-waterbody exhibited minimal changes with non-significant trends (NDVI slope: 0.024, p = 0.92; MSAVI slope: 0.02, p = 0.51; EVI slope: 0.02, p = 0.87) (Table 1, Figure 3a, 3f, 3k). Agriculture showed a significant positive trend (NDVI slope: 1.35, p = 0.08; MSAVI slope: 1.02, p = 0.16; EVI slope: 0.99, p = 0.05), indicating an increase over the study period (Table 1, Figure 3b, 3g, 3l). Forest-agriculture mosaic experienced substantial growth with highly significant trends across all spectral indices (NDVI slope: 2.44, p < 0.01; MSAVI slope: 2.22, p = 0.07; EVI slope: 2.09, p < 0.001) (Table 1, Figure 3c, 3h, 3m). Dense forests exhibited a significant decline (NDVI slope: –2.88, p < 0.001; MSAVI: –2.63, p = 0.01; EVI: –2.36, p < 0.01), indicating a noteworthy reduction in coverage (Table 1, Figure 3d, 3i, 3n). Highly dense forest, however, demonstrated non-significant trends (NDVI slope: –0.47, p = 0.55; MSAVI slope: –0.46, p = 0.20; EVI slope: –0.52, p = 0.09) in this land cover class within the RRF (Table 1, Figure 3e, 3j, 3o).

3.3. Pixel-to-Pixel Changes of Vegetation Indices Over Time

The pixel-to-pixel change expressed as the percentage of area change of NDVI, MSAVI, and EVI values to the 2003-pixel values provided detailed dynamics of land cover change in the study areas over the co-management period. Starting from 2004, all study areas showed a similar pattern until 2006, and they diverged during 2011–2012, where CWS and TWS had net negative change and LNP, SNP, and RKWS had a net positive change in NDVI, MSAVI, and EVI values (Figure 4). After that, the pattern was reversed for all protected areas and turned into a negative change except for TWS. From 2013 onwards, the negative change continued for CWS and SNP, TWS oscillated between positive and negative changes, and LNP and RKWS started to show positive change for the rest of the study period. On the other hand, our non-co-managed protected area RRF exhibit a continuous negative trend from 2003 to 2015.

3.4. Forest Cover Change from Global Forest Watch Data

Historical forest cover loss (ha/year) for five protected areas of Bangladesh included in this study was retrieved from the Global Forest Watch database to assess the general deforestation trend (2003–2015) in the area (Figure 5). Over 10 years, the highest total forest cover loss occurred in CWS (158.77 ha) while the lowest was in SNP (0.49 ha). In LNP, total forest cover loss was estimated as 15.45 ha throughout 2005–2015 with a peak loss of 5.8±0.44 ha (mean±SD) in 2009. In CWS the peak forest cover loss occurred in 2006 (36.2±1.3 ha), 2009 (24.2±1.3 ha), and 2012 (36±1.7 ha). In RKWS, total forest cover loss over the period was 128.42 ha with a peak loss of 68.6±5.36 ha in 2007. In TWS total of 154.42 ha and the highest loss (16.4±1.3 ha) occurred in 2011 (Figure 5). In the control site RRF the total forest cover loss over the study period amounted to 77.8 ha. The most substantial annual loss occurred in 2006, reaching 48.2±2.17 ha. Notably, 2008 and 2009 witnessed significant losses of 34.6±0.8 ha and 37.2±4.38 ha, respectively (Figure 5).

4. Discussion

4.1. Co-Management Approach Reduces Forest Cover in Major Protected Areas

A decline in dense forests was found in SNP, TWS, RKWS, and CWS during the co-management period (2004–2015), with the CWS also showing a significant decrease in highly dense forests (Figure 3). This is consistent with findings from prior studies that reported deforestation during certain co-management periods in the TWS [64], RKWS [65,66], SNP [67], and CWS [33,44,68,69]. One of the major drivers for the reduced forest cover was the expansion of agricultural land, with the overall increase of agriculture and forest-agriculture mix land cover being observed for most protected areas (Figure 3; Table 2). High economic return was likely the major motivation for the conversion of forests to agricultural land. With agricultural land, local communities can generate more annual income by growing cash crops with short-term returns, such as betel leaf [33,64]. Due to the similar economic incentive, Islam et al., [70] found that the co-management approach reduced forest biodiversity because fast-growing exotic valuable timber species were preferred over slow-growing native tree species in tree plantation programs in protected areas.

Urban expansion, forest dependency, and ineffective policy enforcement also play important roles in the decreased forest cover in the protected areas. Land expansion due to increased population may result in land encroachment for human settlement and infrastructure development [64], which require tree removal and land clearance. In addition, local communities, in particular, indigenous people, have a high dependency on forests and forest resources [64] and directly contribute to forest cover loss through a wide variety of activities, such as illegal logging and grazing [33]. The high forest dependency of local households may be due to the lack of alternative non-forest income sources, suggesting the necessity of reducing the forest dependency of local people to enhance forest conservation. Studies on RKWS [38] and SNP [71] revealed that the co-management approach shifted strong forest dependency toward alternative income generation activities, which may potentially reduce the loss of forest cover. Ineffective co-management policy enforcement may also contribute to the failure of forest conservation in the protected areas. The different demands and mistrust among different stakeholders and limited engagement of local communities in decision-making and enforcement processes [44,72] likely resulted in the local communities’ fear of losing the forest land and the resources and unwillingness of protecting the forests [65]. In addition, since long-term land tenure is not available and the benefits of protecting forests are unclear [21], the incentive for forest protection was weak for local people. Other land cover types, such as barren land, buildings, and waterbodies, were found to gradually increase during the co-management period, being consistent with the increased vacant land from 2005–2015 in the study conducted by Islam et al., [44]. The reforestation practice is needed to convert the barren land into forest areas to increase forest cover and biodiversity.

4.2. Co-Management Approach Increases Forest Cover in LNP

In general, co-management approach has been effective in increasing highly dense and dense forests in LNP during 2004–2015 (see-Figure 3; Table 2). The results contradict the findings from Islam et al., [6] that indicated a decreased forest cover in LNP from 2011–2017 due to illegal logging but are consistent with recent studies [68] that revealed the fluctuating but increasing trend of dense forest cover during the co-management period. The increased dense and highly dense forests were mainly due to the conversion of agricultural and forest-agricultural mixed land into forests, indicated by the significant declines of agricultural and forest-agricultural mosaic land areas (Figure 3; Table 2). Such land cover changes had likely resulted from the creation of alternative income generation activities and reduced forest dependency and anthropogenic effects of local communities in LNP [73,74]. As the co-management approach provided alternative income-generating activities, such as training and educational programs, and diverse employment opportunities, such as jobs associated with ecotourism, that generated higher and stable income and/or returns but lower risks relative to intense forest wood extraction and illegal logging, deforestation activities in the protected areas decreased [73]. The increased forest cover was also likely attributed to the experimental plantation in the LNP by the Bangladesh Forest Research Institute (BFRI) in 2015 [45]. In addition, economic incentives (e.g., paid jobs associated with forest conservation) were provided for forest protection and close participation and engagement of local households in forest conservation, and decision-making processes were encouraged in the co-management approach in the LNP; thus, forest cover and biodiversity both improved under the co-management practice [68,75]. However, some studies have shown a decreased biodiversity in the LNP during the co-management period due to the plantation of fast-growing valuable timber species (e.g., Acacia auriculiformis) as a monocultural tree plantation [70]. Overall, the involvement of local communities in the co-management has been found to be the key factor in improving the livelihood of the communities from financial, social, physical, natural, and human aspects and enhancing sustainable forest restoration [74,76]. In addition, the moderate to high scores that the forest conservation outcomes received in the survey of local people indicate high satisfaction with the co-management practices [77]. This indicates the importance of inclusion and/or enforcement of indigenous and mixed-species plantations in reforestation practice in future co-management plans.

4.3. Co-management Approach Affects Total Vegetation Cover in Pilot Protected Areas for Co-management

In general, co-management generally increased vegetation cover in LNP, SNP, and RKWS relative to the pre-co-management period (2003), while a decline of vegetation cover was observed in CWS and TWS relative to vegetation cover in 2003 (Figure 4). The peak vegetation loss in CWS is consistent with the forest cover loss shown in Global Forest Watch data, with the peak vegetation loss being found around 2006, 2009, and 2012 (Figure 4 and Figure 5). The LNP had major vegetation cover declines in 2006 and 2012, approximately consistent with observations from Global Forest Watch. We suspect that the decreased vegetation cover was possibly due to the illegal tree removal in 2011 [45]. The net positive values of vegetation indices (NDVI, EVI, and MSAVI) in LNP, SNP, and RKWS were likely due to the increase of total vegetation cover of trees in the forests and crops in agricultural land. However, enhancing forest cover needs to be prioritized to preserve dense forests and the biodiversity that is associated with the forests [68,75].

4.4. Forest and Total Vegetation Cover in Forests without Co-management Practices

Recent studies have focused on the effects of co-management practices on the temporal changes in land use and vegetation cover in protected areas [65,67,68,72]; however, RRF, which is the forest reserve without the co-management policy and served as the control forest in this study, has received little attention. Our study revealed that the RRF showed a decline in dense forests but an increase in agriculture and forest-agriculture mosaic lands during 2003–2015 (Table 1; Figure 3). In addition, the total forest cover of the RRF decreased during this time (Figure 4 and Figure 5). The findings were generally consistent with the observed and projected decreasing forest cover in the RRF, with the mean forest cover declining from 13782 ha in 1988–2018 to 11986 ha in 2100 [35], Haque et al., [46] reported the rich angiosperm flora in the Rajkandi Reserve Forest (RRF) and showed that the number of angiosperm species in the RRF was higher than in CWS but lower than in RKWS and LNP. However, 25 species were identified as threatened, indicating the urgent need to protect the forest cover and the species to enhance biodiversity [46]. Since the effectiveness of the co-management practice strongly varied with protected areas, establishing policies and practices that are specific to the RRF is essential. The results of our study suggest the importance of including reforestation and local community engagement in forest management policies.

5. Conclusions

Our study examined the impact of co-management on forest cover in five pilot protected areas where co-management begins its journey in Bangladesh. We identified significant decreases in dense forest coverage alongside increases in agriculture and forest-agriculture mosaic using vegetation indices retrieved from remote sensing. These findings emphasized the urgent need for enhanced co-management strategies to address deforestation drivers effectively. Notably, the research filled a critical gap in understanding overall forest cover changes resulting from co-management, providing valuable insights for sustainable forest management in Bangladesh and globally. Targeted conservation efforts are imperative, particularly in areas experiencing substantial forest loss such as CWS and SNP. The study underscores the importance of adaptive forest governance frameworks that prioritize community engagement and ecosystem resilience for effective deforestation mitigation and biodiversity conservation.