Assessment of Forest Cover Changes in Vavuniya District, Sri Lanka: Implications for the Establishment of Subnational Forest Reference Emission Level

Abstract

Assessment of forest cover changes is required to establish the forest reference emission level (FREL) at any scale. Due to civil conflict, such assessments have not yet been undertaken in Sri Lanka, especially in the conflict zone. Here, we assessed the forest cover changes in Vavuniya District, Sri Lanka, from 2001 to 2020, using a combination of the Google Earth Engine (GEE) platform and the phenology-based threshold classification (PBTC) method. Landsat 5 TM data for 2001, 2006, and 2010, and Landsat 8 OLI data for 2016 and 2020 were used to classify forest cover by categories, and their related changes could be assessed by four categories, namely dry monsoon forest, open forest, other lands, and water bodies. With an overall average accuracy of 87% and an average kappa coefficient of 0.83, forest cover was estimated at 57.6% of the total land area in 2020. There was an increase of 0.46% per annum for the entire district between 2001 and 2010, but a drastic loss of 0.60% per year was observed between 2010 and 2020. Specifically, the dry monsoon forest lost 0.30%, but open forest gained 3.62% annually over the same period. Loss and gain of forest cover resulted in carbon emissions and removals of 165,306.6 MgCO₂ and 24,064.5 MgCO₂ annually, respectively, over the same period. Our findings could be used to set the baseline trend of deforestation, based on which, a subnational forest reference emission level can be established as an emission benchmark, against which comparisons of carbon emissions following the implementation of REDD+ activities can be made, and result-based payment can be claimed under the Paris Agreement.

1. Introduction

Forest ecosystems act as a natural remedy to combat climate change and provide ecosystem services to meet the basic needs of global human survival [1,2]. Forest cover is an important indicator in understanding the extent of forest loss or gain, and related carbon stocks in the locations in question [3]; however, its assessment is still challenging. The local circumstance of forests and their characteristics can determine the different roles of forests, from being carbon sources or carbon sinks. Therefore, forest cover is a crucial metric for measuring the performance of implemented activities to obtain result-based financial incentives for “Reducing Emissions from Deforestation and forest Degradation (REDD), conservation of forests, sustainable management of forests, and enhancement of forest carbon stocks” (REDD+) in developing countries; i.e., the REDD+ scheme [4]. The United Nations Framework Convention on Climate Change (UNFCCC) adopted this international framework in 2007, with the aim of carbon-based financial incentives to reduce carbon emissions in developing nations. At the 26th Conference of the Parties (COP26) of the UNFCCC in 2021, world leaders under the “Action on Forests and Land Use” theme committed to halt and reverse forest loss and land degradation by 2030 [5]. To achieve this under the Paris Agreement’s Warsaw Framework for the REDD+, parties (annexed and non-annexed countries) need to submit a national action plan, national forest monitoring system (NFMS), forest reference emission level (FREL), and safeguard information system (SIS) in order to benefit from the REDD+ scheme during the Paris Agreement period from 2020 to 2030 [4,6].

Among these four essential elements, the FREL is a benchmark to measure the performance of implemented REDD+ activities under the business-as-usual scenario, and is defined as the amount of carbon emissions from deforestation and forest degradation in a geographical area in the absence of the project’s activities [7,8,9]. To be eligible for the result-based financial incentives, developing counties need to establish a national or subnational FREL to compare future carbon emissions against historical reference levels [10,11]. According to the recent update by the Food and Agriculture Organization of the United Nations (FAO) in 2020, about 80% of the FRELs were nationwide in scope, using the historical average of carbon emissions. To measure the performance as required by the UNFCCC, the subnational FREL can be used as an interim measure [12,13]. However, previous studies suggest that starting FREL at the subnational level would be a more precise procedure for obtaining reliable data and reducing uncertainty than a national-level generalization [14]. Until recently, only four Latin American and two African countries submitted their subnational FRELs, followed by the national FRELs [13].

FREL at national and subnational levels can be developed using retrospective and prospective approaches [15]. The retrospective approach uses the historical emissions to project future emissions [16]. Development of the subnational FREL is required for monitoring the performance of subnational REDD+ activities. Subnational monitoring can also play a key role in countries’ long-term REDD+ and low-emissions development strategies. According to the United Nations Development Programme (UNDP), research conducted in 2009 estimated that 50–80% of mitigation and up to 100% of adaptation to climate change could be influenced by decisions made at the subnational- and local-level implementation of REDD+ [17]. Apart from national-level preparedness initiatives, such as improving data collection, developing FREL and operational frameworks, and evaluating REDD+ policy options, the majority of REDD+ demonstration activities and project investments are currently concentrated at the subnational and local levels. However, the REDD+ activities and project investments often disregard the crucial role of subnational governments and agencies in forest and land use policy to accomplish REDD+ and sustainable development goals. Pointedly, subnational FRELs and monitoring are recognized by the UNFCCC as a provisional solution; these solutions will also be significant as a permanent measure to enable REDD+ to be implemented under different existing governance frameworks [18].

To establish FRELs using a retrospective approach, developing countries are required to have reliable and consistent long-term forest cover data, because such data are important for establishing baseline carbon emissions. Nevertheless, not every developing country has such data. Recent studies have developed subnational and national FRELs with two data points of forest cover data due to limitations in forest inventory [10,16]. Previous studies indicate that more data points of historical forest cover can reduce the bias and uncertainty in the development of FRELs at different levels [19]. Due to the lack of such data, the inaccurate estimation of historical carbon emissions could result in unnecessary biases in predicting future emissions as part of the requirements for implementing the REDD+ and for receiving carbon benefits. Fortunately, as more satellite images and the latest cloud computing technologies become increasingly available, the use of remote sensing approaches for obtaining historical forest cover data make it possible to assess forest cover changes for countries with otherwise limited data coverage, at any scale [20].

Sri Lanka is a small tropical island in the Asia-Pacific region of the Indian ocean, with a diversity of flora and fauna and different types of natural forests [21]. Based on national statistics, the total forest cover was reduced by 14.5% between 1956 and 2010. Fernando et al. [22] reported that forests covered 29.7% of the total land area in 2010. Therefore, determining FRELs is necessary to reduce deforestation and carbon emissions at subnational and national levels. In fact, Sri Lanka submitted their national FREL and removals to the UNFCCC in 2017 as a benchmark for monitoring and measuring the performance in the forest sector [21]. To develop a subnational FREL, two important variables, forest cover (activity data) and emission factor, are needed [10]. Sri Lanka has forest cover data at the national level. However, there is a lack of forest cover data and inconsistencies in national forest inventories at the subnational level [23]. Chokkalingam and Vanniarachchy [24] stated that Sri Lanka has limited information about the precise extent and nature of existing forests, and recent trends in these forests. The northern and eastern provinces of Sri Lanka were severely affected by the Sri Lankan civil war, and the districts belonging to these provinces lack updated forest cover data [25]. Due to the lack of forest cover data, REDD+ implementation at the subnational level is challenging. However, there is a need to evaluate baseline deforestation to obtain financial incentives under the Paris Agreement by estimating and comparing carbon emissions. A recent technological advancement, Google Earth Engine (GEE), is a cloud-based alternative geospatial platform that can be used for acquiring more than 40 years of historical and current global forest cover data using satellite images [26,27,28]. Thus, we used this platform to derive the past and recent Landsat satellite images to obtain the long-term forest cover data required to overcome the aforementioned challenges at the subnational level. Importantly, a subnational FREL can convey the REDD+ activities and benefits into practical reality with less uncertainty. Therefore, our study focused on assessing the forest cover changes in Vavuniya District, Sri Lanka, between 2001 and 2020, using a combination of the phenology-based threshold classification (PBTC) method and multitemporal Landsat collections available in GEE to establish the subnational FREL.

