Predicting land use/cover changes and its association to agricultural production on the slopes of Mount Kilimanjaro, Tanzania

ABSTRACT Increasing demand for food production results in Land use and land cover (LULC) changes, which afflicts the provision of ecosystem services in high mountain areas. This work used time-series LULC and selected spatial metrics to predict the LULC changes for Kikafu-Weruweru-Karanga (KWK) watershed (on the southern slopes of Mt. Kilimanjaro) for the next decade. LULC maps were generated by classifying time-series satellite images. We further predicted the implications for selected staple crop production over the next decade. The simulated LULC shows expansion in built-up (by 32.55%/27.04 km2) and agriculture (by 39.52%/52.0 km2) areas from 2018 to 2030. These results suggest that urbanization is likely the next biggest threat to water availability and food production. Grasslands and wetlands are expected to decrease by 57.24% and 39.29%, respectively. The forest area is likely to shrink by 6.37%, about 9.82 km2, and 1.26 km2 being converted to agriculture and built-up areas, respectively. However, expansion in agricultural land shows very little increase in staple food crop production records, suggesting that farm size plays a minor role in increasing crop production. Predicting the near future LULC around KWK is useful for evaluating the likelihood of achieving development and conservation targets that are set locally, nationally and internationally.


Introduction
The global and regional expansion in agricultural land reflects the intensified pressure over the land, putting uncertainties in food production at all levels (Mora et al. 2020). This pressure creates a global challenge of increasing food security while ensuring environmental sustainability, increasing pressure to alteration in land use and cover (Lobell and Field 2007), hydrological cycles (Schlesinger and Jasechko 2014), loss of biodiversity (Newbold et al. 2015), and soil degradation (Smith et al. 2016). These are global issues with local effects, the land use decisions are due to and affected by local, regional, and global driving forces. However, the driving decisions are driven by forces from local (McCusker and Carr 2006), regional (van Vliet et al. 2015) and global level (Lambin and Meyfroidt 2011).
Land-use decisions and food security are determined by the changes in the global production of major crops (Lobell and Field 2007). In many parts of the world, drivers of land use/cover changes are those linked to anthropogenic activities and the natural ecological progressions (Liu et al. 2014). Similarly, mountains of tropical East Africa suffer pressure from growing human populations, opening up of more land for agriculture, and the misuse of natural resources; leading to remarkable changes in its landscapes (Nogués-Bravo et al. 2008). Mostly, LULC changes results from direct and indirect consequences of anthropogenic pushing for natural resources utilization (Misana, Sokoni, and Mbonile 2012). Thus, the growing population among other drivers has increased the vulnerability of the tropical mountains to LULC changes (Kidane, Stahlmann, and Beierkuhnlein 2012). Other vital drivers occurring together with population pressure are poverty and government policies, resulting in significant losses, especially in tropical ecosystems (t Mannetje et al. 2008).
Mount Kilimanjaro is a unique ecosystem harbouring natural vegetation, home gardens, agroforest areas dominated by native trees and savannah in the lowlands. Agriculture is the primary activity that influences the livelihood prosperity of communities in villages and suburbs on Mt. Kilimanjaro slopes (Munishi, Lema, and Ndakidemi 2015). Reports show rapid expansion in agricultural land along the mountain Kilimanjaro slopes Sokoni 2003, Soini 2005b;Misana, Sokoni, and Mbonile 2012). However, the efforts of the farming community seem to be focused on expanding the farm area aiming at increasing food production. This practice is a primary factor driving the conversion of most of the other forms of land use to agricultural land. The natural vegetation cover on the southern slopes of Mount Kilimanjaro is a cornerstone of conservation efforts, with their unique biodiversity and ecosystem, doing their best performance in sustaining ecosystem services (e.g. water, energy and food production) to increasing local population.
The past land-use changes may be used to establish a model to predict the future land-use change trend (Liping, Yujun, and Saeed 2018). LU change models may be categorized into; Empirical and statistical models, e.g., Markov chains and Regression model. Further, dynamic models, e.g. Cellular Automata (CA) model, and integrated model, e.g. integrated Grey wolf optimizer (GWO) and Cellular Automata (CA) model (Cao et al. 2019a), and Conversion of Land Use and its Effects model (CLUE-S) (Guan et al. 2011) may be used. Dynamic and integrated models such as CA are best suited for future LU change prediction.
Several LULC change studies have been done on different spatial and temporal scales on the southern slopes of Mt Kilimanjaro (Mbonile, Misana, and Sokoni 2003;Fairman et al. 2011;Misana, Sokoni, and Mbonile 2012;Mmbaga, Munishi, and Treydte 2017). Conversion of different LULC types to agricultural land on the southern slopes of Mount Kilimanjaro is both historical and a continuous phenomenon. Forest degradation on the southern slopes increased 1,606 ha to 5,170 ha between 1973 and 2000. Further, the land under agroforestry decreased by 25%, while that under annual crops increased by 41% (Mbonile 2005). Soini (2005) found that about 49.97 km 2 (32.8%) of shrubs and bushland were converted to agriculture and other uses from 1961 to year 2000 in Southern slopes of Kilimanjaro. Around the same time, there was a conversion of about 39.5% of bushland to agriculture (Mbonile, Misana, and Sokoni 2003) while the same area experienced a decrease in closed and open forests by 56% and 64%, respectively (Munishi, Hermegast, and Mbilinyi 2009), with a large part of this being converted to seasonal cropland. During this period, grasslands in the area increased by 116%, and riverine vegetation decreased by 53% (Munishi, Hermegast, and Mbilinyi 2009). However, none of the previous studies considered projecting the future LULC scenario in the area. This conversion prompts the need to establish a link between increased agricultural land and food production. The continued land conversions to agriculture by the local communities have been on the increase, with assumptions that the more land (from other forms such as a forest) is converted to agriculture land, the higher the production. However, the assessment of the effectiveness of this conversion by relating it to actual yield and their future projection in sustaining production of the main staple food crops remains poorly documented. Here, we address this gap using a combination of Cellular Automata and Markov Chain (CA-Markov) analysis approach by modelling and predicting the spatiotemporal patterns of LULC using time-series satellite images. The results from this were further compared with the crop production data of major crops recorded from the area as the basis for future LULC prediction to agricultural productivity on the southern slopes of Mount Kilimanjaro. CA-Markov model is reported as a suitable for land cover change detection and simulations due to its ability to account for spatial and temporal components of land cover dynamics (Behera et al. 2012;Baysal 2013). Further, CA-Markov is very efficient in simulating multiple land cover and complex patterns, high effectiveness in datascarce areas and simple in calibration (Memarian et al. 2012).
The objectives of this paper are therefore three folds (i) to map LULC in three historical periods (1993, 2006 and 2018) on the southern slopes of Mt. Kilimanjaro, (ii) to project LULC in 2030 and to evaluate the differences from the current LULC conditions, and finally (iii) to study the implication of the LULC changes on the selected staple food production. The information from this study is expected to provide useful information that can help to advance our policies and practices, as well as to inform planners on the extent to which sustainable development goals can be realized in future. Furthermore, the information generated can further inform the action to accelerate progress towards the achievement of the targets of Sustainable Development Goals (SDGs) and Paris Agreement (COP 21).

