The Sigi River watershed is located in northeastern Tanzania, with geographic coordinates spanning from −5°12′S, 38°36′E to −4°48′S, 38°70′E. It is situated within the East Usambara Mountains and their foothills, transitioning into the coastal plain (Figure 1). Our study covers the watershed area defined by the outlet of the Mabayani Reservoir, which serves as a vital water resource for the greater Tanga region, home to Tanzania’s second-largest port. The total watershed area covers 887 km², with elevation ranging from 90 m to 1,200 m asl. The region has experienced substantial environmental changes over decades, particularly due to deforestation driven by the expansion of smallholder agriculture linked to population growth. According to Yanda and Munishi (2007), forest cover in the Sigi River watershed was reduced to 33% (3,896 ha) between 1955 and 1995, while cultivated land increased by 37% (4,341 ha) over the 11,728 ha area of the watershed. This deforestation has led to severe soil erosion and sedimentation, which have adversely impacted downstream water bodies, including the Mabayani Reservoir. The siltation has threatened the reservoir’s operational capacity, reducing its storage potential and posing a risk to urban water supply (NERC, 1994; Yanda and Munishi, 2007). Population figures for the study area have increased from 169,200 in the 1980s to 238,260 in 2022 (URT, 2022).

The climate of the region is humid, characterized by a bimodal rainfall pattern. Rainfall occurs primarily from March to May and from October to December (Hamilton and Bensted-Smith, 1989), with a dry season between January and February, marked by occasional rainfall events. The upland areas receive an average annual rainfall of 1,540 mm and have temperature around 20°C, as observed at the Amani station, while the lowland Lanconi station near the reservoir records an average of 1,173 mm of rainfall and temperature averaging 27°C (Figure 2).

Geologically, the catchment is dominated by crystalline bedrocks in the mountain block and escarpment regions, with sedimentary bedrocks found in the foothills and coastal plain. The soil in the mountainous areas are primarily classified as Acrisols, which are nutrient-poor, acidic, and loamy in texture resulting in low soil fertility that is quickly depleted by continuous cultivation (Hartemink, 1997; Kirsten et al., 2016). In contrast, soils at lower elevations are more fertile and better suited for agriculture (Hartemink, 1997). The alluvial valleys are characterized by Fluvisols, which are nutrient-rich and more conducive to cultivation (Ministry of WaterUNDP and GEF, 2019).

