Hydrological Response of Dry Afromontane Forest to Changes in Land Use and Land Cover in Northern Ethiopia

Abstract

This study analyzes the impact of land use/land cover (LULC) changes on the hydrology of the dry Afromontane forest landscape in northern Ethiopia. Landsat satellite images of thematic mapper (TM) (1986), TM (2001), and Operational Land Imager (OLI) (2018) were employed to assess LULC. All of the images were classified while using the maximum likelihood image classification technique, and the changes were assessed by post-classification comparison. Seven LULC classes were defined with an overall accuracy 83–90% and a Kappa coefficient of 0.82–0.92. The classification result for 1986 revealed dominance of shrublands (48.5%), followed by cultivated land (42%). Between 1986 and 2018, cultivated land became the dominant (39.6%) LULC type, accompanied by a decrease in shrubland to 32.2%, as well as increases in forestland (from 4.8% to 21.4%) and bare land (from 0% to 0.96%). The soil conservation systems curve number model (SCS-CN) was consequently employed to simulate forest hydrological response to climatic variations and land-cover changes during three selected years. The observed changes in direct surface runoff, the runoff coefficient, and storage capacity of the soil were partially linked to the changes in LULC that were associated with expanding bare land and built-up areas. This change in land use aggravates the runoff potential of the study area by 31.6 mm per year on average. Runoff coefficients ranged from 25.3% to 47.2% with varied storm rainfall intensities of 26.1–45.4 mm/ha. The temporal variability of climate change and potential evapotranspiration increased by 1% during 1981–2018. The observed rainfall and modelled runoff showed a strong positive correlation (R² = 0.78; p < 0.001). Regression analysis between runoff and rainfall intensity indicates their high and significant correlation (R² = 0.89; p < 0.0001). Changes were also common along the slope gradient and agro-ecological zones at varying proportions. The observed changes in land degradation and surface runoff are highly linked to the change in LULC. Further study is suggested on climate scenario-based modeling of hydrological processes that are related to land use changes to understand the hydrological variability of the dry Afromontane forest ecosystems.

1. Introduction

Forests provide multitudes of essential ecosystem services [1], such as soil nutrient cycling, build-up of organic matter, and water retention. Forest biodiversity losses may seriously jeopardize the functioning of forest ecosystems, and consequently the ability of forests to provide ecosystem services [2]. In the past few decades, tropical landscapes were facing critical environmental challenges [3,4] that were related to rapid land use and land cover (LULC) change, such as massive deforestation and land conversion. The Forest Resource Assessment (FRA) that was completed by Food and Agriculture organization (FAO) in 2005 and preliminary results of FRA 2015 [5] indicate that, on a global scale, the total forest area continues to decrease by 3.1%, which is about 0.6% per year from 2000 to 2005 [3,6] and 3.3% per year from 2010 to 2015 [5]. Deforestation continues at the highest rates across the tropics, and FRA estimates for 1990–2005 suggest annual rates of deforestation at 0.9% in South and Southeast Asia, 1.2% in Central America, 0.45% in South America, and 0.62% in Africa [3]. This means that, on average, each year from 1990 to 2015, an area of 4128–3999 million ha, which is above 3.1% of the tropical forest, was lost [3]. The change in land use is considered to contribute to 20% of global greenhouse gas emissions [7,8]. In particular, forest degradation in sub-Saharan Africa is the main source of land-based emissions [9,10].

According to Pielke [11], LULC changes influence climatic and weather conditions, ranging from local to global levels. Biophysical factors, such as topography, soil type, and seasonal rainfall, as well as demographic, political, economic, and social factors combined together lead to large changes in land use in East Africa, and these trends are expected to continue into the future [3,5]. Similarly to the above studies, Ethiopia is experiencing significant LULC dynamics, such as land conversion from natural vegetation to agricultural and urban areas [12,13,14,15]. The problem is more severe in the highlands, which account for 44% of the country’s landmass and they have been cultivated for centuries [16,17,18,19,20,21,22,23,24]. Land use dynamics during 1975–2016 in India reported that agricultural land increased by 18%, whereas forest land and rangeland decreased by 7% and 10%, respectively [25]. Specifically, in the northern highlands of Ethiopia, increasing land degradation due to the loss of the remnant forests, which are found in sacred and inaccessible places [26] and cultivation of the sloppy lands [19] aggravates runoff and reduces the infiltration capacity of the area. The change in land use is considered to contribute to 20% of global greenhouse gases emissions [7,8]. Land cover change is defined as the change in land cover and its associated properties over time [27,28]. LULC change constitutes one of the fundamental environmental problems in Ethiopia [29,30,31,32]. For example, a study in the Afar region of northern Ethiopia reported that the invasion of Prosopis juliflora increased at annual rates of 31,127 ha, while grassland and bush-shrub-woodland declined at 19,312 ha and 10,543 ha, respectively [31].

The findings of Palma et al. [33] also indicate that land use changes have the potential to affect landscape hydrology, thereby aggravating runoff. The hydrological response is a combined indicator of the watershed condition and significant changes in LULC, reflecting the overall health and function of a semi-arid watershed [34]. Other studies [18,19,35] showed that the simulated impact of climate change on hydrological processes significantly varies, depending on the climate model and emission scenarios, in addition to geographical variation [36,37].

The alteration of tropical landscapes, primarily the conversion of natural forests to agriculture or pasture [11], directly influence the hydrology of the landscape, and the relative impacts of different potential drivers of these extreme events on the hydrological cycle are still debated. External factors include meteorological variations, which result in periods of increased extreme events or climatic non-stationarity. The internal factors include abstractions and discharges in the river, land use change, and modification of the river flow channel.