2. Materials and Methods

2.1. Study Area

Our study was carried out in the Vavuniya district, Northern Province, Sri Lanka. This district is categorized as a dry zone (rainfall < 1750 mm) based on the rainfall pattern. Its geographical coordinates are 8°45′0″ N and 80°30′0″ E. This district has four Divisional Secretariat (DS) divisions, namely Vavuniya, Vavuniya North, Vavuniya South, and Vengalacheddikulam (Figure 1). The climate is tropical, with an average temperature of 28 °C. Sri Lanka receives rainfall during two monsoon periods. In this study area, the highest rainfall occurs during the Northeast monsoon (December to February). The dry season, experienced from May to September, is defined as the Southwest monsoon [29]. In summer, the mean temperature is high, while the minimum temperature is recorded in winter. The annual rainfall received in this district is 1400 mm. Surface wind circulations influence the seasonal variation in rainfall of the climatic zones, causing differences in monsoonal patterns [30]. Persistent year-round cloud cover significantly affects the availability of remotely sensed data. In this district, the dry monsoon forest (semi-deciduous vegetation) is the dominant type, with a considerable proportion of the open forest. This is one of the war-affected districts in the Northern Province. The internal civil conflict has badly impacted vast areas of the forest cover in the northern part of the district. This district is one of the most important districts in the Northern Province, and the town of Vavuniya serves as a regional transit hub. The district’s economy is mainly based on agriculture, and has been highly developed over the last decade. However, the rural community still strongly relies on the forest to obtain products and services to enhance their livelihood.

2.2. Satellite Remote Sensing Data and Pre-Processing

We acquired satellite remote sensing data to assess the changes in forest cover from 2001 to 2020. GEE is an effective processing tool of global-scale satellite imagery dating back 40 years [26]. The required data were freely collected from the Landsat collections on the GEE platform [27]. The selected Landsat image collections were available from the Landsat 5 and Landsat 8 satellites. Landsat 5 is the longest operating earth observation satellites among the current Landsat series. It was launched in 1989, and was decommissioned, after its 28-year mission, in 2012. Meanwhile, Landsat 8 was launched in 2013, and the latest data were obtained through the operational land imager (OLI) and thematic infrared sensor [32]. In particular, Landsat 5 thematic mapper (TM) top of atmosphere (TOA) and Landsat 8 OLI TOA reflectance image collections with 30 m moderate spatial resolution were obtained for land cover classification [28,33,34] to estimate the forest cover changes [35]. Landsat TOA reflectance refers to the reflectance measured by a sensor from space flying above the earth’s atmosphere. Landsat TOA reflectance data has significant potential to determine threshold values and perform PBTC of land cover categories with greater accuracy due to its higher reflectance-based EVI values [32]. Due to climatological variations, the dry season is a suitable time to collect Landsat images as there is less cloud cover [28,36] for EVI-based land cover assessments. Thus, we collected Landsat 5 TM data for 2001, 2006, and 2010, while Landsat 8 OLI data were obtained for 2016 and 2020 during the dry season, covering the whole district. Therefore, the entire study area was covered by two Landsat 5 collection pathways and three paths of Landsat 8 collections. In total, we gathered 150 images from the Landsat 5 TM and Landsat 8 OLI collections.

Image pre-processing is an important step to produce good-quality maps before performing classification. Initially, the research area’s boundary shapefile was imported as an asset to perform this pre-processing. Then, GEE’s image collection feature was applied to gather the Landsat TM TOA and OLI TOA reflectance collections, and the filter function was employed to obtain the dry season collections. Next, GEE algorithms were applied on Landsat composites to reduce the imagery to median pixel values [37]. Subsequently, the cloud mask function was executed to reduce the cloud cover by less than 50% on all Landsat collections in GEE [38,39,40] to reduce the probability of acquiring bad images [33,41,42]. Subsequently, we applied the image compositing and reducer functions to choose the recent pixel in the Landsat collections. The median reducer function estimated the median value for all collections. This function helps to remove the high-value and low-value clouds and cloud shadows. Finally, the clip function was applied to remove the unnecessary portions of outsized images other than the study area.

2.3. Methodology for Phenology-Based Threshold Classification

2.3.1. Enhanced Vegetation Index (EVI)

Subsequently, our study applied the EVI function on the median Landsat TOA reflectance collections (Landsat 5 TM and Landsat 8 OLI) obtained for 2001, 2006, 2010, 2016, and 2020 [43,44]. This index resembles the normalized difference vegetation index (NDVI); however, it can overcome the saturation problems in dense vegetation cover [45], and can minimize the effects of atmosphere, soil, and solar irradiance [46] due to the use of the blue reflectance region. In addition, research by Gao and Mas [47] showed that the EVI is a good indicator for phenology classification of land cover types. Consequently, Landsat EVI time series data were used to classify the selected land cover categories [45]. The EVI calculation formula is shown below [48,49]:

$$\mathrm{EVI}=\mathrm{G}\times \frac{\mathrm{NIR}-\mathrm{R}}{\left(\mathrm{NIR}+\mathrm{C}1\mathrm{R}-\mathrm{C}2\mathrm{B}+\mathrm{L}\right)}$$

(1)

where G = 2.5 (gain factor), L = value to adjust for canopy background, C1 = 6, and C2 = 7.5 (aerosol correction factors) [50]. NIR, near-infrared; R, red band; B, blue band.

For Landsat 8: Band 5—NIR, Band 4—R, and Band 2—B;

For Landsat 5: Band 4—NIR, Band 3—R, and Band 1—B.

2.3.2. Enhanced Vegetation Index Threshold Values to Classify Land Cover Categories Using PBTC

The PBTC method plays a significant role in land cover mapping because it can differentiate the phenological behaviors of plants throughout the year [28]. This method is sensitive to climatological variation, but it is a powerful tool for acquiring information regarding land cover types using EVI threshold values. Recently, the application of PBTC in forest cover mapping has drawn significant attention. This method has shown a high level of accuracy and higher reliability for classifying and mapping croplands, forest cover, and land cover to assess changes and estimate forest carbon stocks during recent decades, due to many technological challenges being resolved in the GEE [28,35,36].

