Study area
Chure-Terai-Madhesh landscape (CTML) covers the entire elephant distribution range in Nepal. The CTML spreads across 25 districts and covers an area of 42,456 km2 (Fig. 1). CTML comprises five physiographic units i.e., Chure hills (34.4%); Chure narrow gorges (2.2%); Dun/Inner Tarai (8.4%); Bhavar region (14.9%); and Tarai Madhesh (40%). Forty-eight percent of the landscape comprises agriculture and settlement; 47.16% forest, shrub-land, and grassland; and the rest 4.65% river and riverbed 38. CMTL is a part of the global biodiversity hotspot 39 and provides essential environmental services such as groundwater recharge for more than half of Nepal's human population (~ 15 million) 40,41. The major habitat types are Himalayan subtropical broadleaved forests, b. Gangetic plains and moist deciduous forest, and c. Terai-Duar savannas and grassland. Apart from elephants, the study area is also a refuge for several endangered large mammals, including the tiger (Panthera tigris), greater one-horned rhinoceros (Rhinoceros unicornis), Gaur (Bos gaurus), and wild buffalo (Bubalus bubalis arnee). The annual rainfall ranges between 1,138 mm to 2,680 mm, with over 80% of the rain occurring during monsoon months 42. The altitudinal range lies between 60-1500 meters 43. CTML is densely populated with an average human density of 392 persons/km2 40. Sixty percent of the people depend on subsistence agriculture and are involved in farm and off-farm-based livelihood activities (Chaudhary and Subedi, 2019). Paddy (Oryza sativa), maize (Zea mays), wheat (Triticum aestivum), lentils (Lens culinaris) are some major food crops, where jackfruit (Artocarpus heterophyllus), mangoes (Mangifera indica), bananas (Musa acuminata) are some fruit crops farmed in the area 44. Large-scale linear infrastructure projects and mining activities are the major drivers of deforestation and habitat fragmentation in the landscape.
We divided CTML into four regions (Eastern, Central, Western and Far-western) of similar size to assess the extent of forest loss (Table 1). Thus, elephants are distributed in four population clusters with limited connectivity viz. eastern population (Mechi River to Kamala River), (b) a central population (Kamala River to Narayani River), (c) western population (Narayani River to Western boundary of Dang district), and (d) a far-western population (Eastern boundary of Banke district to Mahakali River) 45,46.
Table 1
Four different regions within Chure Terai Madhesh Landscape Nepal, area, forest cover and elephant population status. The elephant population was obtained from Ram & Acharya (2020).
SN
|
Region
|
Coverage
(Districts)
|
Total Area (Km2)
|
% Forest cover
|
Elephant population
|
1
|
Eastern
|
Jhapa to Siraha
|
11,116.96
|
31.92
|
Residential: 27–35; ~100 migratory elephants each year from West Bengal, eastern India).
|
2
|
Central
|
Dhanusa to Chitwan
|
8,169.43
|
46.17
|
Residential: 45–53
|
3
|
Western
|
Nawalparasi to Dang
|
8,777.95
|
56.89
|
Migratory: 8–12
|
4
|
Far-Western
|
Banke to Kanchanpur
|
14,391.54
|
46.94
|
Residential: 80–125; ~45 migratory (from Uttarakhand and Uttar Pradesh, India migrated to far western habitats in Nepal )
|
|
Total
|
|
42,455.88
|
44.92
|
|
Derivation of forest cover
We analyzed forest cover change and fragmentation using both the patch and landscape metrics and considered forest fragmentation as habitat fragmentation 47. We categorized forested areas as natural and plantations with a tree canopy cover of more than 10 percent and an area of more than 0.5 ha 48. We used the hybrid classification techniques to combine high-resolution images, medium resolution images, and digitization of topographic maps. First of all, we prepared a forest cover map of the 1930s by digitizing greenwash areas shown on topographical maps prepared by Army Map Service, U.S. Army, Washington, surveyed during 1920–1940 (http://legacy.lib.utexas.edu/maps/ams/india/) at 1:250,000 scale. Due to the unavailability of multi-spectral satellite images of the study area before the 1970s, we relied on the existing topographic maps to obtain forest cover of 1930 49,50
51,52 found 5–10% inherent errors at various stages of land cover change analysis; while using historical data and topographic maps. The inaccuracy of forest cover mapping was minimized by visual interpretation and overlay analysis in the topographic maps. In addition, we resampled all the digital images at a 30-meter resolution to improve the mapping errors. 51,52 reported the reliability of topographical maps to reconstruct forest cover. We also obtained the forest cover map of 1975 by on-screen digitization of Landsat 1 TM level 1 satellite images.
We produced the forest cover maps of 2000 and 2020 from Landsat imagery scenes of respective years (Table 2; Figure. 2). All the Landsat data processing was conducted using the cloud-computing technology in the Google Earth Engine (GEE) platform (https://earthengine. google.org/). The GEE platform carried out a fast analysis using Google’s computing infrastructure 53,54. We used the pre-processed Landsat imagery available through GEE to assess forest cover change across the study area 55.
Table 2
Sources of data used in this study.
SN
|
Data layer
|
Source
|
Spatial resolution (m)
|
Year
|
1
|
Topographic map
|
Army Map Service, U.S. Army, Washington
|
Scale 1:250,000 m
(based on Arial photo)
|
1920–1940
|
2
|
Landsat 1 TM
|
Earth Explorer (USGS)
|
60 m
|
1975–1976
|
3
|
Landsat 5 Surface Reflectance Tier 1
|
GEE dataset (USGS)
|
30 m
|
2000
|
4
|
Landsat 8 Surface Reflectance GEE dataset
|
GEE dataset (USGS)
|
30 m
|
2020
|
5
|
Administrative boundary
|
Department of Survey, Nepal
|
Scale 1:25,000
(Based on Arial photo)
|
1996–1998
|
We used a cloud screening algorithm to remove cloud contaminated pixels from each Landsat image by applying quality assessment (QA) bands for 2000 and 2020. Then, we produced an annual composite by taking the median value from images from the target year 56. We delineated > 1000 reference points for each period 2000 and 2020 respectively. We used supervised machine learning classifiers, i.e., Random Forest (RF), to classify remotely sensed data 57. Random Forest Classifier creates a set of decision trees from a randomly selected subset of the training set and aggregates the votes from different decision trees to classify the image 58. The classified image was downloaded as raster tiff files. The raster was converted into vector polygons and overlaid with high-resolution google earth images of respective years. The final forest cover map was obtained with the highest accuracy by post-processing (validating) the forest polygons through onscreen digitization to match the forests visible in Google images 57
Data analysis
Analysis of forest loss/gain
Forest cover maps of the four different periods of 1930 (before malaria eradication), 1975 (the initial stage of PA system development), 2000 (well-established PA system) & 2020 (current scenario) were post-processed according to FAO forest definition. These layers were analyzed to understand changes in extent and location of forests using a post-classification change detection technique in Arc GIS 10.4.
We have estimated the conversion of forests into the non-forest area on a grid overly basis and grids were generated on the basis of minimum home range size of elephants i.e., 18 km2 45,59. We generated 5× 5 km2 grids for the time series assessment and analyzed spatial distribution trends of forest cover in these grids from 1930–1975, 1975–2000, and 2000–2020 49,60. We computed the forest cover area (distribution of transitions and persistence of forest) of four different periods in each grid using the zonal statistics tool of ArcGIS software 61. Overall, forest cover change was calculated by combining all the grids and calculating the annual deforestation rate (percentage) using a compound-interest-rate formula 62.
where a1 and a2 are the area covered by forest at times t1 and t2. The region wise rate of deforestation was computed and presented.
Modeling Forest fragmentation
We carried out habitat fragmentation analysis in the four regions of CTML (Fig. 1) and measured fragmentation in terms of core, perforated, edge, and patches. We used 30 m cell resolution for fragmentation analysis for four different periods. We used patch analyst 63 to obtain the patch matrix for each region viz. patch density and size (number of patches, mean patch sizes, patch size standard deviation), edge metrics (edge density, mean patch edge), and shape index (mean shape Index, mean perimeter area ratio, mean patch fractal dimension), (Supplementary table S1).
Similarly, Landscape Fragmentation Tool (LFT V2.0, http://clear.uconn. Edu/tools/lft/lft2/) was used to estimate landscape metrics 64. The change of fragmentation during the 1930 to 2020 periods was carried out by cross-tabulating the fragmentation classes. Landscape Fragmentation Tool (LFT) classifies forests at pixel-level into fragmentation classes: core 1, core 2, core 3, perforated, edge, and patch. Core forests are located far from the forest/non-forest boundary and surrounded by other forest areas. We considered the core forest as 100 m distance from the edge 65. The core forests include three different types – 1) Core1: forest patches area < 250 acres (1.012km2), 2) Core 2: Medium core (forest patches area between 250–500 acres (1.01–2.2 km2), and Core 3: large core (Forest patches area > 500 acres (> 2.2 km2) 49. The peripheral forest was further classified into perforated (i) inner edge: forest pixels on the edge of small interior non-forest, and (ii) edge forest or outer edge: pixels that are between forest and large non-forest areas 66.