LULC changes significantly affected regional hydrology [38,39,40,41], and consequently the delivery of ecosystem services [42,43] and climate processes. The impact of LULC on environmental processes is not universal and it depends on the local context of a particular region [44]. Gebremicael et al. [45] revealed that more than 72% of the northern Ethiopian landscape was significantly altered during 1972–2014 in the head water of the Tekeze-Atbara Basin. A study by Guzha et al. [46] in East Africa has shown that forest cover loss increased the annual discharges by 16% and surface runoff by 45%. Bewket and Sterk [47] showed that changes in land management practices are the main factors of variability in water availability in Ethiopia. In particular, Gebresamuel et al. [13] have reported that the changes in LULC decreased the water-storage capacity of soils by 1.63 mm/y at Gum Selassa and 1.09 mm/y at Maileba of northern Ethiopia, with a corresponding increase in the surface runoff by 2.7 mm/y and 2.3 mm/y, respectively. The other concern is the effect of global climate change, mainly relating to the patterns and amounts of rainfall in the upper catchment areas, which directly affect runoff patterns [48]. Significant variation in evapotranspiration occurs due to LULC accompanied with the leaf area index change [49]. Generally, land use change can lead to a significant change in groundwater recharge and base flow [50], flood frequency and interval [51], peak runoff [52], and total suspended sediment and nutrient concentration [53]. Land restoration has been successful in improving the vegetation cover, and consequently the hydrological processes of degraded landscapes [47]. For example, the restoration of steep hillsides enhances the initial retention capacity of the soil, thereby limiting the runoff curve number [13,35,54,55,56].

Physically based distributed hydrological models have become important for the understanding of fundamental physical processes that underlie the hydrological cycle and developing scenario-based solutions to adverse changes in forest hydrology at different spatial scales [57]. The development of satellite remote sensing is important for modelling capabilities to evaluate and predict the hydrological consequences of land-use change at multiple scales. For example, the conceptual and empirical foundations of the soil conservation service curve number (SCS-CN) method have been applied to a wide range of catchments across the world to assess the effects of land cover change on surface runoff [27,58,59,60,61,62,63,64]. The runoff coefficient forecasts direct surface runoff volume for a given rainfall event and estimates the volumes and peak rates of surface runoff in catchments [58]. The rainfall-runoff relations within a watershed are primarily driven by the relationship of factors, such as climate, land cover, and soil [36,37,38,39]. The SCS-CN model expresses runoff volume as a function of rainfall volume, hydrologic storage, and initial abstraction [65]. The CN value depends on the land surface features and hydrological soil conditions of the entire watershed.

The original SCS-CN method computes direct runoff by only considering the available rainfall on the current day without taking the effect of the moisture available prior to the storm into consideration. In contrast, the curve numbers are sensitive to antecedent conditions [66]. The existing SCS-CN method does not contain any expression for time, and it ignores the impact of rainfall intensity and its temporal distribution. There is no explicit provision for spatial scale effects either. The other demerits of the existing SCS-CN method are the absence of clear guidance on how to vary the antecedent moisture conditions and the fixing of initial abstraction ratio 0.2, preempting a regionalization that is based on geological and climatic settings [67].

Research on runoff patterns in response to land cover dynamics allows for the assessment of sustainability of land use systems in dry areas, as the runoff reflects the ecological state of the entire watershed. Modeling allows evaluating long-term hydrological consequences of LULC change in dry Afromontane forests and provides quantitative information for land use planning and water resources management. However, developing countries face difficulties in different aspects of hydrological modeling research, such as the availability of continuously measured data. Therefore, combining the remote sensing technique, point-based meteorology data, and hydrological modeling can help in understanding and evaluating the existing LULC conditions of the undulating and mountainous dry Afromontane forests. This information can also be employed to forecast the likely effects of any potential changes in land cover on water resource systems. Hence, such a study has practical relevance for devising strategies and policies for sustainable land and water use. No studies have been published to date combining the remote sensing technique and evapotranspiration and modeled surface hydrology information in forest areas of the dry Afromontane forest in northern Ethiopia, to the best of our knowledge.

The dry Afromontane forest landscape in northern Ethiopia is characterized by prolonged dry seasons, under conditions of low humidity, high air temperatures, and drying up of seasonal rivers and streams due to natural and human-induced factors [55,68]. Some studies on streamflow patterns with respect to land cover dynamics enabled the assessment of sustainability of land use systems, because stream flows are reflections of the ecological state of the entire watershed [47,55,69,70]. However, little is known regarding land use dynamics with respect to catchment hydrology [47,69], soil degradation [71], soil properties, and surface runoff [13] in Ethiopia. This study integrated high-resolution temporal LULC maps that were derived from satellite imagery with hydrological modeling to evaluate long-term hydrological consequences of LULC change in dry Afromontane forests and provide quantitative information for land use planning and water resource management. In particular, we aimed to adapt the conventional SCS-CN model by incorporating the variation of daily curve numbers with respect to the variability of antecedent rainfall, potential evapotranspiration (PET) with integrated LULC change images, and grid-based hydrological soil texture classification groups (HSG). The modified SCS-CN-based models provide a hydrologically sound procedure for a more accurate representation of the catchment behavior through analysis of the hydrological response of LULC changes, as indicated by their impact on the PET and runoff coefficient in the dry Afromontane forest landscape of northern Ethiopia.

2. Materials and Methods

2.1. Study Area Description