This study applied the PBTC method to classify four land cover categories, namely, dry monsoon forest, open forest, other lands, and water bodies (

Table 1. Description of land cover categories.

Land Cover Categories	Description
Dry monsoon forest	Altitude 0–500 m, rainfall 1250–1900 mm, areas influenced by the summer monsoon, with less rainfall, 4–5 dry months, and canopy cover above 40% [22]
Open forest	Land spanning more than 5 hectares, tree canopy cover between 10% and 40% [22]
Other lands	Croplands, barren lands, shrubs, and settlements
Water bodies	Areas with fresh water

and Figure 2). However, our research mainly focused on forest cover classes (dry monsoon and open forest). We modified a PBTC model developed by Venkatappa et al. [51] to classify the land cover as per the climatological characteristics and the phenological behavior of vegetation of our study area. Landsat collections were obtained during the dry season in Sri Lanka, from May to September, to differentiate each land cover category from other categories. The collection of images for the dry season was performed because the dry season is the period during which plants behave differently; thus, the satellite technology can capture the EVI, with less non-overlapping values, of the selected land cover categories, thereby avoiding the misclassification of the EVI if wet season Landsat images are collected. The land cover categories during the wet season can appear as green, and yield overlapped EVI threshold values. Note that our study identified new EVI minimum and maximum threshold values for each land cover category during the dry season in Sri Lanka.

Table 2. Vegetation thresholds to classify the land cover categories.

Land Cover Categories	EVI Threshold Values
	Landsat 5 TM		Landsat 8 OLI
	Minimum	Maximum	Minimum	Maximum
Dry monsoon forest	0.5027	0.776	0.650	0.882
Open forest	0.4845	0.501	0.631	0.648
Other lands	0.110	0.421	0.210	0.556
Water bodies	−0.13	0.105	−0.13	0.209

Land cover categories 1 and 2 are considered as forest and the remaining are non-forest.

shows the EVI threshold values for land cover categories. After the image pre-processing stage mentioned in Section 2.2, all image collections from 2001 to 2020 were processed using the EVI function for the median Landsat collections. Thereafter, to perform the PBTC, the identified minimum and maximum threshold values were applied for each land cover (Table 2). Land cover classification and spatial mapping were performed by executing the PBTC function in GEE using the minimum and maximum EVI threshold values [28,51]. Different colors were assigned for the different land categories in the layer palette to visualize the land cover categories in GEE. Finally, our research yielded five land cover maps from which to analyze forest cover changes.

2.4. Estimation of Forest Cover Changes

Multitemporal datasets were used to distinguish areas of land cover changes between remote sensing imaging dates [52,53]. After an accuracy assessment, we applied the following formulas to assess the annual, and rate of, forest cover changes [54,55] during 2001–2006, 2006–2010, 2010–2016, 2016–2020, and 2001–2020. In addition to these time frames, the rate of change in forest cover was assessed from 2001–2010 and 2010–2020 based on the pattern of change. The annual rates of forest cover changes were calculated by comparing the area of forest cover at different time points.

$$\mathrm{Forest} \mathrm{cover} \mathrm{change} \mathrm{in} \%=\frac{{\mathrm{FAt}}_1-{\mathrm{FAt}}_0 }{{\mathrm{FAt}}_0}\times 100$$

(2)

$${\mathrm{FCC}}_{\mathrm{n}}=\frac{{\mathrm{FAt}}_1-{\mathrm{FAt}}_0}{{\mathrm{t}}_1-{\mathrm{t}}_0}$$

(3)

$${\mathrm{R}}_{\mathrm{n}}=\frac{{\mathrm{FCC}}_{\mathrm{n}}}{{\mathrm{FAt}}_0}\times 100$$

(4)

where FCC_n = annual forest cover change (ha/y) in nth period, FAt₁ and FAt₀ = forest area (ha) at final year (t₁) and initial year (t₀), R_n = annual rate of forest cover change in each period n (%/y)

2.5. Accuracy Assessment and Validation

Accuracy assessment reveals discrepancies between the classified maps and reference data [56,57]. Therefore, the accuracy assessment and validation were considered as a foundation [58]. In this accuracy assessment, the Google Earth Pro platform was used to provide reference data to validate the land cover maps [59]. This platform regularly updates when new data are available. Google Earth’s timelapse function was accessed to obtain zoomable historical images from 2001 to 2020 to visualize and monitor land cover changes using GEE [60]. In this study, we acquired, from Google Earth, very high-resolution time series images to validate the accuracy of land cover maps by applying standard protocols recommended by Olofsson et al. [61]. The results must be quantitative to determine their classification accuracy. Therefore, a stratified random sampling technique was applied to validate the maps [51]. We considered each land cover category as a stratum for this sampling.

Consequently, the generated maps were exported to Google Drive using the export function to carry out the accuracy assessment. The ArcGIS Pro is a professional single-desktop geographic information system (GIS) application with a feature-packed software that has replaced ArcMap software [62]. We used the generated PBTC maps in the latest ArcGIS Pro software to transform the land cover labels into accuracy points. The sampling points were created using the stratified random sampling technique tool “create accuracy assessment points” in ArcGIS Pro. A total of 1000 accuracy assessment points were generated, and we then converted the sampling point files into a “keyhole markup language” (kml) file, which visualizes the geographic data of created accuracy assessment points on the Google Earth Pro. Later, the converted kml files were imported to Google Earth Pro for validation. Subsequently, sampling points were visually analyzed using the Google Earth “time slider” tool to compare and validate the maps with reference data. An error matrix is a contingency table used to show the frequency of discrepancies between classified maps and reference data. The classification accuracy of remotely sensed images is widely assessed using the error matrix [63]. Finally, the following equations were used to calculate producer accuracy, user accuracy, overall accuracy, and the kappa coefficient [63,64] using the error matrix method in Microsoft Excel.

$$\mathrm{Producer} \mathrm{accuracy}=\frac{\mathrm{Number} \mathrm{of} \mathrm{accurately} \mathrm{identified} \mathrm{points} \mathrm{in} \mathrm{an} \mathrm{individual} \mathrm{land} \mathrm{cover} \mathrm{category}}{\mathrm{Total} \mathrm{reference} \mathrm{points} \mathrm{of} \mathrm{respective} \mathrm{land} \mathrm{cover} \mathrm{category}}\times 100$$

(5)

$$\mathrm{User} \mathrm{accuracy}=\frac{\mathrm{Number} \mathrm{of} \mathrm{accurately} \mathrm{identified} \mathrm{points} \mathrm{in} \mathrm{an} \mathrm{individual} \mathrm{land} \mathrm{cover} \mathrm{category}}{\mathrm{Total} \mathrm{number} \mathrm{of} \mathrm{classified} \mathrm{reference} \mathrm{points}}\times 100$$

(6)

