Exploring Forest Change Spatial Patterns in Papua New Guinea: A Pilot Study in the Bumbu River Basin

Papua New Guinea is a country in Oceania that hosts unique rain forests and forest ecosystems which are crucial for sequestering atmospheric carbon, conserving biodiversity, supporting the livelihood of indigenous people, and underpinning the timber market of the country. As a result of urban sprawl, agricultural expansion, and illegal logging, there has been a tremendous increase in land-use land cover (LULC) change happening in the country in the past few decades and this has triggered massive deforestation and forest degradation. However, only a few studies have ventured into quantifying the long-term trends and their associated spatial patterns—and have often presented contrasting responses. Herein, we intended to assess the extent of deforestation and the rate of urbanization that happened in the past 33 years (1987–2020) in the Bumbu river basin in Papua New Guinea using satellite imagery—for the years 1987, 2002, 2010, and 2020—and Geographic Information System (GIS) tools. On performing image classification, land use maps were developed and later compared with Google Earth’s high-resolution satellite images for accuracy assessment purposes. For probing into the spatial aspects of the land-use change issues, the study area was divided into four urban zones and four forest zones according to the four main cardinal directions centered in the urban and forest area centers of the 1987 image; subsequently, the rate of urban area expansion in each urban zone was separately calculated. From our preliminary analysis and literature survey, we observed several hurdles regarding the classification of regenerative forests and mixed pixels and gaps in LULC studies that have happened in Papua New Guinea to date. Through this communication paper, we aim to disseminate our preliminary results, which highlight a rapid increase in urban extent from 14.39 km2 in 1987 to 23.06 km2 in 2020 accompanied by a considerable decrease in forest extent from 76.29 km2 in 1987 to 59.43 km2 in 2020; this observation favors the presumption that urban and agricultural land expansion is happening at the cost of forest cover. Moreover, strategies for addressing technical issues and for integrating land-use change with various socioeconomic and environmental variables are presented soliciting feedback.


Introduction
Land-use land cover (LULC) change is a very well-known and one of the most significant research topics in domains of climate change, urban planning, and forestry [1,2]. Herein, land-use denotes how people use a particular piece of land, and land cover is indicative of the physical land type such as forest or barren land. Having an in-depth understanding of land-use planning is vital as it opens up new opportunities to mitigate the detrimental effects of adopted land-use practices and thereby contributes towards sustainable use of resources at our disposal [3][4][5]. Nevertheless, the issue of deforestation, especially due to urbanization, has been a major concern globally due to its multi-spatial and multi-dimensional facets affecting copious amounts of ecosystem services from region-level species richness to landscape-level precipitation anomalies [6][7][8][9]. Moreover, the consequent decrease in the carbon sequestration capacity of forests and associated greenhouse gas emissions-which account for approximately 20% of all annual global emissions-stimulate prolonged distress worldwide given their role in mitigating global climate change [10,11].
Papua New Guinea is a country in Oceania that hosts unique rain forests and forest ecosystems which are crucial for sequestering atmospheric carbon, conserving biodiversity, supporting the livelihood of indigenous people, managing catchment, and underpinning the timber market of the country. However, as a result of urban sprawl, agricultural expansion, and illegal logging, there has been a tremendous increase in land-use change happening in the country in the past few decades and this has triggered massive deforestation and forest degradation [12,13]. Although the United Nations Framework Convention on Climate Change (UNFCCC) has attempted to implement strategies to reduce emissions from deforestation and degradation (REDD) in developing countries, only a few studies have ventured into quantifying the long-term land-use trends, their spatial patterns, and associated carbon emissions in Papua New Guinea. This paucity of validations can be primarily attributed to limited data availability and concurrent socioeconomic constraints. Also, the existing literature presents contrasting opinions, with one faction portraying alarming rates of deforestation occurring due to urbanization and another stating the former compilations to be an overestimation of reality [14,15]. In this regard, remote sensing techniques-especially satellite imagery-in association with Geographic Information System (GIS) analysis tools, can be used to obtain detailed, frequent, and up-to-date information of LULC patterns; the advancements in these technologies have revolutionized the LULC studies.
Since many tropical forests are in the frontier stage undergoing land-use transformations due to logging and agricultural expansion, several recent remote sensing studies-not only satellite-based endeavors, but also initiatives from drones and LiDAR (Light Detection and Ranging) spheres-have focused on these aspects in addition to analyzing the impacts of deforestation from the perspective of increased carbon emissions [16][17][18][19][20][21]. By identifying inherent trends, it was found possible to further investigate the rate of land-use change happening over the years and the motives behind these land-use transformations. Pratihast et al. [22] explored the integration of satellite data and community-based observations for forest monitoring and highlighted the importance of involving local communities in enhancing implementation of the REDD framework in regional areas. By contrast, Popkin [23] emphasized the large-scale significance of satellite imagery for enabling alert-systems to detect illegal logging in near-real time and sending warnings to environmental managers when a forested area is endangered. Among the satellite imagery used, the open-source Landsat imagery series, which provides data on changing landscapes from the year 1972, is one of the most commonly used spatial databases [24]. Banskota et al. [25] offer a review of recent advancements in forest monitoring achieved using Landsat Time Series Data and applaud the steady levels of development that have occurred in terms of science and application capacity. Studies assessing the consequence of land-use change on agricultural productivity, climate anomalies, and biodiversity are also commonplace in the satellite remote sensing realm [26][27][28][29].
In this study, we assessed the extent of deforestation and the rate of urbanization that has happened in the past 33 years (1987-2020) in the Bumbu river basin that runs through Lae, the capital of Morobe Province. Previous GIS analyses of land-use change in this area have provided brief snapshots of a growing metropolitan area and of the impact of this growth on the surrounding forested region [12,30], yet there are no extensive studies to approximate the extent of forest degradation and spatial variability happening in this area. Herein, for analysis purposes, we chose Landsat data, since the repository includes images that are consistent with data from previous missions, thus allowing for long-term change detection at various scales [31].
Through this communication paper, we aim to disseminate our preliminary results, which highlight a rapid increase in urban extent, accompanied by a loss of forest cover, in the past 33 years. In addition, we look forward to receiving feedback on our presented interpretations regarding the land-use land cover change patterns, proposed strategies for improving adopted classification approaches, and plans for integrating land-use change with various socioeconomic and environmental variables.

