2.1 The Study Area: Suluh River Basin
The study area is found in the highlands of the north-eastern part of Tigray region, northern Ethiopia. The geographic location of the river basin extends from 39°24'59.06 ” E to 39°26'22.73 ” E latitude and 13°38'18.27 ” N to 14°13'53.29 ” N longitude (Fig. 1).The total area of the basin is about 930 km2.The elevation of the study area varies from 1700 to 3,298 meters above sea level. The study watershed falls in four districts (Saesie Tsaeda Emba ,Hawuzen, Kilte Awlealo and Degua Tembien ) of Eastern and South East.
The climate of SRB is characterized as semi-arid. The warmest months are May & June &coldest months are December and January. Average annual rainfall and temperature from 1988–2018 is 420.4 mm and 17.5°C respectively (Fig. 2).The rain fall season is between June to early September which is mono- modal rain fall distribution.
Regarding the hydrological condition, drainage pattern of the SRB is dendrite (Zenebe et al., 2019).The longitudinal profile of the SRB main stream has a length of 94.77 km starting from volcanic mountains of Mugulat & ending at the embouchure into Genfel. The geology of the basin is trap basalt accounts 2.8%, granite and shale accounts 1.8%, metamorphic rock accounts 28.9%, limestone accounts 13.9%, sandstone accounts 52.7% (Zenebe et al.,2019 ).The major soil types of SRB area are haplic lixisols covers 41.4%,lithic leptosols covers 22.7%, Eutric leptosols covers 17.8%,Chromic Cambisols covers 15.6% &Vertic Cambisols covers 2.5%. The major soil textures of the study area are sandy clay loam accounts 41.4%, sandy loam accounts 40.5% and clay accounts 18.2% (Zenebe et al.,2019).
Population Census Commission (2015) shows that in 2015 the population density of the river basin was 142 persons per km2. The dominant economic activity in the basin is agriculture. The crops grown vary according to altitude. The main crops in the highlands are barley, wheat, maize, teff & pulses. Sorghum is widely grown in the lowlands. Cultivation is done by means of the traditional ox-drawn plough. The main livestock are cattle, sheep, goats, donkeys and chickens. There is a livestock feed crisis, resulting from crop residue and vegetation biomass reduction (Alemayehu et al., 2009 & Aredehey et al., 2018). Most of the areas are highly cultivated even on steep valley sides and the remnant original highland forests are being overgrazed & deforested, resulting in shallow soils due to serious soil erosion by water. The major strategies for land management in Tigray highlands are building of stone terraces, micro dams including gully treatment, establishment and development of enclosures and community woodlots, enforcement of used rules, regulations for grazing lands & reduced burning activities(Sembroni et al.,2017 ; Zenebe et al.,2019 & Hishe et al.,2020).
2.2. Sources of Data
Sources of data used, data processing image segmentation, nearest neighbor fuzzy classification, accuracy assessment & LULC change modeling were done according to the overall workflow of Fig. 3. Free satellite images (Landsat-5 TM of 1990, Landsat-7 ETM+ of 2002 & Landsat-8 OLI–TIRS of 2018) were used for the LULC change analysis of the study area (Table 1).These datasets were acquired from the National Aeronautics and Space Administration (NASA) through their EOS Data Gateway Database. Land sat images that were used for this study area freely available. The Landsat-7 ETM+ 2002 & Landsat-8 OLI–TIRS 2018 of a 30 meters pixel were resampled to a 15 meters pixel size. Digital Elevation Model (30 m) based on Aster imagery and ancillary data (topographic maps, field thematic layers (roads and towns), & village and district boundaries) were also utilized during analysis. All data were projected using Universal Transverse Mercator projection system zone 370 N and datum of World Geodetic System 84 (WGS84). An intensive pre-processing such as layer-stacking, resolution merge, and sub setting were carried out and to remove disturbances such as haze, noise, steep slope effect & radio metric variation between acquisition dates.
Table 1
The characteristics of land sat satellite data
Sensor | Path / row | Acquisition time | | Spatial Resolution Resolution | Sensor |
Landsat TM | 169/50 | 01/07/1990 | | 30 m | TM |
Landsat ETM+ | 169/50 | 03/08/2002 | | 15 m | ETM+ |
Landsat OLI-TIRS | 169/50 | 03/12/2018 | | 15 m | OLI-TIRS |
Aster DEM | | | | 30 m | |
Ancillary data | | | | | |
Topo map | | | | 1:50000 | |
Field data | | | | Nov.2017-Jan 2018 | |
Roads and towns | | | | | |
District boundary | | | | | |
Village boundary | | | | | |
We applied nearest neighbor fuzzy (Eq. 1) classification in eCognition Developer 9.2. It is used in eCognition-automatically by generating multidimensional membership functions (Foody,1999; Groenemans et al.,1997;Zhang & Foody,1998;Yan et al., 2006;Salman et al.,2008, Zhou et al.,2008 and Kalantar et al.,2017).We designed eight LULC classes (cultivated, bar, built up, grazing, plantation, shrub and bush, water body & forest land) based on the Table 2.
Table 2
Land cover, Land use types and their descriptions
LU/LC classes | Descriptions |
CL | Areas covered with annual crops followed by harvest and bare soil |
FL | Areas covered with a natural forest community with a closed, deep and complex canopy often consisting of several crown layers |
GL | Landscapes that have a ground story in which grasses are the dominant vegetation forms |
SBL | This category includes low woody plants, generally less than three meters in height, usually with multiple stems, growing vertically |
BL | Areas with little or no vegetation cover consisting of exposed soil and/or bedrock |
BUL | Areas for construction sites and towns |
PL | Areas composed of transplanted seedlings of Cactus, Eucalyptus globules and Cupresus spp. |
WB | Includes rivers, reserves, lakes (artificial dams and natural lakes) and so on |
Note: Forest land (FL), cultivated land (CL), shrub-bush land (SBL), built up land (BU), grazing land (GL), bare land (BL), plantation land (PL) and water body (WB) |
A={(X, µA(x));x ϵ X},Where µA →[0,1] (1)
Where A = fuzzy set X = a space of objects X = elements belonging to space X µ – membership
Function
In this study, over all accuracy (Eq. 2) and Kappa coefficient (Eq. 3) were used to assess the accuracy of the classified images. As a result, we found an excellent accuracy in which as to the kappa coefficient results reveals 0.886, 0.883 & 0.852 for the years of 1990, 2000 and 2018, respectively (Table 3).Since the values falls above the cut point of the standard overall classification accuracy level of 85% (Fleiss et al.,2003;Doxani et al.,2008 and Congalton & Green ,2019) with no class less than 70% (Kumar et al.,2018).
Table 3
Summary of error matrixes for the classified images of 1990 and 2002.
LULC classes | 1990 | 2002 |
User Accuracy | Producer Accuracy | User Accuracy | Producer Accuracy |
BL | 56% | 100% | 65% | 100% |
BUL | 88% | 100% | 89% | 100% |
CL | 100% | 55% | 100% | 56% |
FL | 100% | 94% | 100% | 100% |
GL | 100% | 94% | 100% | 100% |
PL | 86% | 100% | 88% | 100% |
SBL | 78% | 100% | 79% | 100% |
WB | 94% | 100% | 100% | 100% |
Overall Accuracy | 87.12121212 | 89.79592 |
Kappa Accuracy | 0.852591473 | 0.883327 |
$$\text{OA}=\frac{X}{Y}*100 \left(2\right)$$
Where, OA is overall accuracy, x is number of correct values in diagonals of the matrix, and y is total number of values taken as a reference point.
K=\(N\sum _{i=1}^{r}xii-\sum _{i=1}^{r}(xi\times x+1)\)/\(N2- \sum _{i=1}^{r}xii-\sum _{i=1}^{r}(xi*x+1) \left(3\right)\)
Where: r = is the number of rows in the matrix
Xii = is the number of observations in rows i and column I (along the major diagonal) Xi + = the marginal total of row i (right of the matrix) Xi + 1 are the marginal totals of column i (bottom of the matrix) N is the total number of observations. K = kappa coefficient
LULC Modeling using LCM
LCM model is based on the artificial neural network (ANN), Markov Chain matrices & transition suitability maps, generated by training multilayer perceptron (MLP) or logistic regression (Ramachandra et al.,2013., Mishra et al.,2014, Roy et al.,2014, Gibson et al.,2018, Dzieszko,2014). Megahed et al.,2015, Ansari & Golabi,2019). This model predicts the LULC changes from the thematic raster images having the same number of classes in the same sequential order (Mas et al.,2012). In this study, the LCM is used to predict the future LULC changes SRB for the next 30 years by following four steps, namely: (1) change analysis, (2) transition potential and determination of explanatory variables, (3) change prediction(4) model validation (Megahed et al.,2015).
I.Change Analysis
In the change analysis panel, the changes between two different time periods time 1(1990) & time2 (2002) LULC maps were calculated. To calculate the annual rate of change Peng et al.(2008) formula(Eq. 4) was adopted:
C=\(\frac{\varDelta f-\varDelta i}{\varDelta i}\times \frac{1}{T}\times 100 \left(4\right)\)
Where C = is the annual change rate of a given LULC type, \(\varDelta\)f and\(\varDelta i\) are the final and initial area coverage of LULC type during the specific period and T = year difference between the initial and final period.
The gain and loss(Eq. 5) was also calculated adopted from LCM in IDIRISI software.
\(\left[\text{P}\text{l}\text{o}\text{s}\text{s}\right(\text{i}),\text{j}\) = (Pj,i- pi,j)/( pi - pi)*100 i#j]
\(\left[\text{P}\text{g}\text{a}\text{i}\text{n}\right(\text{i}),\text{j}\) = (Pi,j- pj,i)/(pi- pi)*100 i#j]\(\left(5\right)\)
Where \(\text{P}\text{l}\text{o}\text{s}\text{s}\left(\text{i}\right),\text{j}\) Is the percentage taken by j in LUCC in total ‘conversion loss’ of category row i\(\text{P}\text{g}\text{a}\text{i}\text{n}\left(\text{i}\right),\text{j}\) Is the percentage taken by j in LULC change in total ‘conversion gain’ of category row i, pi,j and pj,i
II. Transition Potential Modeling and Driving Forces Determination
Transition Potential
The transition potential determines the area of change (Megahed et al., 2015).LULC changes with common driving variables were grouped into sub-models. In addition; evidence likelihood was selected to determine the relative frequency of different LULC types which had occurred within the transitional areas (Megahed et al.,2015).
Selection of Explanatory Variables
Explanatory variables responsible for LULC change were selected (Gibbs et al.,2010 & Li et al., 2015).In the current study, biophysical (rainfall, slope and elevation), socio-economic variables (distance to rivers, distance to roads, distance towns) and demographic variables (population density) were identified as variables. For testing, selection and transition of Model variables, Cramer’s V coefficient was used.All the variables were then add into transition sub-mmodel. For transition sub-model structure and running the model MLP neural network were used. In this study, the transition pixels to be occurred in between 1990 to 2002 were assign randomly to one of two groups: the training set & the testing set. The variables were derived via geographical and geo-statistical elaborations of a geographic information system & were formalized as follows (Eq. 6):
X = x1, x2, …, xn\(\left(6\right)\)
Every variable associated with a neuron in the input layer was normalized using (Eq. 7):
Xi=( Xi-min)/(max-min)\(\left(7\right)\)
In the hidden layer, the signal that is received by neuron j in the hidden layer for pixel k was calculated as follows(Eq. 8):
netj(k)=\(\sum _{i}\text{W}\text{i},\text{j} \text{X}\text{i}\left(\text{k}\right) \left(8\right)\)
where netj(k, t) is the signal that is received by neuron j, and wi,j is the weight between the input layer I and the hidden layer j. The output layer has 2 neurons that correspond to 2 possible significant states (1 = transition, 2 = permanence of the pixel); the neuron l generates a value that indicates the transition probability. Transition probabilities can be calculated using a sigmoidal function using Eq. 9 (a sigmodial function is used to represent the non-linearity of each node):
p(k,1=1)= \(\sum _{i}\text{W}\text{j},\text{i}\frac{1}{1+{x}^{netj(k,t)}}\) \(\left(9\right)\)
III. Change Prediction
Change prediction is the last step in which the future prediction is executed on the basis of MC, and using the historical rate of change and the transition potential maps (Wang et al.,2012). MC analysis were run for this study to determine the amount of change using two LULC maps (1990 and 2002) along with the date specified (2018,2028 and 2048). The steps determines how much land was expected to transition from the later date (2002) to the prediction date (2028 and 2048). A MC(Eq. 10) is made of s, a vector of the distribution of LULC classes at time t, and A(ፐ), a matrix of transition probabilities from land use u to land use u’ in a given time interval(ፐ)
St+ፐ=A(ፐ)St\(\left(10\right)\)
Using hard prediction model, a map for 2018 of study area was simulated in order to compare with the ‘actual’ land cover map of 2018.
IV. Future Scenario
Hard predictions land change modelers was applied in this study.
V. Model Validation
It is needed to assess the accuracy. In this study validation was done by comparing the simulated and the actual LULC maps of 2018. The qualitative information collect using direct observation, focus group discussion and interview were analyzed and interpreted using qualitative techniques, whereas quantitative data collect were analyzed using descriptive statistics ( like percentage).