$$\mathrm{Overall} \mathrm{accuracy}=\frac{\mathrm{Total} \mathrm{number} \mathrm{of} \mathrm{correctly} \mathrm{classified} \mathrm{land} \mathrm{cover} \mathrm{categories} \mathrm{reference} \mathrm{points}}{\mathrm{Total} \mathrm{number} \mathrm{of} \mathrm{reference} \mathrm{points}}\times 100$$

(7)

$$\mathrm{Kappa} \left(\mathrm{K}\right)=\frac{\left(\mathrm{Overall} \mathrm{accuracy}-\mathrm{random} \mathrm{accuracy}\right)}{\left(100-\mathrm{random} \mathrm{accuracy}\right)}$$

(8)

The kappa coefficient is a measurement used to determine the agreement between classification accuracy and the reference data [49,65]. The parameter reflects the difference between actual agreement and the agreement expected by chance. Values of K range from 0 to 1, where 1 denotes perfect agreement, while 0 indicates poor or no agreement in the classification [66]. This indicates how well the categorized map corresponds to the reference data. When compared to a random classification, the kappa coefficient reflects the proportionate reduction in classification error.

In addition, ground truth data were used to validate the PBTC map for 2020. In total, we collected 350 ground truth field points (dry monsoon forest: 139 points; open forest: 68 points; other lands: 84 points; and water bodies: 57 points) randomly, during the dry season in 2021, using the Global Positioning System (GPS) and SW maps. The SW map is a freely available application on Android for recoding data through GPS [67]. The collected points were imported into ArcMap, and the spatial analyst tool was used to extract the values to points on classified maps to compute the error matrix.

3. Results and Discussion

3.1. Land Cover Pattern of the District

We assessed land cover changes using multispectral images from Landsat 5 TM and 8 OLI collections with 30 m moderate spatial resolution from 2001 to 2020. Figure 3 represents the land cover maps for different time frames, and indicates that the classified maps consist of four land cover categories; dry monsoon forests, open forests, other lands, and water bodies. The produced maps are a noteworthy finding in terms of understanding their dynamic pattern for the different periods. Pointedly, in all PBTC maps, dry monsoon forest was the predominant land cover category, followed by other lands, open forest, and water bodies. Over the study period, the extent of other lands increased, while that of the dry monsoon forest declined. The open forest fluctuated, and other lands developed along the road networks. The maps show a considerable degree of variation in the area coverage of these land cover categories during the study period (Figure 3).

3.2. Forest Cover at District Level

The changes in forest area over the 19 years studied are shown in

Table 3. Forest cover areas and their percentages for the total district.

Forest Cover Class	Forest Cover Area (ha) between 2001 and 2020
	2001	%	2006	%	2010	%	2016	%	2020	%
Dry monsoon forest	110,990	55.4	109,560	54.7	109,295	54.5	106,621	53.2	104,627	52.2
Open forest	6338	3.2	10,617	5.3	13,453	6.7	8145	4.06	10,702	5.3
Total forest area	117,327		120,177		122,748		114,766		115,329
District area	200,390 ha

. The total district area was around 200,390 ha; however, more than half of the total land area was occupied by natural forests. In 2020, 52.2% and 5.3% of land area were occupied by dry monsoon and open forests, respectively. The maps describe the spatiotemporal distribution, and its changes, of the district’s forest cover (Figure 3). Overall, the dry monsoon forest was the most dominant and widely distributed forest type in the study area. In contrast, open forest covered a lower proportion of the land area. The open forest was scattered within the dry monsoon forest, and fluctuated throughout the study, as a result of drivers of deforestation and forest degradation looming towards the dry monsoon forest. Specifically, in recent decades, other lands have started to occupy areas that were once dry monsoon forest. Although dry monsoon forests are enriched in the northern part of this district, the forest area in this region declined in 2020, while the southern region gained forest cover. Remarkably, the highest percentage of total district land area was covered by dry monsoon forest.

3.2.1. Trend in District Forest Cover

Figure 4 represents the overall trend in forest cover concerning the total district area from 2001 to 2020. The data show an overall increasing tendency between 2001 and 2010; however, a significant reduction was observed after 2010 [68]. About 58.5% of the total geographic area of Vavuniya District was sheltered by forest in 2001. In 2006, this started to gradually increase, peaking at 61.3% by 2010, and then rapidly reduced to 57.6% in 2020. However, a sudden decline was noted between 2010 to 2016, which slightly improved thereafter. Our findings indicate that forest loss occurred at a low rate over the study period; however, the loss accelerated after 2010.

In contrast to these results, the latest global forest resource assessment from the FAO [69] reported that 178 million ha of the world’s forest has been lost since 1990. Worldwide forest cover declined to 3999 million ha from 4128 million ha between the years 1990 and 2015. Thus, about 420 million ha of forest has vanished worldwide due to deforestation since 1990, with the yearly rate of deforestation estimated at 10 million ha between 2015 and 2020, down from 12 million ha in 2010–2015. In South Asia, 29.62% of forest has been lost in the eight decades since 1930, with annual net deforestation rates of 0.06% and 0.04% during 1995–2005 and 2005–2014, respectively [70]. Gunawardane et al. stated that forest cover in five districts in the Northern Province, including Vavuniya, increased by 3% between 1999 and 2008 due to limited human access and disturbance to the forests. From national statistics, the total forest cover in Sri Lanka reduced by 14.5% between 1956 and 2010, with forests covering 29.7% of the total land area in 2010 [23]. From 2005 to 2014, the annual rate of deforestation was 0.01%, while a high rate of deforestation was recorded between 1976 and 1994 [71]. The deforestation hotspots were discovered in Sri Lanka’s southeast and northern regions [72]. Deforestation is a major driver of forest cover loss in the tropical region [72,73]; our results indicate that this is also a common problem in Sri Lanka. Therefore, smart conservation and management strategies are required to safeguard remaining resources. However, forest cover change assessments in Sri Lanka are complex due to land use transitions that have occurred throughout the country’s history [74,75].

3.2.2. Annual Rates of District Forest Cover Change

We divide the study period into two time frames based on pattern variation (Figure 4), 2001–2010 and 2010–2020, to explain the rate of forest cover change. During the first time frame, the district’s forest area increased by 4.62% (5420 ha). However, a substantial rapid decline of 6.04% (7419 ha) was observed between 2010 and 2020. As shown in Figure 4, considerable deforestation occurred during this later period. Specifically, during 2010–2016, other lands increased in cover (Figure 3) due to rehabilitation, housing settlements, agricultural expansion, and development activities implicated soon after the internal civil war in 2009. It is important to note the forest loss in this period. Fernando et al. [23] noted rapid economic development in the northern part of Sri Lanka in 2009. After 2016, post-war resettlements and advancements in rural areas and infrastructure slightly increased the forest area [23,76].