Study Area
The study focuses on the Bumbu river basin in Papua New Guinea (see Figure 1), located between the Markham river basin to the west, and the Busu river basin to the east, cutting through the city of Lae, the capital of Morobe Province in the South. Lae is an industrial city, and the second largest in Papua New Guinea, located at Lat/Long: −6.66942, 146.97911. The Bumbu river is considered narrow with a flow of medium speed. However, the river flows and erodes at a much higher rate during the flood season [32], as it crosses a basin consisting mainly of sandy loam [30]. Such floods can be detrimental to the floodplain of the affected stream, which for the Bumbu main channel includes the industrial city of Lae. Uplift along the Markham valley within Lae from tectonic plate convergence [33] provides the area with a unique confluence of three major rivers with interspersed plains across which the Lae urban district is now spread.
Land 2020, 9, x FOR PEER REVIEW  3 of 18 forested region [12,30], yet there are no extensive studies to approximate the extent of forest degradation and spatial variability happening in this area. Herein, for analysis purposes, we chose Landsat data, since the repository includes images that are consistent with data from previous missions, thus allowing for long-term change detection at various scales [31]. Through this communication paper, we aim to disseminate our preliminary results, which highlight a rapid increase in urban extent, accompanied by a loss of forest cover, in the past 33 years. In addition, we look forward to receiving feedback on our presented interpretations regarding the land-use land cover change patterns, proposed strategies for improving adopted classification approaches, and plans for integrating land-use change with various socioeconomic and environmental variables.

Study Area
The study focuses on the Bumbu river basin in Papua New Guinea (see Figure 1), located between the Markham river basin to the west, and the Busu river basin to the east, cutting through the city of Lae, the capital of Morobe Province in the South. Lae is an industrial city, and the second largest in Papua New Guinea, located at Lat/Long: −6.66942, 146.97911. The Bumbu river is considered narrow with a flow of medium speed. However, the river flows and erodes at a much higher rate during the flood season [32], as it crosses a basin consisting mainly of sandy loam [30]. Such floods can be detrimental to the floodplain of the affected stream, which for the Bumbu main channel includes the industrial city of Lae. Uplift along the Markham valley within Lae from tectonic plate convergence [33] provides the area with a unique confluence of three major rivers with interspersed plains across which the Lae urban district is now spread.