Study area
This study was carried out on Kikafu, Weruweru, and Karanga rivers catchment (KWK catchment) on the southern slopes of Mount (Mt.) Kilimanjaro. The watershed originates from the thick mountain forest on the southern slopes of Mount Kilimanjaro (Figure 1). The catchment joins the Meru river networks before entering the Nyumba ya Mungu dam. The area is mostly under rainfed and few irrigated agriculture practices. This zone experiences an increasing population since the colonial times (NBS 2012(NBS , 2018. The population density is almost 300 people per square Kilometre (Mbonile, Misana, and Sokoni 2003). Rainfall follows a bimodal regime and is influenced by trade winds (Coutts 1969), which are controlled by Indian Ocean Dipole (IOD) and meridional displacement of Inter-tropical Convergence Zone (ITCZ) (Otte, Detsch et al. 2017b). There are two distinct wet seasons in the hydrological year, namely Masika (March-May) and Vuli (October-December). One short dry period occurs in January-February and June-September (Appelhans et al. 2016, Otte, Detsch et al. 2017. The mean annual precipitation sums to more than 2,500 mm in the southern montane forest belt. The northern mountainside is comparatively dry and receives less than 1,000 mm per year (Hemp 2005). In both mountainsides, the mean annual rainfall increases with elevation.
Mt. Kilimanjaro slopes in northern Tanzania offer fertile volcanic soils and promising weather condition for agricultural production. Thus, its ecosystem supports about three-quarters of the population in the region (Misana, Sokoni, and Mbonile 2012). Agriculture production and its expansion have long been a critical determinant of land-use change, and settlement on the lowlands have contributed significantly to change in land use, ecosystem functions and environmental changes on the slopes of Mt. Kilimanjaro (Peters et al. 2019). Hence, the need to understand land-use change dynamics is an essential requirement for the global assessment of mountain landscape dynamics (Mottet et al. 2006).
Cultivation on Mt. Kilimanjaro slopes, mostly the southern slopes, started when people inhabited the slopes approximately more than 2000 years ago (Masao 1974). From there, there has been a fast expansion in land under agriculture and changes in LULC from time to time (Agrawala et al. 2003; Misana, Majule, and   Chiwa 2012;Mmbaga, Munishi, and Treydte 2017). The common staple food cultivated in the area includes banana, maize, beans, horticultural crops, among others (Mulangu and Kraybill 2013). About 28% of the Chagga home gardens in the highland areas are connected to the irrigation network (Mulangu and Kraybill 2013). In connection to agriculture, population growth and village fragmentation have increased the need for more land for cultivation, residential and infrastructures on the upper slopes, which prompted opening up more land for agriculture in the lowlands.
On Mt. Kilimanjaro slopes, the rapid conversion of bush and grassland to agricultural land for food production is evident in all areas on the southern slopes except in areas with shallow soils on the volcanic hills (Misana, Sokoni, and Mbonile 2012). Overall, changes in land cover and use have increasingly affected microclimate, land management, and biodiversity conservation in the region (Soini 2005, Hemp 2006). These changes affect watershed moisture balance in most of the places (Behera et al. 2012), and hydrologic processes (Fohrer et al. 2001;Mulangu and Kraybill 2013). Hence, it is essential to understand the past vegetation dynamics and their possible future evolution for biodiversity conservation and management in mountainous areas as also previously argued by Tattoni, Ciolli, and Ferretti (2011).

Data acquisition and image pre-processing
The past land-use data (1993, 2006, and 2018) of the catchment were obtained by classifying their respective Landsat (4, 7, and 8) level 1 images. The satellite images were accessible freely from the Global Visualization Viewer of the United States Geological Survey (USGS) website (https: glovis.usgs.gov/). The Landsat images used for land use/cover classification in this study were the Operation Land Imager (OLI) for 2018, Landsat 7 for 2006, and Landsat 4 images for 1993. All images used in this study were of the same season, aiming at minimizing seasonal variability. All images were resampled to a cell size of 30 m x 30 m and also corrected for atmospheric and geometric errors. For the purpose of matching these images with other spatial data of the study area including the ground truthing coordinates for accuracy assessment, their projections were transformed to local projections of the study area spatial data, as previously described by Liang et al. (2018). Individual bands of the images were extracted, and band composites were made before the extraction of the required study area using the polygon boundary. Due to the nature of the study area, the upper part is covered by Mount Kilimanjaro and the thick forest of Kilimanjaro National park. This part is covered by the clouds for most of the time, especially during wet seasons. Hence, it was challenging to get cloud-free images, specifically for 2018 during the wet season. Thus, to overcome this data quality challenge, the cloudy parts of the image (which is mostly rock and snow which do not change significantly over time) were sliced out. The sliced out locations were replaced by the cloud-free observations which were extracted from the satellite image of 23 September 2016 by using the Mosaic function in ArcMap. The Landsat 7 scan line errors were corrected by infilling the missing data with the scan line correction data using 'fill no data' function in Quantum GIS (QGIS) before the images were transferred to ArcGIS for further processing.
In addition to satellite images, other data used in this study include Digital Elevation Model (DEM), Road and stream networks as well as the National Parks map. These data were obtained from different sources. Table  1 summarizes the data types used and their characteristics.

Image classification and accuracy assessment
Land cover was classified into ten (10) major classes (Table 2), which represent both natural protected and non-protected areas. This means both changes due to natural drivers and human influence were captured. Prior to Image classification, ground-truthing points guided the expert to make an on-screen selection of the areas to pick pixels to be used as the training samples. The ground truth points and Google Earth images overlays were used to identify ambiguous features during training of the past (1993 and 2006) images, thus ability to label the training sample and give a specific land class. Then, the signature file was developed from the training samples. The images were then classified  1993, 2006, and 2018 land cover maps of the area. Since supervised classification resulted in some pixels placed in wrong spatial areas; for example, agricultural land on top of Mt. Kilimanjaro, expert knowledge was used to eliminate these wrongly classified features. In this study, a raster calculator function was used to guide at which elevation a particular landuse class should appear. This was based on the expert knowledge of the study area, thus called expert classification. The details of this procedure are as described by Taweesuk and Thammapala (2005). Expert classification is used to integrate remote-sensed data with other sources of georeferenced information such as DEM, land-use data, and spatial texture, to improve classification accuracy. An accuracy assessment was done to ascertain the level of the agreement of the classified images with the ground features. Ground truthing points and other reference points which were identified from the features of the satellite and Google Earth images were used. The examination of the hard to access areas (Remote montane forest, Kilimanjaro National Park, and high mountain areas) and training samples for the past (1993 and 2006) was achieved through the use of Google Earth. ArcGIS 10.6 student version and IDRISI Selva v17.0 were the image processing packages used in this study. Also, interviews with local people were conducted aiming at improving the classification accuracy of historical images (1993,2006). Field surveys were done aiming to collect ground-truthing points using a hand-held GPS receiver. The ground truth points were used together with expert knowledge of the study area for LULC classification accuracy assessment. A total range of 50 to 70 ground-truthing points was collected during fieldwork and was used for accuracy assessment of each land-use class. Both producer's, users and overall accuracy and Kappa coefficient were calculated. Producer accuracy was calculated as the number of samples identified correctly divided by the reference totals, whereas user's accuracy was calculated as the number of samples identified correctly in each individual class divide by the classified totals, respectively (Foody 2008). The overall classification accuracy of each image was obtained as a total number of correctly classified pixels divide by the total number of sample points used in accuracy assessment (Foody 2008). Equation 1 is used to calculate the coefficient of agreement, that is Kappa statistics. Kappa coefficient is a measure of the total statistical accuracy of an error matrix between the classified map and the reference data; this coefficient takes non-diagonal elements into account. The procedure for image classification and accuracy assessment is summarized in Figure 2.
where n is the number of rows, x ii is the number of points correctly classified for a particular land-use category, x i+ and x +i are marginal totals for raw and column i associated with the category, and N is the total number of observations in the entire error matrix (Bishop, Feinberg, and Holland 1975).

Geospatial modelling framework
Generally, the main procedures involved to achieve the objectives of this study include; (a) land-cover mapping of 1993, 2006 and 2018 through the classification of satellite images (Maximum likelihood). (b) Computation of transition area matrix obtained from the Markov process, indicating the number of pixels that is likely to change from each land-cover class to another class over a specified time interval (1993--2006  Built-up area (BU) Tarmac and gravel roads, concrete areas, urban and rural settlements 2 Agricultural land (AGR) All land with crops 3 Water (WAT) Rivers, water in wetlands, fish ponds and water in agricultural areas 4 Forest (FOR) Areas with closed trees, thick canopy, both natural and planted forests. 5 Bare land (BAR) Bare soils, quarry areas, sands, and eroded soils 6 Grassland (GRA) Tall to short grasses, sometimes bare soils in the dry season. 7 Wetland (WET) Areas with lands partially covered with water and grasses 8 Shrubland (SHB) Areas with wood trees but a less closed canopy 9 Rocky surface (ROS) Areas covered with bare rock and less vegetation 10 Glacier  Barredo et al. (2003) described the factors involved in the simulation of land-use change. A factor is a criterion that influences the suitability of the particular land cover in a particular area. Factors vary in values; however, landuse cover change projections are usually stretched from 0 to 255 (Musa, Hashim, and Reba 2018). These factors are grouped into five groups, namely; Environmental factors, i.e., constrained areas, local-scale neighbourhood, i.e., one land use affected by the other, spatial characteristics of the cities, i.e., distance and access to the city centre. Others are land-use zoning status, i.e., urban planning policies and factors related to individual preferences, social-economic growth, and political system such as population growth dynamics, i.e., population growth. These factors were assessed to monitor their behaviour in the study area; however, the future changes in urban plans and policies were not included in the current model. The land-use class suitability maps are essential and play a central role in CA-Markov modelling because they drive the spatial patterns of the simulations. Suitability images infer the suitability of each cell for a particular land cover. Thus, their preparation requires thorough knowledge and technique as the other steps in  modelling (Tattoni, Ciolli, and Ferretti 2011). In this study, the factors and constraints for the year 2018 were used to develop the suitability maps. Factors and constraints were developed for each land-use class ( Table 3). The method used to prepare the factors and constraints has been well described by Hyandye and Martz (2017). The Multi-Criteria Evaluation (MCE) function was used to integrate the information from various criteria to arrive at the individual evaluation index (suitability maps) (El-Hallaq and Habboub 2014).

Criteria development and generation of land-use suitability maps
Factor maps were prepared in ArcGIS v.10.6 software and then imported to IDRISI Selva v.17 software for further processing. All constraint maps were prepared in Idrisi Selva by editing the existing values in the feature class (es) and assigning zero to the targeted constraint of the feature class. Vector layers of proximity to highways, and waters, accessibility to existing settlements, and forest reserves were obtained from Tanzania National Roads Agency, Pangani Basin Water Office (PBWO), classified satellite images and Tanzania National Parks GIS unit, respectively. The slope map was created from the DEM data in ArcGIS before importing it to Idrisi Selva. The annual precipitation records for the stations in the study area were obtained from the Tanzania Meteorological Agency (TMA), while soil types were generated from the soil data collected from PBWO. Population density map was generated using the population figures projected from the 2012 National Population Census using the formula in equation 2.
where N t is the number of people at a further date, P is the present population, e = is the natural logarithm base of 2.71828, r = rate of increase divided by 100, and t = the time.
All the vector files were rasterized in ArcGIS, and imported into IDRISI Selva software. Using a decision support tool in IDRISI Selva Software, factor maps were standardized and assigned a continuous suitability scale from 0 to 255, where 0 represents unsuitable sites, and 255 represents the most suitable sites. This was done using different fuzzy memberships and control points (Table 3). Constraints are the variables that refer to the restricted areas for development; thus, there are no medium values. The values are either 0 (not suitable) or 1 (suitable) (Eastman 2012). The constraints were assigned as Boolean image characters of 0 and 1, in these characters; 0 was assigned to areas not suitable for future development, and 1 was assigned to areas suitable for future development. The control points were used to evaluate the degree of membership of a pixel in a factors layer; these are set considering the association of a factor to urban growth (Eastman 2012).
The fuzzy membership function used in this study was Sigmoidal, J-shaped, and Linear membership functions. These usually decide the shape and trend (Decreasing, increasing, or symmetric) (Eastman 2009(Eastman , 2012. The trend is said to be 'monotonically increasing' when the factor value increases and the suitability of the particular land cover increases. It is said to be 'monotonically decreasing' when the value of the factor increases and the suitability of the particular land cover class decreases. Also, it is said to be 'symmetric' when the value of the factor increases, and the suitability of a given land cover class increases from a given point, attains maximum suitability at a given point, and then begins decreasing. Standardized factors and constraints for each land cover category were combined to produce suitability maps ( Figure 3). Weights were assigned to each factor to indicate their importance concerning one another using Analytical Hierarchy Process (AHP) employing pairwise comparisons. Other approaches can be user defined or assign equal weights (Eastman 2012). The importance of each factor was evaluated by assigning values with the ranging from 1/9 (less important) to 9 (more important). Then, the principal eigenvector was employed to generate the final weight of each factor. A consistency threshold of less than 0.10 was acceptable for evaluation of the pairwise matrix, and re-evaluation of the pairwise matrix was done for anything above the threshold (Lin, Kou, and Ergu 2014).
Then, the MCE module was employed to aggregate all the factor layers to form the suitability maps for each land-cover class (Figure 4). In this study, Weighted Linear Combination (WLC) was used (Figure 3). The WLC employs linear functionality to overlay standardized factors according to their weight of importance, where total factor weights sum to 1 (Eastman 2012).

Prediction of future land use/cover changes and the link with crop production in KWK watershed
The earlier land cover image (t-1) of the 1993 and the current land-cover image (t = 1) of 2006 were used to run a Markov model, which generated two files, transition probability, and transition areas file. The transition probability indicates the likelihood of a pixel to change to a different land use/cover (or remained constant) in the next time. In contrast, a transition area matrix shows a total area with the likelihood of changing the next time (Eastman 2012). In Table 5, the categories of the 1993 (t-1) are shown in rows, whereas the column shows those of 2006 (t = 1). The outputs of the Markov model combined with the suitability maps were used to run the CA-Markov model using a contiguity filter of 5 × 5. The output of this model is the simulated map of the year 2030.
The model validation is an integral part of the modelling process (Memarian et al. 2012). Literature provides a wide range of the methods available for performing model validation; some of these methods are the kappa coefficient and Cramer's V (Baysal 2013). The quantity and allocation disagreements (Brown et al. 2013), chisquare, and F-test of the observed and simulated images (Katana, Ucakuwun, and Munyao 2013), as well as the relative operating characteristics (Pontius and Schneider 2001), are among the available validation methods. This study used quantity disagreement and allocation disagreements through the VALIDATE function embedded in IDRISI selva software. The kappa coefficient of 0.61 and above was considered acceptable agreement level of the two land cover maps compared (simulated and observed) (Eastman 2012). The model was considered suitable for further processing and used to simulate the future land cover map of the year 2030 after achieving the minimum kappa coefficient. According to Eastman (2012), this method entails analysis of both the level of agreements in quantity and the location of cells between the simulated and observed land cover maps through statistical algorithms. Integration of the model results with agriculture production for the selected staple food crops is a crucial procedure to measure the feasibility of expansion of cultivated land for production. Although it was challenging to separate agriculture land into LULC for each of the selected crops due to high fragmentation and complexity of the farming system practiced in the study area, secondary data were used. These data were sourced from the agriculture division, Kilimanjaro region whereby they provided the estimates of the available cultivated land and crop yield for each crop. Thus, the association of the crop production and agriculture land from the model output was done by computing the proportion of land size for each crop, from the known proportion of arable land from the long-term records of agricultural production from the agriculture department of Kilimanjaro region.

Accuracy and model validation results
The overall accuracy assessment of the supervised images for the years 1993, 2006, and 2018 were 85.32%, 87.23%, and 88.77%, respectively (Table 4). These values satisfy the minimum accuracy threshold of 85% needed for efficient and authentic land use/ cover change analysis and modelling (Araya and Cabral 2010;Ahmed, Ahmed, and Zhu 2013). The results of this study are regarded as acceptable because the accuracy values are higher than 80%, as stated by Jensen (1996).
The total agreements of 0.7286 between the simulated and observed maps of 2018 (κstandard) suggest that there are 72.86% better agreements between the two pair maps of 2018 than by chance. High values for the overall simulation accuracy (κ-no) and κ-location were 78.78% and 73.86%, respectively. The disagreements due to quantity and allocation were low (0.1283 and 0.1431, respectively) almost with similar magnitude.
The VALIDATE module in IDRISI Selva generated the kappa coefficients (κ), which were used to assess the predictive power of the simulation model. The overall simulation accuracy (κ-no) was found to be 78.78%; the model's ability to accurately specify the location (κlocation) was found to be 73.86%. Table 6 depicts the agreement and disagreement measures between the observed and simulated LULC of 2018. The total agreement between the simulated and observed maps (κstandard) was found to be 0.7286. The quantity and allocation disagreements were low (12.83% and 14.31%, respectively).

Patterns of land use/cover changes in the KWK watershed
Classified land use/cover image maps of KWK catchment are presented in Figure 5. The results of the analysis of the historical land use/cover changes selected spatial area gain/loss are summarized in Table 5. Generally, the land use/cover change results (Table 5) show the increasing trend in the built-up area from 1993 to 2018,

Projection and validation of land use/cover change on the slopes of Kilimanjaro
The model's ability to accurately simulate the year 2030 land use/cover maps were validated by comparing the observed and simulated land use/ cover maps of 2018 ( Figure 6). There is an acceptable degree of resemblance between the two maps in terms of the spatial distribution of all land-use categories. The corresponding land use/cover classes in Figure 6(d, e) demonstrated a good match in terms of spatial distribution with a minor adjustment in the intensity and alignment.

The transition matrix of land-use/cover change on the slopes of Kilimanjaro
The Markov model generated the probability transition matrix for land use/cover maps of 1993 and 2006, as shown in Table 7. The transition matrix areas show the total cells to be transformed into another category in the next time span (Eastman 2012), it is a stochastic matrix whose (i, j) entry gives the probability that an element moves from the jth category to the ith category during the next step of the process (Andrilli and Hecker 2010). This probability transition matrix was used to simulate LULC for the year 2018. In Table 7, the rows show the land use/cover classes for the earlier date t-1 (1993) whereas, the columns show the probability of the landuse/cover classes of the current date, t = 1 (2006).   Most of the LULC categories were dynamic. However, the probability of most classes to remain unchanged was high for rock surface (0.9791), shrubland (0.7818), and forest (0.6537); agriculture showed moderate probability to remain unchanged (0.4010). Other classes showed a low probability of remaining unstable in the future. The probability transition values in Table 7 showed that there was no possibility that other lands use classes could change to water.
Further, there was a higher probability of wetland to change to grassland (0.3798) and agriculture (0.3639) than agriculture and grassland to change to the wetland. There is also a high probability that water (0.4134) will change to agriculture and 0.2807 chance of water to turn to grassland. The probability transition values showed that the chance for other land-use classes to change to glaciers is almost negligible, except high probability of glaciers to disappear and replaced by rocky surface (0.7460).

Projected land use/cover changes on the slopes of Kilimanjaro
The simulated land cover map of the year 2030 from the CA-Markov model is presented in Figure 7. The total areas of LULC classes of 2030 were compared with the ones of the actual land cover map of the year 2018 to ascertain the quantitative changes in land use (Table 8).
The results in Table 8 show that under the current scenario for the period 2018-2030, the built-up area will expand by 32.55% (29.25 km 2 ), whereas agricultural land is expected to expand by 39.52% (64.18 km 2 ). Of all builtup areas, 10.52, 5.34, and 13.22 km 2 will have been converted from agriculture, forest, and barren land, respectively, whereas a large proportion of grassland (53.92%/57.16 km 2 ) and forest (7.46%/9.82 km 2 ) will be converted to agriculture. Glacier ice is also expected to shrink by 18.79% (0.31 km 2 ). The model predicts that by  2030 there will be a significant decrease in grassland, wetland, and barren land by 57.24%, 39.29%, and 15.8%, respectively.

Food production patterns in relation to expansion in agricultural land
Maize, rice, banana, and beans are everyday staples and cash crops along the Mt. Kilimanjaro slopes. The selected crops are among the six major crops in the world (Lobell and Field 2007). Integrating the expansion in agricultural land results observed in this study to averaged food production statistics in the five districts (Rombo, Hai, Siha, Moshi district council, and Moshi urban district councils) of the study area 2016/17 reported by URT (2018); the maize, rice, banana, and beans production were expected to increase by 28.33% from the current records. However, the available record does not reflect the increment in crops production as the agricultural land increases, as observed in Figure 8(a). The observed trend in Figure 8(a) for banana shows a decreasing trend. The sharp decrease in banana yield is observed between 2005 2008, and continue to decrease up to 2017. A similar trend is observed in Maize, rice, banana, and beans, with production less than 200,000 tons despite the increment in agricultural land described in section 3.2 and Figure 8(b).

Patterns of land use/cover changes
The increase in the built-up area (Figure 8) reflects the growing population trend. Studies show that future growth in the population in the catchments goes together with the urbanization of more land, including areas along sideways, the road network, and urban suburbs (Hyandye and Martz 2017). The evidence of this is the emergence of the Siha district; studies show that many new villages have mushroomed along the road on the slopes of Mt. Kilimanjaro (Soini 2005). It is also worth to say that the growing population goes together with the food demand, hence expansion in agricultural production. From 1993 to 2018, there was a small proportional (7.04%) loss of the forest area initially present in 1993 (2.28 km 2 ). This may be due to the expansion of settlements close to the forest buffer zone area and shrinkage of the riverine forest due to the expansion of agriculture land. The encroachment into the forest reserve and replacement of the natural vegetation by the farmed crops was also reported in some areas of Mt. Kilimanjaro slopes (Mbonile, Misana, and Sokoni 2003). The results from this study complement the idea to revisit the implementation of laws and regulations governing protected areas. Evidence on land-use change especially on the boundaries  of protected areas from other countries suggests directions for research and the need for policymakers to redesign the policies and their implementation in protected areas (Tesfaw et al. 2018). Shrinkage of the wetlands by 1.62% may be due to increased drought, especially in lowland areas. Drought is one of the primary reasons for the invasion of wetlands by farmers and pastoralists (Mbonile, Misana, and Sokoni 2003). Similar studies on Mt. Kilimanjaro slopes reported comparable LULCC trend. Mainly conversion of by 53 km 2 of shrub and grassland and sparse vegetation to maize and vegetable cultivation from 1987 to 2005 in Kahe plains (GITEC 2011) increased forest degradation on the southern slope from 1,606 ha to 5,170 ha between 1973 and 2000 (Mbonile 2005), conversion of about 39.5% of bushland to agriculture between 1973 and 2000 years (Mbonile, Misana, and Sokoni 2003), degradation of more than 41 km 2 of the forest between 1952 and 1982 (Yanda and Shishira 2001), conversion of about 49.97 km 2 (32.8%) of shrubs and bushland to agriculture and other uses from 1961 to 2000 years in Kirua Vunjo division (Soini 2005), and increased cultivated land from 54% (in 1973) to 63% in 2000 on the southern and eastern slopes (Misana, Sokoni, and Mbonile 2012). Replacement of the bushland by cultivated fields is a continual phenomenon. In Kirua Vunjo division on the southern slopes, for example, bushland occupied 40% of the area in the early 1960s, by 2000 this was reduced to 7% through transformation to cultivated fields (Soini 2005).
Clouds patching at high elevation may influence the classified image interpretation and hence discrimination of glacier ice on top of Mt. Kilimanjaro. In the long run, this affects the preparation of training samples hence impacting classification accuracy. In most cases, clouds and cloud shadows obscure the parts of the glacier, snow cover, and glacier surface debris cover (Rastner et al. 2019). Although the results from satellite image classification accuracy pass the minimum requirements of 80%; it is worth mentioning that the interpretation of the glacier ice results in this study might have been affected by this clouds cover challenge.
Historically, highlands were the most populated areas as compared to the lowlands. Expansion in agriculture is reported to have created a new and distinct group of farmers who have settled in the dry lowlands previously considered by the Chagga as unsuitable for permanent settlement due to inadequate rains and malaria (Soini 2005). Regardless of the low rainfall amounts in lowlands, most of the intensive agriculture activities currently are in the lowlands. Periodic floods and running rivers make the lowlands the most productive areas in the season and throughout the year in irrigated areas. Thus, a need to invest in advanced agrarian technologies to increase production per unit area. Also, it is important to focus on alternative activities that will supplement agricultural production (Soini 2005;Mbonile, Misana, and Sokoni 2003).
The growing need for cultivable land increased land scarcity; thus, even the most peripheral areas were cultivated, e.g., cultivation besides the river banks despite the available bylaws. Agriculture expansion took its course to the lowlands with a high rate after the introduction of the agriculture policy in 1983 (Mbonile, Misana, and Sokoni 2003). Annual crops such as maize and horticultural crops replaced perennial crops to a larger extent. From that time, then lowlands and highlands are intensively cultivated.
The changes observed in this study are in a similarity to global reported trends. Reports show that the decrease in the area of tropical savannas by over 50%; mainly due to the alteration into cultivated land (Assessment 2005). About 20% of the Sub-Saharan Africa savanna vegetation was also converted to cropland (UNDP 2000). Also, the Food and Agriculture Organization (FAO) expects a further conversion of natural savanna areas to small-scale subsistence agriculture (Zhou 2011). Further, in upper Suha watershed, Ethiopia; agricultural land, plantation forests and settlement were reported to increase by 18%, 93% and 125%, respectively. The report further shows a decrease in natural forest, shrubland and grassland by 54%, 82% and 16%, respectively (Mekuriaw, Cherinet, and Tsegaye 2019).

Validation of land use/cover simulation
The accuracy values during validation mean that the simulated map did not differ much from the observed map in terms of the quantity and location of cells in each land-use category. The error (disagreements) may have emanated from the transient nature of most of the LULC classes except few, which are persistent and mostly located on the upper part of the Kilimanjaro National park that is under protection by the law. Memarian et al. (2012) reported a comparable high quantity and allocation disagreement. In most cases, errors could be due to the insufficient land-use suitability maps as well as overgeneralization in image classification. Also, the shape of the chosen contiguity filter influences the suitability maps. This could result from inefficient digitization of some of the dataset, and factors and constraints could have influenced the results (Musa, Hashim, and Reba 2018). Another common source of cell-by-cell error is inefficient georeferencing and geometric correction of the comparison map (Eastman 2012). To improve these validation results, high-quality dataset such as satellite images and other spatial data needs to be available. Also, advanced methods need to be integrated into the preparation of more realistic suitability maps; one of the suggested methods is logistic regression. However, the results in this study meet the requirements for use as described by Pontius and Millones (2011), and eligible for use in simulating the future land use/cover for 2030.
Suitability maps can be validated before their use for land-use prediction. Methods such as Relative Operating Characteristic (ROC) may be used. However, ROC does not account for the spatial arrangement of the model's successes and errors. All of the information for the ROC comes from contingency tables that show only cell-bycell analyses of the association between a map of real change and maps of simulated change, where each grid cell is a homogenous category (Pontius and Schneider 2001). There is a spatial-dependent uncertainty in the prediction maps that ROC method is not able to reveal. Thus, ROC requires other statistics such as Cohen's kappa index to reveal the magnitude of the agreement between the maps (Vakhshoori and Zare 2018). Therefore, it requires ROC to be supplemented with visual comparison and additional measures of association that account for the spatial pattern (Costanza 1989). Thus, in this manuscript, the model was validated using the approach proposed by Pontius and Millones (2011), which advocates for quantity and allocation agreement and disagreement for accuracy assessment. However, it is worth stating that suitability maps were modified several times (following the procedure in Figure 3) until the final acceptable Kappa index was achieved.

The probabilities of land use/cover change
The low likelihood of rock surface, shrubland, and forests to change to other classes could be attributed forest, and mostly shrubland is located in the upper part of Kilimanjaro National park which is the conservation area under protection by the law (Soini 2005). Human activities are restricted to tourism and other conservation activities. For agriculture, the notable high probability of remaining unchanged may be due to the reason that some of the lowland areas are practicing irrigation agriculture and that all the images were of the same season. Agriculture activities are historical in the area, and most of the suitable areas have been utilized by the residents mostly residing in the upper elevations. Conversion of more grasslands and barren land areas to agriculture and urbanization denotes the need for more investments in other economic activities apart from agricultural production. The high probability of wetlands to change to grassland and agriculture is probably due to the growth of vegetables and short time crops in seasonal wetlands. The factors that may have accounted for the stability of grassland and glaciers LULC classes could be due to moisture fluctuations due to annual precipitation amount variations, mostly affecting wetlands, grassland, and bare lands (Hyandye and Martz 2017). However, it is worth noting that the probability values reported in this study are a subject of the magnitude and direction of change in land use, future changes in magnitude and direction of changes may result in a different pattern.

Future land use/cover changes
There is no particular method to follow in comparison to the simulated future map with reality in the future (Jokar Arsanjani and Kainz 2011). This study adopted a technique described by Arsanjani, Kainz, and Mousivand (2011). The location of the predicted cells for different land use/covers was assessed based on the feasibility of occurrence on the predicted sites. Just by visual interpretation of the LULC maps of 1993LULC maps of , 2006LULC maps of , and 2018 (Figure 7), it becomes evident that the builtup areas expanded enormously in the KWK catchment. The observed changes observed between 2018 and 2030 are relatively higher than those observed in the past. This could be attributed to the transition probability between 1993 and 2006. Nevertheless, validation results show that the model correctly captured the driving factors. The results can be improved by the addition of more years during model training and validation. Most of the factors deriving the land-use change in Mt. Kilimanjaro slopes were included in this model. However, the model results are based on the assumption that the current drivers will continue to operate in a similar pattern without significant changes. The changes taking place between the base year (2018) and 2030, e.g. changes in government policies or implementation of the major projects are not included in the model simulation.
Also, the results in Figure 8 show that glacier ice will increase in 2030. However, the interpretation of this increment needs to be taken with caution; this is because drivers for de-glaciation and re-glaciation are more likely linked to the glacier type and changes in climate (Geilinger 1936;Kaser et al. 2010) which except precipitation are not part of this model. It is also worth mentioning that this model includes only a small fraction of the entire Mt. Kilimanjaro glacier.
The results in this study are consistent with the global trend reported by Cao et al. (2019b) showing increases in agricultural land from 2010 to 2100 in all three scenarios used. The report further shows shrinking in the areas under forest, grassland and shrubland. Similar results were reported in Doha; the LULC prediction shows expansion in built-up areas by 20.1% from 2010 to 2020. (Hashem and Balakrishnan 2015). Also, in Long Island Sound Watersheds, USA, the reports show a continued loss of forest in the near future (2026) (Zhai et al. 2018).
Impacts of the LULC change in the hydrology of the area are out of the scope of the current work. However, it is worth mentioning that the influence of land use on terrestrial hydrology is stronger than climate change (Taylor et al. 2013). Land-use changes over Kilimanjaro are likely to affect the water balance of the mountain, which is known to be a water tower of the East African region (Viviroli et al. 2007). The mountain supplies water to South-eastern Kenya and North-eastern Tanzania, where demand for freshwater for domestic, hydropower production, industrial, and irrigation is arguable of rapid increase (Agrawala et al. 2003). The mountain plays a vital role in the regulation of regional climate and provides ecosystem services to communities dwelling on the slopes (Hemp 2005). Also, LULC change leads to losses of the most effective water source in the fog interception zone of Mt. Kilimanjaro (Hemp 2005). Also, rapid increases in deforestation and expansion of agricultural activities resulted in grouping the Kahe catchment aquifer among the degraded mountain catchments (Grove 1993). In other catchments such as Lake Babati catchment, agriculture increased from 23% to 58% of the area between 1960 and 1990 years; accounting for a recharge drop from about half to a quarter of the annual rainfall (Sandström 1995). Therefore, future works need to focus on studying the impacts of LULC changes on the hydrology of Mt. Kilimanjaro slopes.
Apart from being a water tower, Mt. Kilimanjaro slopes landscapes and ecosystems provide a unique avenue for both local and global community in terms of wildlife conversation (Msoffe, Nauss, and Zeuss 2020), water source (Agrawala et al. 2003;Hemp 2005) and tourism (Anderson 2015). Reports show that the vegetation cover in the highlands is rapidly being depleted and the snow cover is being lost, posing a threat to losing Mt. Kilimanjaro's tourists attraction for the global community (Mbonile 2005). The analysis land use and cover change analysis for Mt. Kilimanjaro slopes contributes to a better understanding of the future changes of the landscape and associated impacts which affect both surrounding local community and global at large. Thus, this study provides crucial future-based information for the conservation of Mt. Kilimanjaro slopes environment, benefiting both the local and international community. These ecosystem services are not only important for the local community but also have agency of impacting negatively on the realization of the achievement of local, regional and global sustainable development milestones.
Although climate change is out of the scope of this work, it is worth mentioning that the link between land use and climate changes results in changes in precipitation and temperature patterns. Ultimately affects the realization of the United Nations sustainable development goals, such as food security (Goal2), ending poverty (Goal 1), energy for all (Goal 7), and combating of climate change and its impacts (Goal 13), as well as the loss of biodiversity (Goal 15) (UN 2019). Therefore, the results presented in this paper are crucial for better planning and conservation of the ecosystems and landscapes of Mt. Kilimanjaro slopes so as to continue harnessing the life support benefits to both local and global community.

Implication of land-use changes to staple food production in the region
Cultivated land is a major land use that affects food production and ecological systems. The area of cultivated land increases from 2010 to 2100 in all three scenarios, as shown in Figure 6(a). In contrast, the areas of forest, grassland, and shrubland decrease at a global scale, as shown in Figure 6(b)-6(d). Cultivated land is a major land use that affects food production and ecological systems. The area of cultivated land increases from 2010 to 2100 in all three scenarios, as shown in Figure 6(a). In contrast, the areas of forest, grassland, and shrubland decrease at a global scale, as shown in Figure  6(b)-6(d). Cultivated land is a major land use that affects food production and ecological systems. The area of cultivated land increases from 2010 to 2100 in all three scenarios, as shown in Figure 6(a). In contrast, the areas of forest, grassland, and shrubland decrease at a global scale, as shown in Figure 6(b)-6(d).
Food production and ecological systems are chiefly affected by changes in cultivated land (Cao et al. 2019b). In this study, despite the expansion in agricultural land, it is not reflected in the production of the selected staple food crops (Figure 8). There is no linear relationship between food production and an increase in cropland. Although this increment in LULC is under the assumption of the uniform spatial increment, promising weather, pests, and disease condition throughout the region. However, the observed mismatch shows that the expansion of the production area does not guarantee an increase in the crop yield. Further, the results suggest that other potential factors, such as good agronomic practices, suitable weather, fertile soils, and pest and disease management practices, influence crop yield. Increased farm fragmentation due to applicable land tenure system in the region, soil erosion and land degradation (UNDP 2014) is among the essential factors that can largely explain this mismatch. The observed sharp decrease in banana production may have resulted from the conversion of the banana farms into short seasonal crops such as maize and bean fields due to low prices in the world market which rendered banana unprofitable (Soini 2005).
Although it is worth mentioning low rainfall, especially for drier lowlands, and poor agronomic practices can cause nutrient mining and finally decrease per capital production on the southern slopes of Kilimanjaro. However, the survey shows that the dynamics of land use for agriculture on the slopes of Kilimanjaro are more about intensification on the higher slopes and expansion on lower slopes.
An increase in the population goes together with the need for more food. Thus, the expansion of agriculture in the future will not be a strange phenomenon. Meeting high food demand can be attained through the engagement of the best farming technologies to increase production per unit area instead of expansion of the area under agriculture. Most of the areas under grassland, barren land, and lower area forests are projected to be transformed into the built-up and agricultural lands. Rapid population increases and expansion of agriculture, especially in the lowlands are also reported in some of the previous studies (Mbonile, Misana, and Sokoni 2003;Misana, Majule, andLyaruu 2003, Soini 2005;Chiwa 2012) Increasing pressure on the arable land results in increased land degradation, leading to lower farming yields (CEP 1995). Declining productivity of the highland gardens has shifted the balance such that in recent years, the lowland plots have assumed more importance, providing food for both livestock and people. However, the lowlands are much drier and require modern agronomic techniques. The substantial nutrient mining that goes together with high costs of inputs, fertilizer and seeds has resulted in declining soil productivity in the lowland, thus, intensifying the soil degradation in the area (UNDP 2014). Therefore, restoration and fertilizers are required to improve the fertility and productivity in the lowlands, but, most smallholder farmers cannot afford the associated costs. Thus, making the available land less productive and accessible for crop production.
Land use and environmental destruction are among the factors that would potentially reduce available water for food production in the Pangani basin (De Wit and Stankiewicz 2006). A study by Munishi, Lema, and Ndakidemi (2015) reports a decline in precipitation and uneven distribution of rainfall during the cropping season are essential determinants of the maize bean and Banana yield in the area. Hemp (2005) reported decreases in annual precipitation amounts on Mount Kilimanjaro by 600-1200 mm over the last 120 years while since 1976, temperatures have also increased considerably. These would significantly increase irrigation demand and hence affect crop production in the area (Munishi and Sawere 2014). Small scale farming and mostly rainfed agriculture are a dominant production in Kilimanjaro; this would significantly result in the mismatch between cropping area and crop yield. Mmbaga, Munishi, and Treydte (2017) show that elevation in minimum temperatures as a result of climate change caused local extinction in coffee. Apart from impacting food crops (Munishi, Lema, and Ndakidemi 2015), cash crops such as coffee are also reported to be highly affected (Craparo et al. 2015). It is projected that coffee production will continue declining even up to 2060s if the current land-use changes and climate change trends are not reversed (URT 2012). Thus, a more integrative study focused on the generation of high-quality data on food production, food demand, climate, pests, and diseases are necessary to bring issues close to reality. This study is in agreement with a study by Lyle (1996) that the main challenge lies in production practices. In order to create a balance between environmental land use and food production, there is a need to change our means of production. The focus may be directed towards creating more technologies that increase production in small areas than relying on the expansion of farm size. However, sustainable land use represents a need to take into account environmental, social and economic issues, as the key issues of sustainability (Loures 2019).
Global food production is expected to be impacted by land use and climate changes; hence, threatening food security (Lobell and Field 2007). Achieving food security will require research translation programme focused on transferring research-based evidence of locally adapted sustainable enhancement of food production practices into community-led actions to reverse land degradation and enhance and diversify productivity simultaneously. In the face of the global and local changes in LULC, population increase and decline in food production, this study emphasizes a need for a greater agency in integrating this information into the decision-making process. This decision-making process is relevant to local support not only sustainable food production but also the restoration of degraded land. The decision-making process intervention can offset both local and global changes and practices that can jeopardize food production and environmental sustainability. Further, a thorough understanding of community needs, priorities, policies and practices can create a balance between production and environmental sustainability. Thus, this should be at the heart of the local adaptation measures.

Conclusion
This work links model output with food production as a way of generating information that can help managers to find a balance between sustainable food production and conservation. The silence relationship between an increase in agricultural land and food production is revealed. Thus, the study opens an avenue for other biotic and abiotic factors affecting the production of the selected crops is a gap that needs to be addressed. Future research may be focused to explore other factors accounting for the decline in crop production. Such factors are; agronomic practices, pest and diseases, soil fertility as well as technologies focused on intensification of without further expansion in the farm area.
KWK catchment is one of the areas located on the productive plains at the southern slopes of Mt. Kilimanjaro. It, therefore, experiences extensive alteration in the agricultural and urban cover during the last few decades. Similarly, agricultural land and urban growth might continue to expand in the future, which will impact on land and other resources if the situation goes unattended. The CA-Markov model predicts further extension in agricultural land and built-up areas by 39.52% and 32.55%, respectively, by 2030 from the using 2018 as a baseline year. It is expected that the need for food will increase with an increasing population. However, the expansion in cropland is not reflected in food production data. This means further studies need to explore more influencing factors and modern technologies to increase food production per unit area are necessary to assist in reducing future expansion while meeting the desired target. It is worth mentioning that the dynamics of agriculture on the higher slopes of Mt. Kilimanjaro is more about intensification than expansion. Although a large part of the forest falls in the upper part which is under the conservation area, the marginal 6.37% shrink in its area offers an alarming message to the managers of the conservation area to take necessary management strategies.
Although glacier ice and water were included in this model, it is worth mentioning that there was a serious challenge in the preparation of suitability maps and future simulation accuracy, especially with glacier ice. Factors such as annual rainfall were included in this simulation, but, more factors such as solar radiation and humidity could improve the future simulation accuracy, however altitudinal good quality data are scant or absent. Therefore, future studies could focus on the application of some modern techniques to simulate future dynamics with high accuracy.
Impacts of LULC changes on water resources in terms of the recharge as a response to the land-use change were not considered in the current work. Future works looking at the changes in river runoff and groundwater recharge could provide insights into the impacts of future climate and land use/ cover changes. This could be realized through the identification of recharge areas, the contribution of the fog interception zone to the hydrology of the area and substantial protection laws imposed for future water sustainability. The results of this research can then be used as input to future research focused on determining the impact of land-use changes in hydrological processes and conditions against flooding in the KWK catchment. CA-Markov model is based on the initial land cover distribution and transition matrix. It assumes that the drivers, which have created the current situation for the region landcover, will continue to operate as before in the future (Guan et al. 2011). Also, socio-economic drivers of LULC are assumed to remain as per the current situation in future. Thus, simulating different LULC change scenarios with different changes in driving forces can provide alternative future trajectories. However, the variables used were measured over prolonged time, it is believed that it accounted for dynamic regular driving forces that were happening in time series, and during the prediction, except extreme events such as prolonged droughts or civil wars.