Comparisons based on spatial extent revealed significant changes observed during 2001–2010 and 2010–2020. Remarkably, a striking turning point was noted in 2010. Importantly, forest cover expanded in the initial time frame, and then rapidly declined during the recent decade. Previous research in the dry zone conducted by Ranagalage et al. [68] demonstrated that forest cover in the Vavuniya District declined by 8.7% in 2019, from 122,325 ha in 2010. In addition, all districts in the Northern Province resembled an increasing and decreasing pattern of forest cover change from 1999 to 2010 and 2010 to 2019, respectively. We also observed a similar result in this study site. From our statistics, the pattern of the district’s forest cover trend was consistent with this previous study. Moreover, the annual rates of forest cover changes were estimated (

Table 4. Annual rate of forest cover changes.

Annual Rates of Forest Cover Changes
Time Frame	Total Forest		Dry Monsoon Forest		Open Forest
	ha/year	%/years	ha/year	%/years	ha/year	%/years
2001–2006	569.94	0.49	−285.88	−0.26	855.81	13.50
2006–2010	642.68	0.53	−66.45	−0.06	709.12	6.68
2010–2016	−1330.27	−1.08	−445.54	−0.41	−884.74	−6.58
2016–2020	140.69	0.12	−498.67	−0.47	639.36	7.85
2001–2020	−105.18	−0.09	−334.90	−0.30	229.71	3.62

) in the period from 2001–2006 (0.49%), followed by 2006–2010 (0.53%), 2010–2016 (−1.08%), and 2016–2020 (0.12%). In the last 19 years, the annual change in the district’s forest cover was −105.18 ha/year, with an annual rate of change of −0.09% (i.e., forest decreased at a rate of 0.09%). The annual deforestation rate was high from 2010–2016. Although forest cover dynamics in Sri Lanka’s dry zone are complex, non-deterministic, and non-uniform, with decadal changes in deforestation and afforestation. In the dry zone, human resettlement, agricultural encroachment, and infrastructure development have induced forest cover loss in recent decades [23,77].

3.2.3. Forest Cover Changes by Type

Figure 5 illustrates the trend in dry monsoon forest and open forest cover change in different time frames. Our findings reveal that, in 2001, 55.4% (110,990 ha) and 3.2% (6338 ha) of the total land area was covered by dry monsoon and open forest, respectively. In 2020, the coverage of dry monsoons decreased to 52.2% (104,627 ha), and that of open forests increased to 5.3% (10,702 ha). Statistical analysis shows a considerable scale of negative change in dry monsoon forest compared to open forest. During the period 2001–2020, the dry monsoon decreased by 6.13%.

3.2.4. Annual Rate of Dry Monsoon and Open Forest Cover Change

We estimated the annual rates of dry monsoon and open forest changes over different time frames (Table 4). The annual deforestation rates of dry monsoon forest were 0.26% per annum from 2001 to 2006, and −0.06% per annum from 2006 to 2010. We computed that the highest annual rate of forest loss was from 2010–2016 (0.41%) and 2016–2020 (0.47%), compared to the previous decade. Between 2010 and 2020, it declined by −0.88% per year. However, from 2001–2020, the annual change in dry monsoon forest area was −334.90 ha/year, with a rate of change of −0.30% per year.

In contrast, the open forest cover increased by 28.95% over the study period. The annual rates of open forest cover changes were 13.50%, 6.68%, −6.58%, and 7.85% for the periods 2001–2006, 2006–2010, 2010–2016, and 2016–2020, respectively. In addition, open forest showed a rising trend of 229.71 ha/year, with an annual rate of change of 3.62% between 2001 and 2020. Overall, the open forest area showed an increasing trend, with fluctuations and dry monsoon forest area, followed by a continuous decreasing trend. Thus, our analysis reveals the conversion of dense forests to open forests and other lands, when deforestation accelerated, at different scales (Figure 3).

The dry monsoon forest plays a crucial role in long-term forest cover dynamics. Interestingly, Sri Lanka has the larger extent of monsoon forests [78], with dry monsoon forests covering over 55.3% of the total land area in 2010 [22]. Noticeably, this forest is widely distributed throughout the dry zones, and can offer a diverse range of products and environmental services [79]. Meanwhile, the open forest can be subjected to land encroachment for various purposes [68]. In the dry zone, illegal extraction of timber, firewood, non-timber product collection, cattle foraging, and mining for gravel and granite, are localized drivers of deforestation [80]. If major forest reduction continues, the area in the dry zone could be further subjected to negative consequences of prolonged drought due to lack of rainfall [81], soil erosion, extreme heat, agriculture loss, human and wildlife conflict, forest fires, and climate shifts [79,82,83,84]. In the Nationally Determined Contributions (NDCs) report in 2021, Sri Lanka declared to increase forest cover to 32% by 2030. To support this, our analysis indicated that open forest in Vavuniya District contributes about 3.62% annually towards achieving the NDCs target. In developing nations such as Sri Lanka, limitations to obtaining forest inventory data and lack of historical records are the fundamental barriers to determining multidecade changes in forest area. The satellite-based inventory using GEE can provide long-term information on forest cover for the REDD+ implementation, FREL development, national-level decision making, policy re-enforcement, and monitoring and planning for forestry and other sectors, at the subnational and local levels [85,86].

3.3. Forest Cover Changes at the Divisional Secretariat (DS) Level

We assessed the proportion of forest area with respect to the total land area of each DS division. Thus, this study identified the overall forest cover reduction and enhancement on different DS divisions of the district. In 2001, forest covered about 42.5%, 51.67%, and 40.34% of the Vavuniya, Vavuniya South, and Vengalacheddikulam divisions’ land area, respectively, increasing to 45.50%, 56.37%, and 57.83% in 2020, respectively. In contrast, 83.78% of the Vavuniya North division was covered with forest in 2001. However, this declined to 67.69% over the same period. Notably, an earlier study by Pathmanandakumar and Thasarathan [87] mentioned that forest cover change in the Vavuniya DS division increased by 8% in 2017 from 39% in 1997. However, we found a slight forest cover increase in this division’s land area. Noticeably, the categorization of our land cover classes differs from that of this previous study. We categorized open and dry monsoon forests, while grasslands and shrubs were included in the land cover class referred to as other lands. Figure 6 presents the annual changes in forest cover by DS division. We observe forest area expansion in the Vavuniya, Vavuniya South, and Vengalacheddikulam divisions. In contrast, the Vavuniya North division underwent a severe reduction of 19.2%. Between 2001 and 2020, the annual rates of forest cover gain were 0.42%, 0.48%, and 2.28% in Vavuniya, Vavuniya South, and Vengalacheddikulam DS divisions, respectively. Meanwhile, the deforestation in Vavuniya North division was 1.01% per annum. Significantly, we observed a dramatic decreasing trend in forest cover in Vavuniya North, while the Vengalacheddikulam division showed a gradual increase because of the expansion and reduction in other land areas.

Annual Rate of Dry Monsoon and Open Forest Changes by DS Division