Data Collection/Preprocessing
The primary data used in this study were derived from a series of Landsat satellite images, path

Data Collection/Preprocessing
The primary data used in this study were derived from a series of Landsat satellite images, path forested areas. Google Earth images were used for accuracy assessment and to aid in choosing the training samples; watershed mapping was used to delineate the study area boundary. The Bumbu river main channel was digitized based on ESRI (Environmental Systems Research Institute) imagery for the year 2020 and used at a later stage for illustration. The images were subjected to atmospheric correction by conversion to Top of Atmosphere (ToA) reflectance, and solar correction was accounted for according to Young et al. [34], and USGS [35,36]. Images were subsequently clipped to the watershed boundary. Geometric correction and calibration were deemed unnecessary for the purposes of this study as the Landsat images were geometrically corrected by NASA (National Aeronautics and Space Administration) [34]. The overall workflow is presented below (Figure 2). Land 2020, 9, x FOR PEER REVIEW 4 of 18 forested areas. Google Earth images were used for accuracy assessment and to aid in choosing the training samples; watershed mapping was used to delineate the study area boundary. The Bumbu river main channel was digitized based on ESRI (Environmental Systems Research Institute) imagery for the year 2020 and used at a later stage for illustration. The images were subjected to atmospheric correction by conversion to Top of Atmosphere (ToA) reflectance, and solar correction was accounted for according to Young et al. [34], and USGS [35,36]. Images were subsequently clipped to the watershed boundary. Geometric correction and calibration were deemed unnecessary for the purposes of this study as the Landsat images were geometrically corrected by NASA (National Aeronautics and Space Administration) [34]. The overall workflow is presented below (Figure 2).

Pre-Classification
The images were duplicated to visualize the true colors and the false colors using three spectral indices, i.e., Normalized Difference Vegetation Index (NDVI) (which helps in discriminating vegetation), Normalized Difference Built-up Index (NDBI) (to assist in identifying urban area presence in imagery), and Normalized Difference Water Index (NDWI) (to enhance open water appearance in remotely sensed imagery) [37]. These indices were derived for each image to simplify mapping of forest, urban area, and open water classes, respectively [37]. The images were later used to aid manual interpretation during the reclassification process.

Supervised Classification
A supervised classification using the Maximum Likelihood Classification (MLC) algorithm was performed for image classification purposes [38,39]. The MLC method is known for its strong theoretical foundation, and for its ability to accommodate varying data. This method was applied in our study using ArcGIS Pro 2.5 [40]. The four main classes of the classification scheme and their descriptions are provided in Table 1. The forest class included two subclasses, namely, dense forest and forest, which were merged after running the MLC. Thirty training sample polygons were created and evaluated in each class/subclass for each image. The training sites were selected using visual interpretation of Google Earth's high-resolution satellite images and the NDVI, NDBI, and NDWI images to compare with Landsat images.

Pre-Classification
The images were duplicated to visualize the true colors and the false colors using three spectral indices, i.e., Normalized Difference Vegetation Index (NDVI) (which helps in discriminating vegetation), Normalized Difference Built-up Index (NDBI) (to assist in identifying urban area presence in imagery), and Normalized Difference Water Index (NDWI) (to enhance open water appearance in remotely sensed imagery) [37]. These indices were derived for each image to simplify mapping of forest, urban area, and open water classes, respectively [37]. The images were later used to aid manual interpretation during the reclassification process.

Supervised Classification
A supervised classification using the Maximum Likelihood Classification (MLC) algorithm was performed for image classification purposes [38,39]. The MLC method is known for its strong theoretical foundation, and for its ability to accommodate varying data. This method was applied in our study using ArcGIS Pro 2.5 [40]. The four main classes of the classification scheme and their descriptions are provided in Table 1. The forest class included two subclasses, namely, dense forest and forest, which were merged after running the MLC. Thirty training sample polygons were created and evaluated in each class/subclass for each image. The training sites were selected using visual interpretation of Google Earth's high-resolution satellite images and the NDVI, NDBI, and NDWI images to compare with Landsat images.

Re-Classification
Manual interpretation was used to enhance the accuracy of the resulting classified images, using the re-classifier tool provided in ArcGIS Pro, with the assistance from Google Earth high-resolution images and the aforementioned three spectral indices. In addition, 100 random sampling points (50 forest and 50 non-forest) were evaluated using Google Earth and subsequently used to enhance the reclassification process of the 1987 and 2020 images, and to assist in the assessment of the extent of forest gain versus forest loss.