Topographically, northern Ethiopia is characterized by an undulating to steep terrain that is frequently divided by stream incisions. Thus the northern Ethiopian climate exhibits extreme climatic variations within short distances, from very dry tropical to sub-humid and subtropical to highly-moist tropical climates. The average annual rainfall varies from 200 mm in the arid lowlands to over 2200 mm in parts of the southwestern highlands. The mean annual temperatures vary from above 35 °C in the lowlands to less than 15 °C in the highland areas.

The estimated area of Ethiopian drylands is over 75 million ha, which accounts for about 66–72% of the total area of the country [72,73,74]. The study area lies between 39°10′E–40°02′34.08″ E and 12°53′29.76″N–12°′2.88″N (Figure 1). The study was carried out in two neighboring fragmented state forests of Hugumburda Grat Kahisu and the Tabotat natural forest areas of northern Ethiopia. This priority state forest has a minimum mean annual temperature range between 6.3 °C and 20.6 °C. The maximum mean annual temperature ranges from 15.1 °C to 29.7 °C. The mean annual rainfall ranges between 350 mm to more than 1000 mm (

Table 1. Mean climate data from 1981–2018 years of the study area.

Month	Average Temperature (°C)	Min. Temperature (°C)	Max. Temperature (°C)	Precipitation (mm)
January	15.5	6.8	24.2	22.98
February	15.7	7.4	24.1	13.6
March	17.1	8.8	25.4	60.2
April	17.8	10.3	25.4	81
May	19.4	10.3	28.5	61.13
June	20.6	11.5	29.7	46.6
July	19.3	13.3	25.3	231.4
August	18.5	12.3	24.7	280.4
September	18	10.6	25.4	82.12
October	16.3	7.5	25.1	47.7
November	15.9	6.4	25.4	20.6
December	15.1	6.3	24	44.3

) in the lowland of Raya Azebo and the highlands of Korem districts, respectively. Mixed farming systems dominate the local agricultural activity. The dominant food crops grown are Zea mays, Sorghum bicolor, Triticum durum, Eragrostis tef, Hordeum vulgare, Pisum sativum, and Cicer arietinum. High-value tree crops, such as Mangifera indica, Persea americana, Carica papaya, and Psidium guajava, in the lowland and Malus domestica in the highlands of the Endamohoni and Offla districts are produced in large quantities. Cattle, sheep, goats, equines (donkey, horse, mule), and camel are the major livestock species reared.

The indigenous dry Afromontane forests of Northern Ethiopia are dominantly characterized by Juniperu sprocera, Olea europaea subsp. cuspidata, Maytenus arbutifolia, Maytenus senegalensis, Dodonaea angustifolia, Hagenia abyssinica, Bersama abyssinica, Acacia abyssinica, and Podocarpus falcatus in the highlands and midlands, while various Acacia species dominate the lowlands.

2.2. Methodology

2.2.1. Remote Sensing Data Collection and Processing

Satellite images for the years 1986, 2001, and 2018 were acquired to analyze the spatio-temporal patterns of the LULC of the fragmented state forests of Hugumburda Grat Kahisu and Tabotat natural forest areas. Landsat 5-Thematic mapper (TM) and Landsat-8 Operational Land Imager (OLI) images for the selected years were obtained from the US Geological Survey (USGS) Centre for Earth Resources Observation and Science (EROS), found in http://glovis.usgs.gov/. A brief description of the collected satellite images, including satellite and sensors used, resolution, acquisition date, 168/051 path/row of each images for each period is summarized in