We further examined the dry monsoon forest cover gain and loss in DS (

Table 5. Annual rates of dry monsoon and open forest cover changes, from 2001–2020, by DS division.

DS Division	Dry Monsoon Forest	Open Forest
Vavuniya	0.18	3.36
Vavuniya South	0.003	3.49
Vengalacheddikulam	0.02	2.43
Vavuniya North	−0.01	4.95

). Overall, the Vavuniya North division contained the highest dry monsoon forest area in the entire district, Vavuniya division had the second highest forest area, followed by the Vengalacheddikulam and Vavuniya South divisions. Among the four DS divisions, about 14,168 ha (23%) of dry monsoon forest was significantly deforested in Vavuniya North due to post-war rehabilitation and ongoing development projects. This would have produced a great impact on carbon stocks. Conversely, a substantial rise in forest cover by 43% was observed in Vengalacheddikulam between 2001 and 2020, while Vavuniya and Vavuniya South divisions showed a slight rising trend of 3.5% and 5.7%, respectively.

Between 2001 and 2020, the annual forest cover change rates were 0.18%, 0.003%, 0.02%, and −0.01% in Vavuniya, Vavuniya South, Vengalacheddikulam, and Vavuniya North divisions, respectively. Conversely, a growing trend in open forest areas was noted in all DS divisions. Annual rates of open forest cover gain in DS divisions were 3.36% in the Vavuniya division, 3.49% in Vavuniya South, 2.43% in Vengalacheddikulam, and 4.95% in the Vavuniya North division. We observed the highest percentage of open forest gain per annum in the Vavuniya North division.

In Sri Lanka, a civil war begun in 1989 between the Government of Sri Lanka (GOSL) and the Liberation Tigers of Tamil Eelam (LTTE); this war ended in 2009. The Northern and Eastern Provinces of Sri Lanka were hard hit by this conflict. The several districts in the Northern Province, namely, Kilinochchi, Mullaitivu, Mannar, and Vavuniya districts, collectively called “Vanni”, were captured and controlled by the LTTE [88]. Cultural, economic, and environmental impacts were caused by this civil war at national and international scales [88]. Remarkably, Chardran [89] mentioned that this district suffered several repetitive displacements throughout the conflict, until 2009. The district was divided into two parts, namely, cleared (under the control of Sri Lankan government) and uncleared area (the civil administration was fully controlled by the Liberation Tigers of Tamil Eelam (LTTE)). These two areas were separated by the Omanthai Sri Lankan Army checkpoint [90]. The uncleared area was mainly the entire area of the Vavuniya North division. The whole of the Vavuniya North division, and parts of the Vavuniya and Vengalacheddikulam divisions, were affected during the conflict, and resettlement and rehabilitation programs were implemented in this district. However, during the past decades, large portions of the forest cover in the northern part of the country had been severely affected by the civil war and post-war consequences. In addition, Sri Lankan forest cover is fragmented; therefore, maintenance and management are challenging. Furthermore, after the end of the civil war, regions in the north and east have been rapidly altering [91]. Previous studies identified that the underlying deforestation drivers for land use and land cover change were population growth and high demand for urban land, urbanization, development projects, and rehabilitation during the post-civil war situation in the northern part of Sri Lanka. This situation has increased agricultural expansion, encroachment of forest lands for human settlements and agriculture, and illegal logging [87,92]. The rate of deforestation and pressure on forests in developing countries are directly and indirectly determined by human population density and the economic crisis [93,94,95,96].

Notably, the number of abandoned lands for agricultural and settlements was high in Vengalacheddikulam, followed by Vavuniya, Vavuniya North, and Vavuniya South divisions; due to the previous conflict, landowners migrated to foreign countries, water scarcity, and threat from wild animals [97]. The abandoned lands became forest areas due to lack of accessibility. As a result, forest cover may increase in Vengalacheddikulam and Vavuniya divisions. Based on the Sri Lanka United Nations Reducing Emissions from Deforestation and forest Degradation (UN-REDD) report in 2016, internally displaced citizens and refugees were returned to the districts of Mannar, Jaffna, Vavuniya, Mullaitivu, and Kilinochchi in 2015. In certain situations, people have encroached and claimed the ownership of Forest Department lands without following the legal processes. Additionally, misappropriation by local officials and infrastructure development are further reasons for forest destruction in the northeast, with overexploitation of forest resources within reserved areas. A noticeable malfunction was observed, that the forest sector is not included in the Climate Change Secretariat activities, since it is not considered a significant component of Sri Lanka’s climate change mitigation and adaptation plan [25]. This exemplifies the poor coordination between sectors in Sri Lanka to conserve forest resources and mitigate climate change. To overcome the insufficiencies at the national level, subnational-level studies can provide the forest cover data and carbon storage potential of forests, which may develop the interest in forest conservation and sustainable management through the collaboration of various organizations. These types of studies can help to achieve the REDD+ goals and receive result-based financial incentives through the policies and measures in the National REDD+ Investment Framework and Action Plan (NRIFAP).

3.4. Accuracy of Land Cover Maps and Significance of PBTC Mapping

Table 6. Summary of accuracy assessment for dry monsoon and open forests.

Year	Dry Monsoon Forest		Open Forest		Overall Accuracy	Kappa Coefficient
	Produce Accuracy	User Accuracy	Produce Accuracy	User Accuracy
2001	100%	83%	66%	94%	86%	0.81
2006	88%	90%	73%	88%	87%	0.83
2010	86%	86%	80%	87%	86%	0.81
2016	90%	88%	80%	82%	88%	0.83
2020	84%	91%	78%	80%	88%	0.83

shows the producer accuracy, user accuracy, overall accuracy, and kappa statistics of the subsequent land cover maps. The overall accuracy of the generated land cover maps for the years 2001, 2006, 2010, 2016, and 2020 were 86%, 87%, 86%, 88%, and 88%, respectively. The kappa coefficients for the selected years were 0.81, 0.83, 0.81, 0.83, and 0.83, respectively. Our study obtained acceptable accuracy and average kappa values, with strong agreement, above 80% [72,98]. The average overall accuracy of the PBTC maps was 87%, while the average kappa coefficient was 0.82. The overall accuracy and kappa coefficient for the 2020 map using ground truth points were 86% and 0.80, respectively. The accuracy assessment findings revealed that there was an almost acceptable agreement between phenology-based land cover maps and reference data. The overall accuracy and kappa coefficient of maps showed higher reliability and accuracy, as this method primarily relied on the seasonal phenological behavior of vegetation.

We attempted to compare accuracy assessment results of PBTC maps with the findings of earlier studies. For example, a similar study concluded that the average accuracy and kappa coefficient of bamboo mapping in Cambodia from 2015 to 2018 were 89.4% and 0.82, respectively [36]. Similarly, the accuracy of mangrove forest mapping in Brazil from 1985 to 2018 was 87% [99], while another study reported the average accuracy and kappa coefficient of forest cover mapping in Ethiopia during 1991–2019 were 82.2% and 0.89 [2]. Another study also yielded an accuracy range from 86.6% to 89.80% for the long-term deforestation assessment in Sri Lanka between 1976 and 2005 [71]. Remarkably, it can be concluded that the classification performance was highly accurate and in line with previous studies [2,66,100]. Therefore, these results are valid for assessing forest cover changes, and are considered as a best-fit model for our study.