Accuracy Assessment
To characterize image classification accuracy, a confusion matrix, which is a table that depicts the correlation of ground truth points and the result of the classification, was used to calculate and compare the overall percent accuracies, the user's accuracies, and producer's accuracies [38][39][40][41]. Herein, the user's accuracy portrays the accuracy levels from the point of view of a map user and is a complement of the commission error; producer's accuracy reflects the map accuracy from the point of view of a map maker and is a complement of the omission error; and the kappa coefficient gives a measure of the overall accuracy of a classification. For each classified image, 200 random points were generated using ArcGIS Pro (50 points in each class), and the pixel values were extracted from the classified images to the points. Google Earth's high-resolution satellite images and Google Earth Pro software were used to examine the ground truth status of reference points and to assign the ground status to the matrix. After developing the matrix, the kappa coefficient (which is a statistical multivariate technique used widely in accuracy assessment) was calculated for each classified image; respective omission errors and commission errors were also determined.
Additionally, a GIS-based analysis using ArcGIS Pro 2.5 was conducted on the classified images to better describe and quantify the change of the built-up area and forest area over the study period. The LULC change maps were overlapped, and the center of urban area in 1987 was located and used to divide the urban area into four zones; Urban NE (northeastern), Urban NW (northwestern), Urban SW (southwestern), and Urban SE (southeastern). Following the same approach, the center of the forest area in 1987 was calculated using GIS and was used to conduct in-depth change analysis of the four forest zones (NE, NW, SW, SE).

Results
The overall classification accuracies for the 1987, 2002, 2010, and 2020 maps were 90.5%, 87%, 93%, and 90.5%, respectively, while the associated kappa coefficients were 0.87, 0.83, 0.91, and 0.87. Table 2 shows the accuracy assessment results along with the error matrix and the producer's and user's accuracy for each classified image. The resulting classified images showing the pattern of land-use/land cover in the study area over time are represented in Figure 3, while the spatial distribution of each class can be observed in Table 3. The results show a progressive decrease in the forest area over the studied period as it dropped from 62.33% in 1987 to 57.26% in 2002 to 54.72% in 2010 and then to 48.55% in 2020. On the other hand, rapid increases in the "Built-up" class and the "Other" class occurred over the same period, indicating the expansion of urban areas and other land type classes at the expense Land 2020, 9, 282 6 of 18 of a reduction in forested areas. The built-up area increased from 11.76% in 1987 to 18.84% in 2020, passing 13.65% in 2002 and 15.71% in 2010, while the other types of land use including agricultural land expanded significantly over the study period, from 24.81% of the study area in 1987 to 31.86% in 2020. This change coincides with overall trends of land-use change in Morobe Province between 1975 and 2000 as found by Ningal et al. [12], and as confirmed by Samanta and Paul [42], who studied LULC change in Lae district between 1992 and 2013.
Land 2020, 9, x FOR PEER REVIEW 6 of 18 including agricultural land expanded significantly over the study period, from 24.81% of the study area in 1987 to 31.86% in 2020. This change coincides with overall trends of land-use change in Morobe Province between 1975 and 2000 as found by Ningal et al. [12], and as confirmed by Samanta and Paul [42], who studied LULC change in Lae district between 1992 and 2013.    Table 4. The southeastern zone (SE) and the northeastern zone (NE) experienced the largest share of the expansion in built-up areas over the period 1987 to 2020, while the expansion in the northwestern zone (NW) was significant in the overall studied period, and the second sharpest increase, after the southeastern zone (SE), in urban classifications occurred in the last ten years (as shown in Figures 4-6). These significant sprawls in the northern zones of urban areas, where the forest cover is dense, along with the overall observed expansion of other land type class encouraged us to substantiate the contention that the forest losses are most likely due to urban and agriculture expansion. Moreover, the overall aerial extent of urbanization increased rapidly and peaked in the period from 2010 to 2020. The lowest rate of expansion was experienced by the southwestern (SW) zone, which was already densely urban in its southern part prior to 1987. The analysis indicates expanding the urban area to the north puts the lands in these northern zones under stress due to this urban expansion ( Figure 5).  largest share of the expansion in built-up areas over the period 1987 to 2020, while the expansion in the northwestern zone (NW) was significant in the overall studied period, and the second sharpest increase, after the southeastern zone (SE), in urban classifications occurred in the last ten years (as shown in Figures 4-6). These significant sprawls in the northern zones of urban areas, where the forest cover is dense, along with the overall observed expansion of other land type class encouraged us to substantiate the contention that the forest losses are most likely due to urban and agriculture expansion. Moreover, the overall aerial extent of urbanization increased rapidly and peaked in the period from 2010 to 2020. The lowest rate of expansion was experienced by the southwestern (SW) zone, which was already densely urban in its southern part prior to 1987. The analysis indicates expanding the urban area to the north puts the lands in these northern zones under stress due to this urban expansion ( Figure 5).