This study used Landsat imagery from Landsat 5 (TM) for 1983, 1996, and 2009, and Landsat eight for 2022, obtained from the USGS via Earth Explorer (http://earthexplorer.usgs.gov). The selected time intervals were chosen to capture long-term LULC trends, minimizing short-term variability and seasonal fluctuations. Images were acquired during December to February, coinciding with the short dry season, which is typically characterized by reduced rainfall and lower cloud cover, increasing the likelihood of cloud-free imagery suitable for analysis. For visual interpretation, field verification, and the assignment of LULC labels for the selected historical years (1983, 1996, and 2009), high-resolution imagery from Google Earth and SPOT 4 (10 m resolution) were also used. These images were crucial for accurate LULC classification during the training data collection. Additionally, a 30-m-resolution Digital Elevation Model (DEM) from the Shuttle Radar Topography Mission (SRTM) was downloaded from the USGS and processed using ArcGIS Pro (version 3.1.0) to derive topographic attributes. Before extracting vegetation indices (VIs), the Landsat Level 2 Science Product (L2SP) data underwent several preprocessing steps to ensure high-quality inputs for analysis. Atmospheric correction was applied to convert the data to Top of Atmosphere (TOA) reflectance. Radiometric correction, including adjustments for TOA reflectance and brightness temperature, was performed to minimize radiometric distortions. Geometric correction was applied to ensure precise spatial alignment of the imagery. Additionally, topographic correction was carried out using the C-correction algorithm to account for the effects of terrain-induced illumination variations (Yin et al., 2022). These preprocessing tasks were executed using ArcGIS Pro and QGIS software to prepare the data for vegetation index extraction and subsequent land use/land cover analysis.

The LULC classes was identified using a combination of field surveys and expert knowledge, which were instrumental in defining 13 distinct LULC classes (Table 1). The study area was stratified into four elevation zones: Upland (>850 m), Escarpment (350–850 m), Foothill (175–350 m), and Lowland (<175 m). Additionally, the slope was classified into four categories based on terrain steepness: gentle (<5°), moderate (5°–15°), steep (15°–35°), and very steep (>35°) (Figure 3). These classification criteria align with approaches used in previous studies (e.g., Wondie et al., 2012; Birhane et al., 2019; Wang and Cheng, 2022), which underscore the influence of both elevation and slope on land cover patterns.

To ensure representative sampling for LULC classification, a stratified random sampling method was applied within each elevation zone following the approach of Olofsson et al. (2014). Specifically, 100 random samples were selected for each LULC class within each elevation zone, resulting in a total of 5,200 samples per year. Historical training data for the years 1983, 1996, and 2009 were obtained from high-resolution imagery (Google Earth, SPOT 4), and these samples were cross validated with field data to ensure the correctness of the LULC labels. To address potential bias from equal allocation of samples across classes, we additionally applied weighted stratified random sampling, where the number of validation samples was proportional to the actual area of each class. The accuracy metrics obtained through this approach showed only minimal differences compared to the original method, confirming the robustness of the classification results. Precision metrics were also included to provide a more comprehensive evaluation of model performance. To select the most appropriate classification model, we evaluated multiple machine learning algorithms, including Random Forest (RF) and Support Vector Machines (SVM), during the model development phase. Random Forest consistently produced the highest overall accuracy and class-specific precision, especially in complex, topographically variable areas. Given its robustness and reliability, RF was selected for final classification.

LULC classification was performed using the Random Forest algorithm in ArcGIS Pro. The model incorporated Landsat spectral bands (B1–B6), along with several derived indices including the Normalized Difference Vegetation Index (NDVI), Normalized Difference Built-up Index (NDBI), Modified Soil-Adjusted Vegetation Index (MSAVI), Soil-Adjusted Vegetation Index (SAVI), and Modified Normalized Difference Water Index (MNDWI). A supervised classification approach was also employed, using expert-driven sampling points and SPOT four imagery to refine results. The Random Forest classifier, was configured with 200 trees and a maximum depth of 500, ensuring strong predictive performance. The dataset was divided into 70% for training and 30% for validation, enabling a thorough assessment of classification accuracy.

Post-classification corrections were carried out using the ArcGIS Pro Pixel Editor to rectify unrealistic transitions between LULC classes, such as implausible changes (e.g., from agricultural cultivation to forest) that are inconsistent with natural land cover dynamics (Figure 4).

In addition to the LULC classification, temporal land cover changes were analyzed using the OpenLand package in R. Sankey diagrams were created to illustrate transitions between LULC classes across the years 1983, 1996, 2009, and 2022, following the methodology proposed by Aldwaik and Pontius (2012). Such diagrams provide a comprehensive visual representation of land cover dynamics, highlighting significant changes over time. To further explore the spatial distribution of these changes, the final LULC maps were overlayed with elevation and slope layers. This spatial analysis offered valuable insights into the relationship between LULC changes and topographic factors, enabling the identification of spatial patterns and trends in the study area.

The performance of the classification model was evaluated for the years 1983, 1996, 2009, and 2022 using four key metrics: Accuracy (ACC), recall, precision and the F1-Score, based on a 30% sample of the validation data. Accuracy was calculated as the ratio of correctly predicted pixels (both true positives and true negatives) to the total number of pixels, as shown in Equation 1: ACC = TP + TN TP + TN + FP + FN (1) where:

TP (True Positive): Number of features correctly predicted for the observed category; TN (True Negative): Number of features correctly predicted as not belonging to the observed category; FP (False Positive): Number of features incorrectly predicted as belonging to the observed category; FN (False Negative): Number of features incorrectly predicted as not belonging to the observed category.

Recall was calculated as the proportion of correctly predicted features within a given category, as expressed in Equation 2: R e c a l l = T P T P + F N (2) where:

TP: True Positives (correctly predicted instances of the observed category); FN: False Negatives (incorrectly predicted instances that do not belong to the observed category).

Precision was calculated as the proportion of correctly predicted positive instances out of all instances predicted as positive. It indicates how many of the predicted positive cases are actually correct. As expressed in Equation 3: P r e c i s i o n = T P T P + F P (3)

The F1-score was computed as the harmonic mean of precision and recall, expressed in Equation 4: F 1 − Score = 2 * T P 2 * T P + FP + F N (4)

Where, TP: True Positives (correctly predicted instances of the observed category); FN: False Negatives (incorrectly predicted instances not belonging to the observed category); FP: False positive (incorrectly predicted instances belonging to the observed category).

The LULC distribution in the Sigi River catchment for 2022 (Figure 5; Table 2) reveals that small-scale cultivation is the dominant land use class, covering 360.4 km², primarily in the foothills and lowland regions, with minimal coverage in the upland. Shrubland (213.4 km²) is the second most prevalent class, concentrated in the foothills, followed by the escarpment and least in the upland areas.

Forest (primary and secondary), covering 154.3 km², is predominantly located in the upland regions, followed by the escarpment and with the least area in the lowland. Tea plantations (34 km²) and spice agroforestry (21.1 km²) are concentrated in the escarpment and upland, while sisal plantations (31.5 km²) are most common in the lowland and foothills, and teak plantations (27.9 km²) and citrus agroforestry (14.0 km²) are predominantly found in the foothills. The builtup class, covering 12 km², is primarily located in the foothills and lowlands. These patterns reflect the topographic and climatic suitability of different regions in the watershed, with agricultural activities and plantations being most concentrated in the lower elevation zones (foothills and lowland), while the upland and escarpment areas are more dominated by forests (primary and secondary) and specialized agroforestry systems.

The performance of the classification model was evaluated using four metrics: recall, precision, accuracy, and the F1-score, calculated for each LULC class across different time periods. Recall measures the model’s ability to correctly identify positive cases, accuracy reflects the overall correctness of the classifications, and the F1-score provides a balance between recall and precision, particularly in cases where there is class imbalance or costly false predictions. As detailed in Figure 6, Forest (primary and secondary), along with small-scale cultivation, consistently exhibited high recall, precision, accuracy, and F1-scores throughout the study period. Shrubland demonstrated moderate to high performance across the years, although with some variation in its metrics. Teak plantations, although generally classified with high accuracy, showed moderate recall. Other LULC categories, such as tea plantations, sisal plantations, spice agroforestry, and citrus agroforestry, exhibited more variable performance across the time slices. Built-up areas, rocks, grasslands, and surface water consistently achieved high accuracy, although their recall and F1-scores varied.

Significant LULC changes were observed in the foothills and escarpment regions, while changes in the upland areas were relatively minor. The primary changes involved deforestation and the expansion of small-scale cultivation. Earlier studies in the Sigi River catchment have noted that substantial LULC changes in the upland occurred prior to 1983, beginning in the early 1900s with the introduction of coffee and tea cultivation by German and British colonizers, respectively (Milne, 1937; Hyytitiinen 1995). After World War I, the British replaced poorly performing coffee with tea, which then became the dominant crop in the upland areas (Milne, 1937; Hamilton and Mwasha, 1989). In addition, Sikh Sawmills, a subsidiary of the Tanzanian Wood Industry Corporation, conducted logging operations in the East Usambara Mountains until the mid-1980s. These activities were criticized for causing habitat destruction and damaging non-harvested trees, leading to severe ecological impacts, including significant soil erosion (Newmark, 2002). However, from 1983 onwards, LULC changes in the upland have been less pronounced, likely due to the implementation of various regulations and external interventions, such as the designation of nature protection areas and the promotion of sustainable agricultural practices.

Due to the high rate of deforestation in the past, significant efforts have been made to rehabilitate and conserve the area through the establishment of protected areas conservation agriculture since 1988. The gazettement of protected areas, such as Amani, Kambai and Nilo Nature Reserves, has restricted local communities from encroaching on forestlands, thereby reducing the rate of deforestation and expansion of small-scale cultivation (Kessy, 1998; Hall et al., 2014). Newmark (2002) and Hall et al. (2014) observed that most reserve adjacent to protected areas did not experience high rates of deforestation or small-scale cultivation expansion, particularly in the upland and headwater areas of the Sigi River catchment. However, the foothills and escarpment experienced a higher rate of deforestation and small-scale cultivation expansion, as much of the land, particularly primary and secondary forests, is public lands. This made forests vulnerable to deforestation, as land can be cleared for cultivation without special permits. Consequently, this has led to a high rate of deforestation in these regions. The presence of protected areas has also strengthened conservation efforts to reduce LULC fragmentation. Over the years, various projects have been launched to counteract trends in LULC change (Hamilton and Bensted-Smith, 1989; Kessy, 1998; Yanda and Munishi, 2007; Hall et al., 2014). However, these projects primarily focused on the upland areas, promoting agroforestry practices to curb the formerly rapid expansion of small-scale cultivation (Bullock et al., 2014). The effects of these efforts are evident in our study, where a distinct portion of small-scale cultivation areas in the upland have been transformed into spice agroforestry, indicating that the goals of these conservation projects have been partially achieved. On the other hand, there are reports that cardamom agroforestry, in particular, has led to encroachment on natural forests (Reyes et al., 2006).

Population growth appears to be a significant driver of the observed land use and land cover (LULC) changes in the lowland, foothill and escarpment areas. These regions, with higher population densities than the upland areas, face considerable pressure on arable land, which has predominantly been converted from primary and secondary forests to small-scale cultivation. The expansion of small-scale agriculture in these areas is largely driven by population growth and the increased demand for land to sustain livelihoods (Bwagalilo et al., 2015; Basche and DeLonge, 2019; Clement et al., 2021). In this context, the frequency of man-made fires used to clear land for cultivation has also risen. Newmark (2002) reported that in 1998, over 600 ha of forest (<600 m asl) in the East Usambara Mountains were destroyed by fires that spread from nearby agricultural fields into adjacent forested areas. Additionally, much of shrubland, which is often abandoned after cultivation, has been converted into small-scale agricultural land. These areas are relatively easy to clear due to the absence of strict regulations governing land development. As a result, the rate of small-scale cultivation expansion is higher in the lowland, foothill and escarpment areas compared to the upland regions.

External market demands on agricultural products and government policies have played a significant role in shaping the differential land use and land cover (LULC) changes across elevation zones. Crops are often cultivated in response to high economic returns (Kihiyo, 1992; Kimaro et al., 1994; Bayliss, 2008). The decline in global markets for certain commercial crops such as sisal and tea, combined with shifts in government policies towards privatization, have led to increased unemployment. In response, local communities have turned to alternative income sources, often encroaching on forested areas to cultivate high-return crops (Kihiyo, 1992; Reyes et al., 2006; Kimaro et al., 1994). This dynamic has driven significant changes in land use and land cover (LULC) in the lowland, foothill, and escarpment areas, which are characterized by large areas of public land, less protected forests, and easier accessibility compared to the uplands. Our study found a notable decrease in sisal cultivation following the shift from communal farming to private ownership in 1996 (Figure 8). This transition coincided with substantial forest loss and an increase in small-scale agriculture in both the foothill and escarpment regions.

Figure 10 illustrates that in foothill, lowland, and escarpment areas, forest (both primary and secondary) were first converted to shrubland before being transformed into small-scale cultivation. It appears that shrubland acted as a transitional land cover during the deforestation process from forest (primary and secondary) to small-scale agriculture. In contrast, the upland areas followed a different path, where forest (primary and secondary) was directly converted to small-scale cultivation without passing through a transitional land cover. This distinction underscores the varying dynamics of land transformation across different elevations, driven by a range of factors. For example, in the foothill, lowland, and escarpment areas, the transition typically begins with the opening of primary and secondary forestland through activities like timber harvesting, charcoal production, and fuelwood collection. As the forest becomes degraded, the land is gradually repurposed for agricultural crops, and over time, the remaining trees are cleared, making way for small-scale cultivation (Newmark, 2002; Yanda and Munishi, 2007). In contrast, in the upland areas, forest (primary and secondary) is initially encroached upon, with the understory cleared to make space for the cultivation of spice crops like cardamom (cf. Reyes et al., 2006; Reyes et al., 2009).

The pronounced loss of both primary and secondary forests on steep slopes may be attributed to the earlier clearing of forests in flat and gently sloping areas before 1983. Population growth in recent decades has driven settlements into steeper terrain, necessitating agricultural expansion in these areas. Additionally, land tenure patterns have played a critical role in shaping land use dynamics. While gentle and moderately sloping areas are predominantly designated as forest reserves for government teak plantations, cash crop estates, and large-scale citrus agroforestry, village-owned lands on steeper slopes have become the primary areas for small-scale cultivation (Yanda and Munishi, 2007). Local households have increasingly encroached on steeply sloped forestlands due to the fertility of these soils, which offer higher organic matter and nutrient content, enhancing crop productivity (Hamilton and Bensted-Smith, 1989; Mwanyoka, 2005; Yanda and Munishi, 2007; Reyes et al., 2009). Notably, our findings indicate that upland areas with steep slopes exhibit the highest coverage of small-scale agriculture compared to other slope classes. This trend may be linked to historical land-use patterns, where tea plantations established primarily by the British after World War I, often replacing abandoned German coffee plantations were preferentially located on gentler slopes (Milne, 1937; Hamilton and Mwasha, 1989). In contrast, in the foothill and lowland areas, small-scale cultivation is more concentrated on gentle slopes, likely due to greater accessibility, ease of management, and technological applications (Newmark, 2002; Yanda and Munishi, 2007). Our results highlight that steep slopes have experienced significant LULC transitions, primarily through the conversion of forests to shrubland then small-scale agriculture. Given the susceptibility of these areas to soil erosion and water resource degradation, they should be prioritized for targeted soil and water conservation measures. While afforestation efforts have been implemented in the foothills, agroforestry may offer a more sustainable approach for steep slope areas, mitigating erosion while supporting local livelihoods.

The performance of our models in identifying various land-use/land-cover (LULC) classes reveals important insights into the strengths and limitations of the classification process. Our results show that the models effectively distinguished between forestland (both primary and secondary) and small-scale cultivation yet faced challenges in accurately classifying shrubland and teak plantations. This discrepancy likely arises from the visual similarities between shrubland, small-scale cultivation, and young teak plantations, which complicates the differentiation process and can lead to reduced classification precision (cf. Zhao et al., 2014; Feng et al., 2022). Furthermore, the varying performance metrics for more complex LULC categories, such as tea plantation, sisal plantations, spice agroforestry, and citrus agroforestry reflect the inherent diversity and complexity within these land-use systems. These classes represent a wide range of vegetation types, management practices, and environmental conditions, contributing to significant variability in their spectral signatures and spatial patterns. The diverse nature of these land-use systems presents additional challenges for accurate classification, highlighting the need for more refined techniques and additional data sources. To address these challenges and enhance classification accuracy, we incorporated high-resolution, multi-temporal imagery, complemented by field observations. His integration not only improved the model’s ability to differentiate between similar land-use types but also provided valuable insights into the nuanced characteristics of these complex land-use systems (cf. Feng et al., 2022). These strategies proved effective in refining the model’s ability to differentiate between land-use types and underscored the importance of leveraging both remote sensing data and ground-truth observations in capturing the diverse dynamics of land-use change.