The accuracy assessment for land cover classification can be validated by different methods, such as cross-validation approaches and bootstrapping [101]. In the cross-validation, the independent references were used to validate the generated maps. Different training and validation partitions of the data are used in different cross-validation procedures. The final accuracy results are then provided as the average or best fit of all iterations or folds in the cross-validation approach [101]. For example, a previous study from the eastern part of South Africa investigated the accuracy of classification methods in support vector machines and random forests using bootstrapping and k-fold cross-validation methods. The support vector machine with k-fold cross-validation yielded producer and user accuracies above 80% [102]. However, the k-fold cross-validation and bootstrapping methods have many merits and demerits in terms of sampling size, time consumption, linearity, and randomness. In our study, we classified the land cover categories using EVI threshold values, and then validated them using reference data. The selection of the validation approach depends on the objectives and methods applied in a particular study.

The present study confirms that the dry season is more suitable for land cover classification owing to less cloud cover and unique EVI values for each land cover class; these key characteristics can improve accuracy by avoiding overlap among land cover categories. In contrast, the rainy season or wet season is not suitable due to the possibility of high cloud cover. All land cover categories can appear green during the wet season. Therefore, it is difficult to distinguish one land cover type from another. Due to this, the overlapping of EVI values can result in misclassification among the land cover categories. The possibility of misclassification can be reduced using high-resolution satellite images, field-based reference points, and drone technology. A similar study also mapped Cambodia’s 12 land cover categories using the PBTC method during the leaf shedding season [36]. Therefore, phenology-based classification can be the most prominent method for land cover classification. This method could help clearly differentiate land cover classes with minimum errors to obtain improved results. Our study delivers significantly better results due to the selection of a suitable time for collecting Landsat images. Additionally, the time consumption for producing maps in GEE was shorter than for other classification methods. Our findings clarify the importance of the PBTC method for land cover mapping, and confirm that this is a good and sensitive method based on seasonal variations.

However, open forest classification is quite challenging, and is a limitation of our study, because open forest coverage was smaller and mixed with croplands. Furthermore, there is also the possibility of misclassification due to soil reflectivity. As we obtained the Landsat collection during the dry season, the temperature was high. This temperature increment can cause the shedding of leaves and drying out of vegetation. As a result, the soil becomes exposed to sunlight, and this may affect the EVI thresholds, leading to misclassification [33]. In dry zones, the distinction between trees and shrubs is critical for land use and land cover classification [103,104]. In particular, climatological factors play a significant role in determining the phenological cycles of vegetation throughout the year [104,105]. Our analysis highlights the fact that increases in spatial resolution from moderate to high can improve the classification performance and differentiation of permanent trees from seasonal vegetation.

GEE, like most commercial image processing tools, provides user help in the form of readily available online documentation and tutorials, with a discussion forum to learn, clarify user questions, and share experiences. Unlike most commercial image processing applications, GEE does not require any special hardware. It does, however, necessitate a constant internet connection, which may not always be available [106]. Despite this, remote sensing has become more accessible through this open-access platform. Remarkably, GEE has made it feasible to interpret large amounts of satellite images quickly and reliably, for various purposes, including forest cover mapping and assessment.

Although GEE is a cloud-computing platform, it has a wide range of applications incorporating the monitoring and management of natural resources. However, this platform is highly beneficial for developing countries to collect required historical data where data are lacking. The combination of GEE and PBTC can be used to perform meaningful spatial data analysis. It has many advanced features to process big data for spatial assessments. The forest cover change assessment can be used to estimate the carbon stocks, baseline carbon emissions, and removals needed for measuring and monitoring the performance of REDD+ and obtaining result-based financial incentives. Furthermore, estimating forest cover change can be used to better understand the proximate and underlying drivers of deforestation and forest degradation at the district level. The trends in forest cover changes found in this study can be utilized to quantify baseline forest loss and carbon stocks, and compare carbon emissions with project actions, and are also beneficial for conserving and sustainably managing forests at various levels. Our findings can greatly help existing policies and measures (the national forest policy, national climate change policy, forest planting and restoration, and NRIFAP) to achieve their goals to obtain carbon-based incentives with the reality of subnational level consideration [107].

3.5. Carbon Loss or Gain Due to Forest Cover Change

Using an average carbon stock (aboveground, belowground, and litter were incorporated in this assessment) of 134.6 MgC ha⁻¹ for dry monsoon forest and 28.6 MgC ha⁻¹ for open forest [22,108], total carbon stocks in dry monsoon forest were estimated at 14,941,360.9 MgC in 2001, but declined to 14,084,772.1 MgC in 2020, representing a loss of 45,083.6 MgC/year, or about 165,306.6 MgCO₂/year (=45,083.6 × (44/12), the molecular weight ratio of CO₂ to carbon) (

Table 7. Changes in forest carbon stocks by forest type between 2001–2020.

Forest Cover Class	Total Forest Carbon Stocks (MgC)					Loss or Gain (2001–2020) (MgCO2)
						Total	Annual
	2001	2006	2010	2016	2020
Dry monsoon	14,941,360.9	14,748,936.9	14,713,157.3	14,353,291.8	14,084,772.1	+3,140,827.7	+165,306.6
Open forest	181,065.3	303,320.1	384,359.6	232,696.1	305,763.0	−457,224.9	−24,064.5
Total	15,122,426.2	15,052,257.0	15,097,516.9	14,585,987.9	14,390,535.1	+2,683,600.8	+141,242.1

Note: Plus (+) represents carbon loss or emissions and minus (−) represents gain or removals.

) . The dry monsoon forest occupied a higher area in land extent than open forest, and this forest type is a major carbon reservoir in this district. Surprisingly, the carbon stocks stored in dry monsoon forest gradually decreased, with an annual reduction of 5.7%. Our findings clearly show that this forest type becomes a carbon source. There was a more than two-fold reduction in carbon stocks in dry monsoon forest and discharged carbon emissions between 2010 and 2020, compared to the loss in carbon stocks from 2001 to 2010. If a retrospective approach is adopted for establishing the subnational FREL, the FREL for Vavuniya District is equivalent to 165,306.6 MgCO₂/year. Between 2001 and 2020, carbon stocks in dry monsoon forest showed a continuous decline due to forest cover loss. This result shows a need for appropriate conservation and management strategies through the implementation of REDD+ to reduce the loss in carbon stocks and emissions due to the changes in dry monsoon forest cover. These strategies under the REDD+ can yield reductions in carbon emissions, and could also contribute towards developing countries receiving economic incentives under the element of reducing emissions from deforestation as part of the REDD+ scheme.