Urban areas expansion between 1987 and 2020
Urban NE zone Urban NW zone Urban SW zone Urban SE zone

Urban areas expansion between 1987 and 2020
Urban NE zone Urban NW zone Urban SW zone Urban SE zone The in-depth analysis of forest area changes shows an overall rapid loss in forest extent in all zones over all the time intervals. The forest extent in the southeastern forest zone decreased significantly, especially in the last ten years. This decrease could be attributed to the expansion of built-up areas and to the increase in the extent of other types of land cover. Forest loss in the northeastern and northwestern forested zones is relatively lower in comparison to the southern zones; however, both northern zones have experienced forest loss. Figure 7 shows the spatial pattern of forest loss in all the studied time intervals. Zonal and total changes in the extent of forested areas along with the yearly rate of change for three time intervals within the study period 1987-2020 are presented in Table 5.

(km 2 )
year −1 ) (km 2 ) year −1 ) (km 2 ) year −1 ) (km 2 ) year −1 ) Urban NE The in-depth analysis of forest area changes shows an overall rapid loss in forest extent in all zones over all the time intervals. The forest extent in the southeastern forest zone decreased significantly, especially in the last ten years. This decrease could be attributed to the expansion of built-up areas and to the increase in the extent of other types of land cover. Forest loss in the northeastern and northwestern forested zones is relatively lower in comparison to the southern zones; however, both northern zones have experienced forest loss. Figure 7 shows the spatial pattern of forest loss in all the studied time intervals. Zonal and total changes in the extent of forested areas along with the yearly rate of change for three time intervals within the study period 1987-2020 are presented in Table 5.   The rates of forest loss in the period between 2010 and 2020 in all zones were significantly higher (−0.517, −0.035, −0.072, −0.131 km 2 year −1 for the zones SW, SW, NW, NE, respectively) in comparison with other periods, with the highest rate of loss occurring in the southeastern forest zone. It should be noted that the northwestern zone experienced a modest gain in forest extent from 1987-2002. On the other hand, there was no gain in the forest extent in any other zone in any other time period; nevertheless, it should be borne in mind that the absence of forest gain in the total area may not account for the marginal forest gain that appear in the maps and charts (Figures 7-10) as the lost area was larger than the gained one; the total change in Table 5 appears only as a loss. Figure 9 shows the gain-loss rate in all zones and it can be noticed that the southeastern zone has witnessed the highest rate of loss and the largest loss of forested area during the study period. The rates of forest loss in the period between 2010 and 2020 in all zones were significantly higher (−0.517, −0.035, −0.072, −0.131 km 2 year −1 for the zones SW, SW, NW, NE, respectively) in comparison with other periods, with the highest rate of loss occurring in the southeastern forest zone. It should be noted that the northwestern zone experienced a modest gain in forest extent from 1987-2002. On the other hand, there was no gain in the forest extent in any other zone in any other time period; nevertheless, it should be borne in mind that the absence of forest gain in the total area may not account for the marginal forest gain that appear in the maps and charts (Figures 7-10) as the lost area was larger than the gained one; the total change in Table 5 appears only as a loss. Figure 9 shows the gain-loss rate in all zones and it can be noticed that the southeastern zone has witnessed the highest rate of loss and the largest loss of forested area during the study period.      The rates of forest loss in the period between 2010 and 2020 in all zones were significantly higher (−0.517, −0.035, −0.072, −0.131 km 2 year −1 for the zones SW, SW, NW, NE, respectively) in comparison with other periods, with the highest rate of loss occurring in the southeastern forest zone. It should be noted that the northwestern zone experienced a modest gain in forest extent from 1987-2002. On the other hand, there was no gain in the forest extent in any other zone in any other time period; nevertheless, it should be borne in mind that the absence of forest gain in the total area may not account for the marginal forest gain that appear in the maps and charts (Figures 7-10) as the lost area was larger than the gained one; the total change in Table 5 appears only as a loss. Figure 9 shows the gain-loss rate in all zones and it can be noticed that the southeastern zone has witnessed the highest rate of loss and the largest loss of forested area during the study period.     A graphical portrayal of the spatial distribution of gains and losses of forested areas across the watershed from 1987 to 2020 is presented in Figure 10. This map was derived from the classified images of 1987 and 2020 by applying pixel-to-pixel comparison with subsequent reclassification of the result into four different classes showing (a) forest loss, (b) forest gain, (c) consistent forest extent, and (d) unchanged other land types ( Figure 10). Moreover, the areas were calculated and the results are presented in Table 6. . Figure 10. Overall forest gain-loss between 1987 and 2020.
A graphical portrayal of the spatial distribution of gains and losses of forested areas across the watershed from 1987 to 2020 is presented in Figure 10. This map was derived from the classified images of 1987 and 2020 by applying pixel-to-pixel comparison with subsequent reclassification of the result into four different classes showing (a) forest loss, (b) forest gain, (c) consistent forest extent, and (d) unchanged other land types ( Figure 10). Moreover, the areas were calculated and the results are presented in Table 6.