Table 2. Landsat images used for the analysis of land use and land cover (LULC) (Source: http://glovis.usgs.gov).

Sensor	Month/Day/Year	Resolution	Path /Row
Landsat 5 TM	12/15/1986	30 m	168/51
Landsat 5 TM	01/28/2001	30 m	168/51
Landsat 8 OLI	01/13/2018	30 m	168/51

. Data acquisition was entirely conducted during the dry season to avoid phenological effects and ensure cloud-free images.

Remotely sensed data were pre-processed while using ERDAS Imagine 2015 software. Image rectification, restoration, enhancement, classification, and accuracy assessment were conducted using this software. ArcGIS 10.6 and ENVI 5.3 software were employed for managing, analyzing, combining, and mapping spatial data.

The analysis was performed on freely available satellite imagery products from Landsat 5 and Landsat 8 OLI at 30 m resolution. Spatial and temporal data from remote sensing integrated to World Geographic Coordinate System (GCS WGS 1984) was used to illustrate the Ethiopian administrative boundaries. The study aimed to determine the hydrological response of the fragmented state forests of northern Ethiopia to the dynamic LULC. Satellite images for a dry month have been considered for analysis to avoid seasonal factors, such as phenological effects. LULC information was then extracted while using the supervised maximum likelihood classification method by collecting 350 ground control points (GCP) for training signature generation and another 350 GCPs for accuracy assessment. Soil data were derived from the initial soil map, which are available online from the International Soil Reference and Information Centre (ISRIC) SoilsGrid250 [75] and modified based on the study area characteristics.

2.2.2. Data Acquisition and Processing

Stow [76] stated that accurate per-pixel registration of multi-temporal remote sensing data is essential for change detection analysis, since possible registration errors might be interpreted as LULC changes, leading to an overestimation of the actual change. In this study, image pre-processing (geometric and atmospheric corrections and topographic and temporal normalizations) was performed for all LULC maps. All of the satellite images and digitized ancillary paper maps were georeferenced to a common coordinate system while using the topographic map at 1: 50,000 scale, the Universal Transverse Mercator (UTM) map projection (Zone 37), and WGS 84 datum in ArcGIS 10.6 and Erdas Imagine 2015.

The digital elevation model (DEM) of the Hugumburda Grat Kahisu state forest was downloaded from the SRTM (http://gdem.erssdac.jspacesystems.or.jp/) with a resolution of 30 m. The downloaded tiles were merged using mosaic capabilities of Arc GIS 10.6 to form a single DEM of the study area since 30 by 30 DEM resolution of SRTM is usually stored as tiled datasets. DEM was used to simulate the stream network and delineate the watershed into a series of interconnected sub-basins. The soil map, available online from ISRIC SoilsGrid250 [75], was modified based on the study area characteristics that were provided by the National Meteorological Agency (NMA). The climatic data, such as daily precipitation, as well as minimum and maximum temperature, were collected from Korem meteorology stations located within the state forest.

Climate data were provided by the National Meteorological Agency (NMA). The temporal climatic data used for this study were daily precipitation and daily minimum and maximum temperature collected at Korem and Hashengie meteorology stations located within the state forests.

2.2.3. Land Cover Dynamics

Supervised classification was used to classify the LULC change of the areas in the northern Ethiopian mountainous ecotones. A total of seven LULC classes, namely shrubland, grassland, bare land, forest, cultivated land, water bodies, and built-up areas were specified. All seven classes were identified in all the images for the selected study years in a consistent manner. For training points, more than 250 sample plots per LULC class were randomly assigned. Simple random sampling was employed to generate reference data while using high resolution Google Earth images and expert knowledge.

Accuracy assessment was conducted on the resulting classified imagery. This process includes generating a set of points in the classified imagery and comparing them with actual points on the ground through field work using GPS point-based sampling. GPS points for each land cover class were collected at the field level to complete the accuracy assessment. The LULC classification assigned to each pixel was then compared with the same location on the reference sources to determine whether the classification result was accurate. Field visits were conducted from December 2018 to February 2019 for collecting data on the existing land use type. A transect line was laid across the gradient with a distance of 2 $\times$ 2 km² between the plot and randomly among the transect lines. More than 750 sample plots were gathered while using GPS with a plot size of 30 $\times$ 30 m² by assuming the pixel size of the Landsat data for the undulating mountainous areas of HGK state forest in northern Ethiopia. These ground truthing data were used to train the maximum likelihood classifier (MLC) classifier. This gradient classification method allocates all of the pixels in a dataset to clusters that are defined by the mean value of the training datasets. The actual classification was carried out after the training data had been established (half of the sample data) and the MLC algorithm that was selected for the classification. MLC uses a parametric statistical approach to prepare the probability density distribution functions for each individual class [11,77]. Hence, MLC considers not only the cluster center, but also its shape, size, and orientation [78,79]. The assumption of most MLCs is that the statistics of the clusters have a ‘normal’ (Gaussian) distribution. The data that were collected during ground truthing were stored in Excel and converted to shape files using Arc GIS 10.6 software, after which they were used for LULC classification, analysis, and accuracy assessments in ERDAS 2015 version software.

Figure 2 illustrates the overall conceptual frame work we followed to assess and estimate the annual runoff using SCS-CN curve number model.

2.3. Estimation for Runoff Response

Long-term hydro-meteorological data is required to understand the hydrological response of a catchment. However, hydrological data is often lacking in the developing world, and particularly in the study area in Hugumburda Grat Kahisu and Tabotat fragmented natural forests of northern Ethiopia. In the absence of sufficient hydrological data, empirical models, such as SCS-CN (USDA, 2004b, 2004a), are used as a substitute [13,54,80]. The SCS-CN method is based on the water balance equation and two fundamental hypotheses. The first hypothesis equates the ratio of actual amount of direct surface runoff Q to the total rainfall P (or maximum potential surface runoff) to the ratio of actual infiltration F to the amount of the potential maximum retention S. The second hypothesis relates the initial abstraction Ia to the potential maximum retention S. This study applied the same method to estimate the direct runoff response. The model quantifies the effect of changes in rainfall and land cover on the hydrological response of the dry Afromontane forest-based catchments.

Although direct runoff estimation is usually carried out at two abstraction ratios i.e., 0.05 and 0.2, studies on the highlands of northern Ethiopia [54,80,81] recommended λ = 0.05 as an optimum initial abstraction ratio that is based on least squares fitting for most experimental plots. The model’s statistical algorithm can be written as:

$$Q_d=\frac{{\left(P-I_{a}S\right)}^2}{\left(P+\left(1-I_a\right)S\right)},\mathrm{when}\mathrm{p}>{\mathrm{I}}_{\mathrm{a}}\mathrm{S}$$

(1)

$$Q_d=0\mathrm{when}\mathrm{P}<{\mathrm{I}}_{\mathrm{a}}\mathrm{S}$$

(2)

$$I_a=\lambda S$$

(3)

Where Q_d is the estimated runoff (mm), P is the measured daily rainfall (mm), I_a is the initial abstraction (mm), λ is the initial abstraction ratio, and S is the maximum water retention parameter (mm) determined from the weighted CN value. S is related to the dimensionless runoff curve number (CN), which is given by the equation below.

$$S=\frac{25400}{CN}-254$$

(4)

Here, S is the maximum water retention parameter (mm) and CN is the weighted curve number that is calculated while using the equation below, based on the storm-event method [82], which is a data-derived value that varies according to the rainfall [58,81].

$$CN=\frac{A_{1}C{N}_1+A_{2}C{N}_2+A_{3}C{N}_3\dots +A_{n}C{N}_n}{A_1+A_2+A_3\dots +A_n}$$

(5)

Here, A₁, A₂, A₃, …, and A_n are the areas of hydrological groups that a given land cover falls into, and CN₁, CN₂, CN₃, …, and CN_n are the corresponding curve numbers.

2.4. Hydrologic Soil Group (HSG) Classification

The HSG of the study area was adopted from the classification that was employed in the USA developed by Cronshey [83], which is based on the infiltration rate controlled by the soil profile (

Table 3. Classification of hydrologic soil groups.

Soil Group	Characteristics
A	Low overland flow potential, high minimum infiltration capacity even when thoroughly wetted (> 0.76 cm/h) Deep, well, to excessively drained sands and gravel.
B	Moderate minimum infiltration capacity when thoroughly wetted (0.13–0.76 cm/h) Moderately deep to deep Moderately to well drained Moderately fine to moderately coarse grained (e.g., sandy loam).
C	Low minimum infiltration capacity when thoroughly wetted (0.13–0.38 cm/h) Moderately fine to fine grained soils or soils with an impeding layer (fragipan).
D	High overland flow potential Very low/minimum infiltration capacity when thoroughly wetted (<0.13 cm/h) Clay soils with high swelling potential Soils with permanent high-water table Soils with a clay layer near the surface Shallow soils over impervious bedrock.

Source: Dingman [84].

). Technical Release 55 (TR-55) of urban hydrology for small watersheds presents the hydrologic soil group type according to the surface soil texture. Sand, loamy sand, and sandy loam belong to type A of HSG, silt loam and loam belong to type B of HSG, sandy clay loam belongs to type C of HSG, and clay loam, silty clay loam, sandy clay, silty clay, and clay belong to type D of HSG (

Table 4. Runoff curve numbers for antecedent moisture condition (AMC II) for the study area.

Land Use Type	Texture	HSG	CN	Area	LU%
Bare	Clay	D	89	160
Bare	Clay loam	D	89	370
Bare	Loam	B	82	2
Sub Total				532	0.96
Built-up	Clay	D	98	71
Built-up	Clay loam	D	98	257
Sub Total				328	0.59
Cultivation	Clay	D	91	10,117
Cultivation	Clay loam	D	91	12,262
Cultivation	Loam	B	81	72
Cultivation	Sandy clay	A	72	5
Cultivation	Sandy clay loam	C	88	6
Sub Total				22,462	40.53
Forest	Clay	D	83	2677
Forest	Clay loam	D	83	9470
Forest	Sandy clay loam	C	77	6
Sub Total				12,153	21.93
Grass	Clay	D	80	532
Grass	Clay loam	D	80	1011
Grass	Sandy clay	A	39	3
Sub total				1546	2.79
Shrub	Clay	D	89	4770
Shrub	Clay loam	D	89	13,490
Shrub	Loam	B	79	2
Shrub	Sandy clay	A	68	13
Shrub	Sandy clay loam	C	86	4
Sub Total				18,279	32.98
Water	Clay	D	100	37
Water	Clay loam	D	100	82
Sub Total				119	0.21
Total				55,419	100

).

2.5. Data Analysis

Model input data of historical daily rainfall and temperature data of 1981–2018 were used to generate climatological data sets for input into the model. The mass curve method [85,86] was employed to verify the consistency of rainfall values. As indicated by Garg [86], all input data were rendered into a grid format, and the isohyetal method was used to interpolate the areal rainfall. Table 1 summarizes the data sources and data collection methods for the input variables. Furthermore, the Mann–Kendall test and linear regression test have been applied to analyses of historical trends, while Pettitt’s test was applied to detect the change point of the runoff response dataset. Evapotranspiration-runoff, rainfall-runoff, and land cover-runoff relationships were also tested while using the linear regression method.

3. Results

3.1. Topographical Classification and Elevation Map of the Catchment

The study area has been classified into six basic slope classes: 0–15% (flatland), 15–30% (gentle), 30–45% (intermediate), 45–60% (slightly steep), 60–75% (steep), and 75–90% (very steep). Figure 3a,b shows that the catchment elevations vary between 1444 m in the eastern part to 3699 m in the northern part of the catchment. The minimum slope is between 0% and 71.88%, and the average slope change is from 17.47% to 38.62%, the highest being 71.88% (Figure 3b). Figure 3a indicates that the areas of the low elevation (1410–1452.22 m) may be flooded during heavy rainfall periods.

3.2. LULC Dynamics of the Dry Afromontane Forests of HGK State Forest

The most dominant LULC types in 1986, 2001, and 2018 were shrubland, cultivated land, forest land, grassland, water, and built-up area, respectively (

Table 5. Land use and land cover (LULC) dynamics of the dry Afromontane forest.

LULC class	Percentage of Land Cover (%)			Change
	1986	2001	2018	1986–2018	1986–2001	2001–2018
Bare	0.00	0.24	0.94	+0.94	+0.24	+0.70
Built-up	0.05	0.30	0.58	+0.53	+0.25	+0.28
Cultivation	42.02	36.79	39.56	−2.46	−5.24	+2.77
Forest	4.72	23.87	21.40	+16.68	+19.14	+2.47
Grass	2.25	2.73	2.75	+0.49	+0.48	+0.02
Shrub	48.51	33.49	32.18	−16.33	−15.02	−1.31
Water	2.44	2.59	2.60	+0.16	+0.15	+0.01
Total	100	100	100

). In 1986, the highest portion of the total area was covered by shrubland (48.51%), followed by cultivated land (48.54%) and 4.72% of high forest (Figs. 4). Grassland (2.25%), the water body (lake Hashengie: 2.44%), and built-up area shared 0.05% of the total area (Table 5). In the year 2001, cultivated land (36.79%), followed by shrubland (33.49%), were the dominant LULC types, in 2018, the largest portion of the land was cultivation area (39.56%), followed by shrubland (32.18%) and forest land (21.40%) (Table 5). Grasslands, bare land, and built-up area accounted for 2.75%, 0.94%, and 0.58%, respectively (Figure 4). Forest land cover increased from 4.7% to 21.4% during 1986–2018 in the study area at the expense of shrubland, owing to the community-based environmental protection effort in the region and the national state policy. This finding is supported by Berhane et al. [39], who reported that forest cover of 35.1%, followed by shrubland (30.14%), occupied the largest portion of the land between 1985 and 2015, indicating an increment of forest cover (714 ha) and grassland (75 ha) in the Hugumburda forest of northern Ethiopia [13,39,47]. However, Haregeweyn et al. [44] reported the major increments of cultivated land by 15.4% and settlements 9.9% at the expense of shrubland and grazing lands over the period of 1976–2003 in the Gilgel Tekeze catchment in the highlands of Northern Ethiopia [44,45].

The overall statistical accuracy assessment was performed on the resulting classified imagery. This process involves generating a set of points in the classified imagery and comparing them with actual points on the ground through field work with an accuracy of 83–90% and a Kappa coefficient of 0.82–0.92 of these LULC maps.

3.3. Hydrological Soil Group (HSG)

Hydrologic soil groups (HSG), along with land use, management practices, and hydrologic conditions, determine the soil cover complexes and their associated runoff curve numbers. The study area consists of ‘A’, ‘B’, ‘C’, and ‘D’ hydrologic soil groups (HSG) (Table 3 and Table 4, Figure 5). HSG ‘D’ covered the highest portion throughout the area, followed by ‘B’, ‘C’, and ‘A’ groups comprising 99.84%, 0.13%, 0.02%, and 0.01% of total area, respectively. This result indicates that most of the soils in HGKF have very low or minimum infiltration capacity when thoroughly wetted, which makes them vulnerable to runoff.

Similar to this, the study by Berhanu et al. [87] also confirmed that HSG classified HSG-A dominated the areal coverage (with 48.2%), followed by HSG-B (30%) and HSG-D (21.6%) for all Ethiopia. In line with these findings, a study from the Kharadya mill watershed, India, by Ningaraju et al. [88], found that HSG of ‘A’, followed by ‘D’ and ‘B’ type, consisted of 58.63%, 38.26%, and 3.11% of the total area in the watershed, respectively.

3.4. CN Number of the Area Based on Their DEM

Analysis of the CN number (Table 4 and Figure 6) indicates that the type of land use will determine the nature of the relationship between runoff CN and land use change. CN is a relative measure of retention of water by a given soil vegetation complex and it takes values from 0 to 100. The runoff generation relies on the CN values, which are a function of AMC, slope, soil type, and land use. The CN values of the HGK state forest, varying from 39 to 100, reflect the runoff potential. Low CN values mean that the surface has high potential to retain water, whereas high values indicate that the rainfall can only be stored to a limited extent. The curve number (CN) is a hydrologic parameter that is used to describe the stormwater runoff potential for drainage areas, and it is a function of land use, hydrological soil type, DEM, and soil moisture. The weighted values of curve number for the AMC conditions of the study area are calculated as per SCS. The CN values that are shown in Figure 6 and Table 4 were used to estimate the runoff for 38 years of the watershed. The monthly and annual runoff values are obtained from the daily runoff results. Land use and land cover has influenced surface runoff generation in a given area to a greater extent. Land-use type explains the difference in surface runoff CNs significantly better when compared with the hydrological soil group. The increased surface runoff causes different types of erosion, as rill and sheet erosion intensify, while gullies expand, leading to reduced soil depth and water retention capacity of the soil. This finding is in line with previous findings [13,35,54] of studies that were conducted in northern Ethiopia. Moreover, it verified that areas with less vegetation cover are exposed to high surface runoff and low retention capacity. In the study area, there is an obvious increase of both direct surface runoff and runoff coefficient over time from 1986 to 2018.

3.5. Annual Runoff Coefficient Potential

The temporal distribution of annual runoff depths is displayed in Figure 7 and Figure 8, where the mean annual potential discharge of 31.6 mm per year was estimated during the period from 1981 to 2018. A large variation from the minimum annual potential runoff of 17.6 mm/year in the forest to the maximum amount of 40 mm/year in some areas that are usually occupied by rainfed agricultural activities was observed. This means that the runoff derived by SCS-CN method is a function of the runoff potential, which can be expressed in terms of the runoff coefficient (ratio between the runoff (Q) and rainfall (RF)) categorized in three basic classes of severity to runoff and soil erosion, i.e., high (> 35%), moderate (20–35%), and low (< 20%). This result is consistent with the study at [89] Kali Watershed in India. Furthermore, another study by Blokhuis [90], noted that the lowest annual runoff potential was observed in the forest area, which dominates the reddish sandy clay soils and it produced the minimum potential of runoff, i.e., 13 mm/year to 32 mm/year. This was due to the fact that sandy soil has high infiltration rates and also because the dense canopy of trees increases the rainfall interception and water losses through evapotranspiration.

3.6. Impact of Rainfall and PET on Runoff Discharge

Table 6. Impact of rainfall and potential evapotranspiration on runoff coefficients.

Variable	Coefficients(B)	Std. Error	T	Sig.
Constant	0.290	0.80	3.640	0.001
Rainfall	0.837	0.061	10.754	0.000
PET	−0.240	0.095	−2.533	0.016
Regression (Sum of Squares)			814	0.00
Residual (Sum of Squares)			133	0.00
F-test				63.482

presents the result of linear regression models. The F-test estimate of the linear regression models is significant at less than the 1% probability level, which indicates the sufficient modeling accuracy (F = 63.48; p < 0.000%). Of the hypothesized three explanatory variables, two of them were found to significantly determine the probability of precipitation.

Rainfall: The rainfall-runoff relationship in the HGK forest area was analyzed while using the SCS-CN method. The regression result shows that the coefficient of rainfall is positive and significant at less than 1% probability level for discharge. All other factors being kept constant, the increase of rainfall results in the probability of discharge increase by 83.7%. This implies that rainfall has a higher probability of influencing the discharge. This finding is similar to those of Refs. [13,35,54,55] in northern Ethiopia.

PET: The regression result shows that the coefficient of PET is negative and significant at P < 0.05 probability level for discharge. Maintaining all other factors as constant, an increase in PET causes a 24% decrease in the probability of discharge. This implies that PET has a lower probability of influencing discharge than rainfall. This finding is similar to that of Ref. [91] in the US watershed.

3.7. Relationship of Runoff with Climatic Factors

The correlation between runoff and climatic factors, such as RF, PET, Min T, and Max T, of the fragmented forest catchment attributes was performed (Appendix A;

Table 7. Pearson’s coefficient of correlation between discharge and climatic parameters.

		Year	Rainfall	Q	PET	Min To	Max To
Year	Pearson Correlation	1
	Sig. (2-tailed)
Rainfall	Pearson Correlation	.361*	1
	Sig. (2-tailed)	.026
Q	Pearson Correlation	.263	.863**	1
	Sig. (2-tailed)	.110	.000
PET	Pearson Correlation	.032	−.076	−.264	1
	Sig. (2-tailed)	.848	.648	.021
Min T	Pearson Correlation	−.044	−.271	−.061	−.385*	1
	Sig. (2-tailed)	.795	.100	.715	.017
Max T	Pearson Correlation	.032	−.076	−.264	1.000**	−.385*	1
	Sig. (2-tailed)	.848	.648	.109	.000	.017
N		38	38	38	38	38	38

* Correlation is significant at the 0.05 level (2-tailed); ** Correlation is significant at the 0.01 level (2-tailed).

). The results show that rainfall is positively and strongly correlated (P < 0.000), while PET and Max T are negatively and significantly (P < 0.05) correlated with discharge. The correlation analysis revealed a strong relationship (R² = 78.4; P = 0.000) between the modeled point-data derived rainfall and runoff estimated by the SCS-CN method.

3.8. Annual Rainfall Dynamics in Relationship to Discharge and PET

The temporal variability of runoff due to climate variability is presented in Figure 8, Table 7, and Appendix B, showing the extent of annual rainfall dynamics in relationship to Q and PET. The average change in PET increased by a 1% during 1981–2018.

3.9. Linear Regression Analysis Between Runoff and Rainfall Intensity

The runoff varied from 25.68 mm/h to 56.4 mm/h for events of high intensity (storm intensities: 52.11–81.77 mm/h) and from 11.19 mm/h to 24.69 mm/h for events of low intensity (storm intensities: 26.07–45.38 mm/h). With high-intensity rainfall, the runoff coefficients ranged from 25.3% to 47.21% (Appendix B). Figure 9 shows the strong positive correlation between runoff (dependent variable) and rainfall intensity (independent variable). The regression equation allows for an estimate to be made regarding the runoff.

4. Discussion

Land use change has a significant impact on the hydrological and ecological processes of the study area of watershed. Our findings indicate a change in LULC between 1986 and 2018, of about 32.33% of the study area, where an increase in the forest land by 16% and bare land by 0.94% was estimated, totaling about 9000 ha, mainly at the expense of the shrubland. Similar studies in northern Ethiopia [24,39,55,70] also showed significant LULC alterations and transformations since the late 1950s.

The dynamics of LULC change depicts a crucial environmental alteration, which has pronounced impacts on human livelihoods, particularly the developing countries, such as Ethiopia. This research is the combination of an empirical land use change model and an event scale, rainfall-runoff model to quantify the impacts of potential land use change on the storm-runoff generation in the HGK state forest area of northern Ethiopia. Land-use and land-cover changes direct impacts on the hydrological cycle [7,92], by causing floods, droughts [93], and changes in river and groundwater regimes [21,22], and they can affect water quality. Therefore, in this study, covering the time period from 1986 to 2018, we employ the SCS-CN model either static (one land use map) or dynamic (land use updates) representation of LULC. The two scenarios provided annual data on LULC changes that served as an input for the dynamic model runs. In the HGK state forest watershed, dynamic LULC changes were revealed during the study period. The major change was the increase of forest cover, which was due to the reforestation and forest conservation program that was implemented in northern Ethiopia since 2000. The bare land area also showed a slight increase at the expense of the area of open shrubland. This LULC dynamics appears to have affected the stream flow of the watershed. During the period between 1981 and 2018, the total discharge of the land scape increased at a rate of 0.78 mm per annum, whereas rainfall only increased at a rate of 0.27 mm per annum. This increase in the discharge was caused by increasing the bare land area at the expense of slight increase of rainfall, increased transpiration losses due to the increased tree cover and a decreased contribution from the base flow, as revealed by the analysis of extreme low flows (

Table 8. Characteristics of monthly and annual rainfall P and stream discharge flow (Q) in Hugumburda Grat Kahisu state forest Northern, Ethiopia (1981–2018).

		January	Februay	March	April	May	June	July	August	September	October	November	December	Annual Average
Discharge flow (Q)	Mean	0.053	0.153	0.518	0.888	0.784	0.15	2.89	2.99	0.382	0.26	0.07	0.07	0.78
	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Maximum	6.36	13.13	32.50	45.75	40.06	27.03	56.40	37.69	18.63	27.04	13.18	15.29	56.40
	Std. Error of Mean	0.008	0.03	0.075	0.11	0.09	0.05	0.19	0.12	0.04	0.05	0.018	0.0208	0.026
Rainfall	Mean	0.51	0.76	1.56	2.44	2.11	0.65	6.42	7.07	1.70	1.02	0.47	0.38	2.11
	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Maximum	18.92	30.59	53.06	68.04	61.69	48.68	81.77	60.82	37.26	48.68	30.67	31.97	81.77
	Std. Error of Mean	0.054	0.09	0.17	0.23	0.200	0.098	0.37	0.29	0.1242	0.12	0.064	0.061	0.057

). The detected increase in the surface runoff can be expected given the significant LULC changes that were revealed in the watershed. Girmay et al. [13] also showed that changes in land use/cover decreased the water storage capacity of soils by a factor of 1.63 with a corresponding increase in the surface runoff by a factor of 2.7 at Gum Selassa, northern Ethiopia. The increased evapotranspiration losses and the decline in base flow are both associated with changes in the land cover of the watershed and/or watershed degradation. The other contributing factor for the decreased stream flow, particularly during the dry season, has been the increased water abstraction to be expected from the increased human and livestock populations in the watershed. The study by Gebremiceal et al. [94] in the Geba catchment of northern Ethiopia that with an increase of agricultural land by 42% and a decrease of natural vegetation cover by 36%, the average median monthly flow during the wet season increased by 4% and the dry months decreased by 23%, after continued LULC changes.

In the Chemoga watershed in northern Ethiopia, Bewket and Sterk [47] showed that the total stream flow decreased at a rate of 1.7 mm/y during the dry season and attributed development partially to changes in LUC and/or to land degradation of the watershed [47,69]. This could explain the decrease in base flow in the Chemoga watershed during the dry season [47]. Similar to this finding, in China, a 3% increase in streamflow for the whole watershed from 1984 to 2010 was observed due to the long term impact of LULC [95]. Appendix B describes that, during the 1981–2018 period, the stream flow increased in accord with a slight increase in rainfall, however the rate of increase in the stream flow was less than that of the rainfall (Figure 9 and Figure 10; Table 8). This suggests the significance of exhaustion of the groundwater, which is to be explained by increased water use by the increased vegetative cover in the watershed that occurred between 1981 and 2018, and/or the increased water abstraction. A similar study [47,69] observed that a large volume of surface runoff occurs during storm events.

5. Conclusions

Our study provides evidence of the significant LULC change in the watershed between the years 1986, 2001, and 2018. The decrease in shrubland and grassland was accompanied by an increase in forest land, bare land, and built up areas. These changes may influence the overall ecosystem functioning. Our findings indicate that these LULC changes have a significant effect on the generation of surface runoff, namely direct surface runoff, runoff coefficient, and storage capacity of the soil in the study watershed in northern Ethiopia.

The study area was classified into three hydrologic soil groups based on the results of soil classification and land use types, being texturally dominated by 99.84% of clay loam type HSG of “D” class. This indicates that the study area has a minimum infiltration capacity when thoroughly wetted, and thus it is increasingly vulnerable to runoff. The runoff for different LULC classes revealed that the wet season flow increased for the most recent year, which is attributed to the conversion of shrublands to bare land and agriculture.

Rainfall was shown to have a positive and strong correlation (P < 0.000) with runoff, while PET and Max T were negatively (P < 0.05) correlated with discharge. Northern Ethiopia dry Afromontane ecosystems are characterized with an erratic rainfall and high potential evapotranspiration, which is very susceptible to drought. Therefore, this finding concludes the long-term LULC change influence the hydrology of the entire dry Afromontane forest landscape. There is a need to devise mechanisms for increasing its resilience to adverse hydrological changes emanating from LULC changes since the HGK state forest is characterized as a semiarid land ecotone with high evapotranspiration climatic regime, scarcity of surface water, and rainfed agriculture. In this regard, the results of the SCS-CN model provide some relevant information for land use planning, and watershed and water resource management. In particular, considering the high runoff potential revealed, soil water conservation structures (e.g., percolation ponds) for the upper catchment are recommended. Continued change in the LULC is becoming a serious threat to the Hugumburda Grat Kahisu reservoir forest area. The LULC change should be controlled, and measures need to be likewise taken for the stabilization of the negative land cover changes.

In addition to the runoff response focused in this study, the response of sediment dynamics to LULC change and its relevance to the landscape management is useful to consider in future studies. Soil water conservation structures, such as dip trench on farm land, percolation ponds, series pond, small scale check dam, subsurface ditch, and percolation tanks, are suggested policy implications in the watershed for sustainable water resource development of the Hugumburda Grat Kahisu reserved state forest of northern Ethiopia. Further investigation is also suggested on the climate scenario-based modeling of hydrological processes that are influenced by land use changes, which can improve the understanding of hydrological variability of dryland forest ecosystems.