In contrast, total carbon stocks in open forest were estimated at 181,065.3 MgC in 2001, but increased to 305,763.0 MgC in 2020, representing a gain of 6563.0 MgC/year, or sequestering (removal) about 24,064.5 MgCO₂ annually. This gain in carbon removal could be eligible for financing support under the element of enhancement of forest carbon stocks of the REDD+ scheme of the UNFCCC. In total, open forest lost about 731,891.1 MgC (0.3%/year), causing net emissions of 2,683,600.8 MgCO₂ between 2001 and 2020, or about 141,242.1 MgCO₂/year (Table 7). Overall, the district’s forests were a source of carbon emissions. The net emissions and annual emissions were about 6.9% and 0.36% of Sri Lanka’s total emissions in 2018, reported in the climate analysis indicators tool [109], indicating that reducing emissions from deforestation and enhancing forest carbon stocks through forest restoration could have an important role in achieving the national emission reduction target in Sri Lanka.

3.6. EVI Thresholds and their Impact on Land Cover Assessment and Carbon Stocks

In addition to the EVI thresholds in Table 2, we applied three different EVI threshold values (Model 1 (M1), Model 2 (M2), and Model 3 (M3)) for each land cover category, except water bodies (

Table 8. M1, M2, and M3 EVI threshold values.

Land Cover Categories	EVI Threshold Values
	Landsat 5 TM			Landsat 8 OLI
	M1	M2	M3	M1	M2	M3
	Min–Max	Min–Max	Min–Max	Min–Max	Min–Max	Min–Max
Dry monsoon forest	0.502–0.776	0.503–0776	0.515–0.776	0.648–0.882	0.652–0.882	0.666–0.882
Open forest	0.485–0.493	0.460–0.502	0.435–0.501	0.569–0.615	0.581–0.650	0.615–0.627
Other lands	0.110–0.460	0.110–0.405	0.11–0.421	0.210–0.516	0.210–0.569	0.210–0.461

), to show the impacts of different EVI threshold applications on the land cover assessment. Notably, the overall average accuracy and kappa coefficient of M1, M2, and M3 were less than 80% and 0.80, respectively. However, we reinstate that the land cover maps produced using the Table 2 EVI thresholds yielded overall accuracy and kappa coefficient above 80% and 0.8, considered as the best-fit model for this study. There was a misclassification among the land cover classes due to the different threshold values of M1, M2, and M3. By applying different EVI threshold values, a range of biases in dry monsoon and open forest cover and carbon stocks occur, as shown in Table 8, when compared with our best-fit model (Table 3) results. Similarly, the percentage bias of each model was estimated for both forest covers (Figure 7). The EVI threshold range for open forest was narrower than for dry monsoon forest. Therefore, the selection of EVI range for open forest cover estimation is challenging because the EVI range was comparatively narrow. The highest percentage of bias range was observed in M1, while the lowest percentage was in M3, for open and dry monsoon forest covers, respectively (Figure 7). Additionally, this study estimated that the root mean square errors (RMSEs) of M1, M2, and M3 for open forest cover were 31,738 ha, 15,901 ha, and 14,709 ha, while RMSEs for dry monsoon forest cover were 22,199 ha, 6715 ha, and 929 ha, respectively. These three models incorporating EVI thresholds resulted in higher error, and also did not satisfy the overall accuracy above 80%. Overall, the dry monsoon and open forest covers, and carbon stocks, were either underestimated or overestimated, depending on the chosen EVI values. Accordingly, these biases can lead to a significant change in the degree of carbon emission and removal (

Table 9. Range the estimated forest cover and carbon stocks affected by the chosen M1, M2, and M3 models by comparing with the best-fit model for forest cover and carbon stocks.

Year	Range of Bias in Forest Cover (ha)		Range of Bias in Carbon Stocks (MgC)
	Open forest	Dry monsoon forest	Open forest	Dry monsoon forest
2001	−323 to +2265	−45,889 to +340	−9048.9 to +64,8016.1	−6,178,856.4 to +43,652.3
2006	−540 to +35,258	−9989 to +562	−15,144.0 to +100,8694.6	−1,346,584.8 to +73,497.4
2010	+373 to +44,575	−15,930 to +1382	+11,080.2 to +127,5251.0	−2,146,234.3 to +188,154.3
2016	+15,989 to +20,785	−9328 to +1420	+457,535.0 to +594,684.0	−1,257,540.4 to +189,098.3
2020	+20,863 to +29,346	−11,571 to +1574	+596,984.7 to +839,621.5	−1,559,399.3 to +209,824.0

Note: Minus (−) represents an underestimation and plus (+) represents an overestimation in forest cover and carbon stocks.

).

The following previous studies found similar results to ours. For example, the variability in EVI thresholds underestimated the soybean cultivation area by 11.4%, while areas of rice field and silviculture were overestimated by 12.8% and 14.3%, respectively [110]. Similarly, three EVI values (0.400, 0.450, and 0.500) used for classifying forests in Guyana, produced biases of 10.0% to 30.0% from lower to higher EVI thresholds [111]. These studies support the notion that, if threshold values are changed, biases can occur in land cover classification and assessment.

4. Conclusions

This study is the first to attempt to classify forest cover categories using the PBTC method. We used multitemporal Landsat 5 TM images for land cover mapping in 2001, 2006, and 2010, while Landsat 8 OLI collections were acquired in 2016 and 2020, with 30 m spatial resolution. This is the first study to have found EVI minimum and maximum threshold values for dry monsoon and open forests for PBTC classification in GEE to classify the land cover classes during the dry season. This season is preferable due to the low cloud cover, and an appropriate selection of EVI threshold values for each land cover class can reduce bias in the land cover classification. Notably, the average accuracy and kappa coefficient of land cover maps were 87% and 0.83, respectively. This shows that a combined method of PBTC and GEE can be used to assess the forest cover changes with higher accuracy, which is especially important in regions with limited historical data, mainly in developing countries. The annual forest cover change was −105.18 ha, and forest area decreased by 0.09% per annum between 2001 and 2020. From 2001 to 2010, forest cover increased by 4.62%, and district forest declined by 6.04% during 2010–2020. The annual rate of loss and gain in dry monsoon and open forest were 0.3% and 3.62%, respectively. Among four DS divisions, we observed the highest forest loss in the Vavuniya North division, at 19.2%. This can hasten negative consequences because of the continuous loss of dry monsoon forests. Changes in forest cover resulted in carbon emissions of 165,306.6 MgCO₂, and removals of 24,064.5 MgCO₂ annually, between 2001 and 2020.

To the best of our knowledge, this study is the first to develop a subnational FREL in Sri Lanka. Developing countries such as Sri Lanka, where data limitations occur, can utilize the satellite-based inventory through GEE to obtain long-term forest cover data. Therefore, the information on district-level forest cover, and carbon emissions and removals could provide the basis for the Sri Lankan government to determine the subnational FREL, against which, activity-based emission performance during the implementation period of the Paris Agreement between 2020 and 2030 can be compared, and the emission reductions could be claimed for financial support in future under the result-based payment of the REDD+ scheme of the UNFCCC. This study could provide support to environmentalists, forest conservationists and scientists, and decision and policymakers, to integrate the forest sector into future climate change plans to maximize the benefits. Appropriate strategies could provide a district-level contribution to achieving the NDC targets, while gaining benefits for national development.