Discussion
On conducting LULC analysis, we observed a progressive decrease in the forest area in the Bumbu river basin for all the time intervals considered, with the maximum change happening between 2010 and 2020. These findings are crucial as there have been very few studies that track the change in deforestation patterns in Papua New Guinea in the past decade, even though around 63.4% of the country consisted of forested areas as of 2010; and endeavors exploring spatial patterns of forest loss within particular study areas are even more scarce. In our case, the study area had 62.33% forest cover in 1987 which dwindled down to 57.26% in 2002 and then to 54.72% in 2010 and afterwards to 48.55% in 2020. Therefore, there exists a need to evaluate if a similar pattern occurs at a national scale.

Discussion
On conducting LULC analysis, we observed a progressive decrease in the forest area in the Bumbu river basin for all the time intervals considered, with the maximum change happening between 2010 and 2020. These findings are crucial as there have been very few studies that track the change in deforestation patterns in Papua New Guinea in the past decade, even though around 63.4% of the country consisted of forested areas as of 2010; and endeavors exploring spatial patterns of forest loss within particular study areas are even more scarce. In our case, the study area had 62.33% forest cover in 1987 which dwindled down to 57.26% in 2002 and then to 54.72% in 2010 and afterwards to 48.55% in 2020. Therefore, there exists a need to evaluate if a similar pattern occurs at a national scale. Our results are also in tandem with the trend observed in a previous study by Ningal et al. [12], who reported a decline in forest area in Morobe Province, from 9.8 ha/person in 1975 to 4.4 ha/person in 2000. However, two major differences we have in our study are i) our study states urbanization as one of the core reasons for the reduced forest cover in addition to agricultural land expansion, whereas in the former case, subsistence agriculture was listed as the primary cause and ii) instead of solely depending on satellite imagery for performing LULC analysis, possibly due to the lack of available data, in the former case, the authors relied on topographic maps for the year 1975. Notwithstanding the minor distinctions, the primary underlying reasons in both cases were attributed to the rise in the population. Another study by Samanta and Pal [42] derived a change detection matrix of 1992-2013/14 from Landsat-based observations and reported the occurrence of extensive urban sprawl in Lae, the capital of Morobe Province, which was followed by a decline in natural forests by around 6977.25 hectares within a span of 20 years. Given that our study area is a watershed region, concerns with respect to a change in the course of Busu river in Lae city and concomitant destruction of extant land cover patterns, as expressed Land 2020, 9, 282 13 of 18 in Sekac and Jana [32], should also be taken into consideration, and we will be weighing this factor while executing the task mentioned in Section 6.3.
A few of the major issues we came across during our pilot study were related to the cloud cover. For the years of interest it was not possible to acquire cloud-free satellite imagery for similar time periods of the year. In the current study, for the years 1987, 2002, 2010, and 2020, the months of October, November, February, and May were considered, respectively. This mismatch in time periods might have placed constraints on inter-annual comparisons as the vegetation zones with a net phenological cycle can be distorted in this case. Consequently, the "Other" class, which includes cultivated lands, could have also been prone to underestimations due to the aforementioned reason; for instance, the cases where post-harvested crop areas might be tagged as bare soil; similarly this might result in some built areas having no real link with urbanization. Therefore, for our upcoming comprehensive study, we intend to adopt a standard LULC classification system with the following classes: forest land, barren, wetland, water, agriculture, rangeland, and built-up areas. In this way, we expect to be capable of teasing out small vegetation areas, isolated small buildings, and similar areas and to be a better judge of LULC activities happening in the study area. Also, for combating loopholes arising from the low resolution of Landsat sensors, Sentinel-1 and Sentinel-2 data-or a fusion of Sentinel and Landsat data-are considered as an alternative. The 10 m resolution open-source images offered at a high revisit time by Sentinel sensors would be able to help in selecting cloud-free images and to identify the best periods to observe different vegetation covers in the years post 2014 [43][44][45]. Additionally, employing semi-automatic approaches in land transformation analyses-which would allow the construction of models based on predefined classification-can be deemed beneficial for performing land-use change monitoring tasks more effectively; these paradigms can be extrapolated to include other formats of data as well, as listed in Wang et al. [46], where the authors presented a framework to analyze path-dependent regional industrial land transition processes using multi-year vector data and underscored opportunities for reducing data collection efforts while monitoring land-use transitions and their impacts on the surroundings. Few other drawbacks regarding this pilot study and potential workarounds are mentioned in Section 6.
For tracking deforestation and urban sprawl in real-time, venturing into the field of drone analytics has proved to be beneficial [17,[47][48][49][50]. With the advancement in image processing techniques and availability of state-of-the-art machine learning algorithms, it is now possible to perform drone-based forest monitoring operations more frequently and at semi-automated levels. Similarly, LiDAR (Light Detection and Ranging) technology also offers multiple benefits along this line; being an active sensor, it is capable of extracting high-definition information of vertical forest structures-unlike satelliteand drone-based remote sensing techniques-and hence plays a crucial role in studies focused on estimating biomass levels and forest structures in dense canopy areas. Nevertheless, the primary constraint while employing LiDAR is the associated cost, which limits its applicability for larger areas, and for experiments needing recurring data [51][52][53][54]. Additionally, combinations of drones and LiDAR are also gaining momentum [55,56]. However, for large-scale studies, satellite imagery can be considered to be the optimal choice.
Based on our results, we hope to attract more international attention, public participation, and local government and NGO involvement for addressing exigent requirements of a sustainable community. This will allow Papua New Guinea to track its domestic emissions of greenhouse gases in a more effective way for supporting REDD+ initiatives and to keep a check on the growing demand for biofuels, which poses a serious threat to existing forested landscapes of the country [57]. Currently, due to a lack of human resources and financial aid, there is a lag in the implementation of state-of-the-art forest monitoring techniques, and most of the existing forest distribution maps are out of date [58]. Therefore, by highlighting the urgency with respect to curbing deforestation and mitigating climate change, we champion the need for revamping the existing forest management paradigms and associated policies. Future directions of this project are described in Section 6.

Conclusions
On examining forest loss and gain, and urban area changes, in the Bumbu river watershed between 1987 and 2020 using GIS and remote sensing, our preliminary results indicated a progressive decrease in the total forested area during the study period, with the rate of decrease peaking in the last ten years between 2010 and 2020. The loss of forested land was most intense in the area to the northwest of Lae. Simultaneously, a rapid increase in the urban extent and other land types was observed across the study area at the expense of a decrease in forested areas. However, given the coarse resolution of Landsat imagery, we had issues with consistent representation of land areas. For subsequent studies, we intend to collaborate with Pacific-Australia Climate Change Science and Adaptation Planning (PACCSAP) and acquire LiDAR data for our study area as these approaches will help us gather more information about the tree populations in various stages of development and will be a valuable complement to the optical data from satellites for LULC studies. Additionally, applications of data from higher resolution imagery acquired using Unmanned Aerial Vehicles (UAVs) in a timely manner is perceived as an alternative from an economic point of view. Regardless of the data acquisition method, we can unanimously agree that up-to-date and site-specific, yet scalable approaches are a must for developing and implementing good management practice for assessing LULC given the transient changes happening in response to global climate change.

Annual Land-Use Land Cover (LULC) Change Analysis
In future analyses, as a part of this task, we will be utilizing annual time series satellite data-as opposed to imagery of the four time periods used in this preliminary study-to analyze land-use land change patterns. Moreover, two images with slightly different time periods will be considered for each year to minimize cloud disturbance. Adding to that, we look forward to evaluating how fusion of Sentinel data with Landsat imagery helps increase our chances of obtaining cloud-free imagery for the years post 2014 to overcome the presence of phenological uncertainties. From the preliminary analysis, we noticed a considerable decrease in forest extent at fluctuating rates among zones and over the studied period. Nevertheless, several changes and a sequence of multiple consecutive events can happen very rapidly and might be untraceable in 5+ year intervals; for instance, excessive tree mortality due to sudden drought or hurricane followed by the transition of that land into industrial structures. We also noticed that there were modest increases in forested landscapes in some areas, though it is still unclear whether this was due to misclassification of pixels or a real phenomenon. Nonetheless, with annual Landsat data we might be able to detect the exact changes happening year by year and have more options to overcome cloudy imagery issues; this could also reduce the error rate through better interpretation capability. Also, the study area itself gives significance to this work as the previous studies were done at a much larger scale. Identifying the urban sprawl directions and forest loss directions at this scale would give essential information to decision-makers and hopefully support implementation of preventative measures in regards to forest loss. For accurate quantification of forest regeneration, the annual information would be critical and it would be able to provide unique perspectives; the literature review reveals very few studies focusing on this subject, and moreover, some have expected continuous loss of forest without accounting for regenerating forests [12].
We also noticed increasing rates of urban expansion to be acute in the past ten years, especially in the northeastern and southeastern zones of the urban area. Annual data will allow us to evaluate if the urban expansion was spread out evenly or if there were critical extreme events that led to the varying rates at which expansion took place. Also, we plan to further investigate anomalies such as the negative rates of change of urban classifications in the urban northwestern zone between 1987 and 2002 which indicates a marginal disappearance of built-up areas; we expect this to be a result of misclassified pixels as reported in the error matrix, rather than a result of demolished urban areas. This will also help us explore issues with classifying transitions in forest lands within the Bumbu watershed-for instance, inference on how much of the dense forests transitioned into regenerative forest and vice versa-and shed more light on where Papua New Guinea should concentrate efforts to preserve forested land.

Improved Classification of Regenerated Forests
Given the relatively moderate resolution of Landsat satellite images, and given the small scale of the study area, one of the major challenges faced in this study was how to distinguish regenerating forested land from sparsely vegetated urban land. Despite the use of NDVI images, the attempts did not prove useful to delineate these differences. To overcome this issue and to likewise not mix the clustered trees in urban areas with forest, a generalization approach was developed during the manual reclassification process. Herein, the focus was to delineate the extent of the forest area and the extent of urban area rather than the footprint of the buildings. Given this focus, we have currently limited classification to the four classes chosen. However, in follow-up studies, we look forward to re-classifying imagery to include additional zones and to providing more detailed mixed pixel classification information by assigning thresholds or extra masks-as well as by adopting the use of improved spectral indices. In this way, it may be possible to evaluate whether harvested areas are lost forever or are in a state of regrowth that will eventually regenerate into mature forest; in the current study, such nuances are subsumed by the "other" category. For comparison purposes, neural networks will also be applied to analyze LULC change trends within the Bumbu river basin. As mentioned in the previous section, understanding whether new forests are generating or not is crucial as this helps answer many questions on forest degradation and assists in developing forest management policies and improved classification techniques.

Integration with Water Quality Data
In future tasks, we plan to integrate the results from the aforementioned annual LULC analysis into a parallel study focusing on the development of an analytic protocol for estimating the relative importance of socio-economic and environmental variables on water quality and their spatial and temporal variation within the Bumbu watershed. The protocol utilizes a 1 arc-second Digital Elevation Model (DEM) to extract the catchment area and the stream drainage pattern within the Bumbu watershed, thereby facilitating estimation of the runoff accumulation of precipitation, natural geologic and geochemical inputs, and human-generated detritus and pollutants in Bumbu waters. In this regard, having LULC data will be valuable for identifying emerging pollution zones-which can be attributed to the newly urbanized areas-and through ANOVA (analysis of variance) between LULC and water quality data, we will be assessing the possible impacts of deforested and built-up areas on deteriorating water quality. Funding: There was no external funding.