Intensity-Based Quantification of Land Use Dynamics in Biodiversity Hotspot Area of Southern Rift Valley Basin, Ethiopia

1. Introduction

Among the global challenges in the 21^st Century, LUC is a rapidly growing environmental problem [1], negatively altering the entire functions of Earth [2]. Recent studies [3,4,5,6] confirm that the pressures of LUC inflicted on Earth’s ecological functions have been pervasive and intensive, resulting in unprecedented deterioration in global ecosystems and biodiversity over the last few decades. In this study, LUC refers to the conversion, modification, and management of landscapes/natural ecosystems/ primarily by human activities [7]. It is a long-lasting phenomenon derived from the effects of human and natural factors; essential for socioeconomic development and ecological balances [8]. Nevertheless, the widespread and intensive LUC, predominantly, massive transitions from natural ecosystems to human land uses have been observed over the last few decades because of increasing land and natural resources demand for agricultural activities [7,9], urbanization and industrialization [10] accompanied by environmentally harmful development strategies [3]. Despite the short-term contribution to socioeconomic enhancement [11]; these land transition processes are strikingly jeopardizing the earth’s ecosystems and biodiversity [4,12]; with estimated loss of >63% of ecosystem services [13]. As a result, people worldwide, largely in developing countries have suffered from catastrophic environmental and socioeconomic challenges: climate change [2] and chronic poverty and food security [1].

Following the multitudinous impacts of human-triggered LUC and degradation; sustainable conservation and use of biodiversity have become the agendas of researchers, planners and governments and measures have been undertaken from local to global levels [3]. Designating biodiversity hotspot sites as protected areas(PAs) is one of the measures taken [6]. Biodiversity hotspots are natural ecosystems, characterized by a significant concentration of flora and fauna diversity with a high level of endemism but susceptible to anthropogenic impacts [14]. PAs including National parks are conceived as the most important strategy for sustainable biodiversity conservation with marked improvement in ecosystems and the services they provide [6,15,16,17]. They have been also recognized as productive places for restoring critically/endangered species and maintaining ecological integrity [12]. Generally, cognizant the potential advantages, the area and number of PAs have shown a radical increment; today there are over 155,584 PAs worldwide [18]. However, due to rapidly growing human encroachments, PAs worldwide are under pervasive LUC and multifaceted ecological problems including massive habitat conversion and degradation, biodiversity loss, and species extinction [14,16]. The magnitude and adverse impacts are more deleterious in PAs found in developing countries because of cumbersome demographic, socioeconomic, political, and administrative challenges [12,17], specifically in Sub-Saharan Africa region where more than 82% of the PAs are under the state of deterioration/failures with 10% only in success/recovery [10].

In this region, Ethiopia has been recognized as one of the countries rich in biodiversity with a wide range of ecosystems [6]. Nevertheless, over the last few decades, due to rapid and intensive land conversion from forest, grass and woodlands to agricultural and settlement areas [9,19] and massive wood production for multitudinous socioeconomic benefits [20] coupled with unsustainable land management practices [21], the country has faced a substantial cascading of ecosystems and biodiversity; subsequently with serious land degradation and socioeconomic crises: prolonged poverty and food insecurity [22]. These situations coupled with aggressive population growth have forced anthropogenic pressures toward biodiversity-enriched/hotspots/areas and markedly worsened the deterioration of biodiversity and ecosystem services [6,23]. To tackle the daunting impacts of environmental degradation and preserve biodiversity-enriched areas from anthropogenic pressures, the government of Ethiopia has established different PAs including twenty national parks since the 1960s [24]. However, like other developing countries, PAs in Ethiopia including NNP (where this study was undertaken) have been under extensive LUC and degradation due to significant habitat conversion to agricultural and settlement areas, livestock grazing, and fuel and construction wood production [25,26,27]. Moreover, LUC and the resulting impacts such as habitat and biodiversity degradation are aggravated by the expansion of invasive species, recurrent wildfires and weak administrative measures [6]. There by, currently, almost all PAs have lost their functionality, the survival of many wildlife has reached a critical stage and they stand on the verge of failure [28].

Among gazetted PAs, NNP is the most biodiverse [flora and fauna} with a large number of endemic and critically/endangered species in Ethiopia [29] and multifold ecological and socioeconomic advantages [6]. Considering the huge potential for biodiversity conservation, it has been registered as one of the major African floristic areas and global biodiversity sites [27]. Nevertheless, NNP is currently one of Ethiopia’s Pas that extremely vulnerable to human-induced LUC and degradation where natural resources are the primary source of livelihood, economic and sociocultural practices for rapidly growing millions of semi/pastoralists and urban dwellers living in and around the park [30].

In general, the pressing loss of functionality of Ethiopia’s PAs is mainly attributed for anthropogenic-driven habitat conversion and degradation. Thus, to monitor the ongoing changes, restore degraded ecosystems and direct the future management and anthropogenic use plans toward sustaining the functionality of PAs; analyzing and mapping of the spatiotemporal dynamics of LUC is imperative. Although several studies have been conducted to address the problems of PAs in Ethiopia, the attention given to analysis and mapping LUC dynamics is scant [25,26,27,30,31,32,33,34]. Moreover, the LUC observed at any spatiotemporal level is characterized by interdependent and complex steps, transition patterns, and processes and also has a relationship with changes in other spatiotemporal levels [35]. Analyzing the entire LUC processes comprehensively and concurrently via a robust model with rigorous procedures is, therefore, necessary for sound conclusions, informed decisions and enhanced management plans [35,36]. Specifically, a detailed assessment and precise measurement of land transition patterns and process for PAs is imperative to show the extent of the change to local people and improve their awareness, identify vulnerable ecosystems, examine key transition and stationary patterns, and develop ecosystem-targeted restoration and protection strategies. However, the previous studies in PAs including NNP did not assess the complete LUC characteristics and have failed to provide full information about change processes and patterns. Even at the national level most of the studies have focused on some aspects of LUC [8]. Thus, the change analysis methods applied in these studies were not enough to accommodate all the necessary steps, patterns and processes of LUC and considerably excluded the effects of length of time intervals, areal extent, and strength of the changes, which markedly discourage, planners, policy and decision-makers from informed decision and policy-making processes for sustainable management and use of biodiversity.

Furthermore, due to considerable deviations in anthropogenic pressures, and ecological setups of PAs, the magnitude and transition patterns of LUC have deviations as shown by earlier studies in Ethiopia and other developing countries [17,37,38]. Thus, for a better understanding of the change extent and also to make actionable strategies considering conditions at local level including, demographic and socioeconomic variables; analysis of LUC for each protected area is extremely crucial and recommended by many researchers [25,26,37,39]. Therefore, based on the mentioned gaps above, this study aimed to quantify the entire LUC characteristics in NNP at three levels(interval, category and transition) using a compressive mathematical model; intensity analysis [35] from 1986 to 2002 (First Time Interval /FTIS/), and 2002 to 2020 (Second Time Interval/STIS/). The study also answered the following questions; in which the study’s time interval, the overall change intensity is fast/slow? which land types’ gain and loss are active/dormant? Why? Which land types are most intensively increased or decreased? Why? which transitions are most intensively avoided/targeted by a given land use type? which transitions in which land types are stationary across the study time interval? Finally, what is the implication of LUC observed in NNP for biodiversity, wild animals’ habitat, and livelihoods?

2. Methods and Materials

2.1. Study Area

The study area, NNP is located East of Arba Minch town at the upper part of Segen river catchment in the Southern part of Ethiopian Rift Valley Basin. With the area of 41400 ha and altitudinal range from 1,086 to 1642 m asl; NNP is situated between 5⁰ 51'- 6⁰ 05' N and 37⁰ 32'-37⁰ 48' E (Figure 1). The climate condition is a hot tropical type climate mainly characterized by inadequate and unpredictable trend of rainfall [29]. According to the data from the Ethiopian Meteorological Agency(EMA, 2022); the total annual rainfall ranges from 622 to 1177 mm with an average value of 888.38 mm. The average annual minimum and maximum temperature are between 16 and 20 ⁰C, 30 and 35 ⁰C, respectively (Figure S 1).

NNP is occupied by markedly diversified terrestrial and aquatic ecosystems. The terrestrial ecosystem ranges from evergreen natural riverain forest at the Western and Eastern edges; woodland and bushlands on the volcanic hill (locally named ‘God’s Bridge’) to savannah grassland on the plain of Eastern part [30]. NNP is also known for its hosting large number of endemic, and endangered species in Ethiopia [40,41]; protecting these species, specifically Swayne’s Hartebeest and plain Zebra is the main objective while it was primarily gazetted in 1974. As estimations show the park accounts for nearly 20% of the national biodiversity with more than 90 mammals, 350 bird and 800 – 1000 plant species [24,29]. Apart from ecological advantages, for a long time, the park has served as vital asset for the socioeconomic activities of many urban and rural populations. However, recently, due to population growth, the rapidly expanded and intensified socioeconomic interventions: cultivated and settlement areas, livestock grazing and wood production jointly with administrative and political complications exposed the park to widespread LUC and degradation, lead to significant habitats and biodiversity loss [29].

2.2. Data Source, Image Pre-Processing and Land Use Classification

In this study, a medium resolution(30 m) Landsat imagery of TM (1986), ETM+ (2002) and OLI (2020) that projected to UTM projection, zone 37 and WGS84 datum, were downloaded from USGS’s database(C2L2) (http://glovis.usgs.gov) for land use classification and change analyses. For the detailed descriptions of imageries used in this study (see Table S 1). In addition, ancillary data such as Google Earth, SAS Planet, and toposheet were utilized in this study.

Using ERDAS Imagine_ 2014; the necessary image preprocessing activities: geometric correction with a 1:50,000 topographic map (1986) and atmospheric correction were performed before the land classification process was undertaken. To identify the major land use types: forest land (FL), woodland (WL), bush/shrubland (BS), grassland (GL); water area (WA), and cultivated land (CL); interviews with elders living in and around and workers of NNP, intensive field surveys, and review of previous studies [27,30,34] were made. The reference and training points for 2020 were collected via GARMIN-GPS during field observations, while, the topographic map, Google Earth and NDVI values were utilized for 1986 and 2002. The Image classification process was undertaken in ArcGIS_10.8.1 using unsupervised (ISODATA) and Support Vector Machine (SVM) methods. SVM is a roust for land classification with less sensitivity to the effects of noise and unequal number of training points compared with other Supervised classification algorithms [42,43]. During the classification process, misclassification was observed between bush and shrubland, and cultivated and settlement areas. This problem was solved by amalgamation of shrubland with bushland, and cultivated land with settlement areas; which was applied and suggested by several researchers [44,45,46]. The nomenclature and descriptions for land use types are presented in (Table S 2)

2.3. Sampling and Accuracy Assessment(AC) Methods

AC is necessary and essential to check the quality of classified images [25,47]. Accordingly, the Kappa coefficient (Khat) was applied to check the accuracy of all classified land use maps in this study (equ. 1); the most widely used index for accuracy assessment in LUC research [48]. In addition, the total sample size(number of ground truth points) for AC was specified based on equation (2), and sample size for each land type was also obtained based on the proportional stratified sampling method with consideration of the rule of thumb set by [48] and presented in(Table S 3).

$$\mathrm{Khat}={\displaystyle \frac{N{{\displaystyle \sum }}_{i=1}^r{X}_{ii}-{{\displaystyle \sum }}_{i=1}^r{X}_{it}\ast X_{ti}}{N^2-{{\displaystyle \sum }}_{i=1}^r{X}_{it}\ast X_{ti}}}$$

(1)

Where; r is the number of rows in the matrix, x_ii is the number of observations in row i and column i, x_it and x_ti are the marginal totals of row i and column i respectively, and N is the total number of observations.

$$n={\displaystyle \frac{B{\prod }_i \left(1-{\prod }_i\right)}{b_i^2}}$$

(2)

Where; n is the total sample, B is determined from a chi-square table with 1 degree of freedom and 1 - $\mathit{\alpha }$/k, (for this study, χ2 _{(1, 0.01)} = 6.63. k is the number of land types (k = 6); a is the precision error tolerance (0.05), b is desire precision (96%) and ∏_i is the proportionate size of the i^th land type out of k types.

The results of AC are presented in (Table S 4). The OA and Khat levels for the classified land use map of 2020 were 96.22 and 95.20%, respectively. Similarly, 92.92 and 94.97% of OA and 90.84 and 93.52% of Khat coefficients were attained for land use maps of 1986 and 2002, respectively. These results indicate the agreement across all land types between the ground truth points and classified maps for all years is strong enough(Khat > 80%), and the classified maps are adequately sufficient to post classification operations including change detections [48].

2.4. Intensity Analysis (IA)

The transition matrix was used to examine the gross gain, gross loss, persistent, swap and net changes. Although this matrix widely applied for land use and land change analysis, it doesn’t consider the necessary steps and spatiotemporal relationship of land change processes during the analysis [8] and indicate the change intensities and reasons why the size of changes are observed among the land types [36]. Considering its ability to characterize and show the entire LUC characteristics compressively, IA was applied to examine the LUC processes at interval, category and transition levels. [35]. It also indicates stationery trends and rational reasons for each land transition in all levels [36,49]. The detailed procedure and steps of IA framework located in Figure 2.

Interval level of intensity analysis(ILI) is 1^st level of intensity analysis that used to computes the magnitude and rate of LUC for each time interval. At this level, IA helps us to get information about in which time interval the annual rate of overall change is relatively slow versus fast by comparing the observed change intensity with a uniform change intensity that would exist if the annual changes in the study area were uniformly distributed throughout the study period [35]. Accordingly, equation (3) was used for calculating the observed annual rate/intensity/ of LUC for each time interval and Equation (4) for quantifying the uniform annual change intensity that would exist if LUC process were uniformly distributed over the spatiotemporal extent of the study. In other words, if values of St were equal for all t, then those values of St would equal U [35].

Category level of intensity analysis(CLI), 2^nd level of IA which applied to analyze the intensity of annul gain and loss among the land use types. It indicates which land types exhibited dormant/active change in a given time interval comparing with the intensity of uniform the change that could observed among the land types over the study area [35,36]. The intensity of annual gain and loss for each land type were calculated (equ.5) and ( equ. 6) [35], respectively.

Transition level of intensity analysis(TLI) is the 3^rd level of IA used to examine the intensity of transitions(gain and loss) of a land type across other land types available for each time interval [8]. It is important to analyze the transition from a particular land type m to different types of n and compare with a uniform intensity of transition that would exist uniformly among the land types available in the transition and to identify which land types are intensively avoided/targeted during the transition period [35]. Equation (7) was applied to calculate the intensity of the annual transition from type i to n during the interval [Yt, Yt + 1], where I # n and equation (8) for the uniform intensity of the transitions from all non-n types to type n (Yt, Yt + 1). If R_tin is greater than W_tn, then the gain in n displaces i; if R_tin is less than W_tn, then the gain in n does not affect i [35]. Likewise, the losing component of TLI is important to quantify the size of land transitions from the losing land types. Equation (9) was employed to compute the observed intensity of the annual transition from type m to type j during the given time interval (Yt, Yt + 1) relative to the size of type j at time t+1, where j # m and equation (10) used for the uniform intensity of annual transition from type m to all other non-m types (Yt, Yt+1) relative to the size of all non-m types at time t+1. If Q_tmj is greater than V_tm, then j targets the loss in m; if Q_tmj is less than V_tm, then j avoids the loss in m [35]. The following equations are the statistical formulas of IA, and the nations for symbols presented(Table 1) adopted from [35].

$$S_t={\displaystyle \frac{{{\displaystyle \sum }}_{j=1}^J\left[\left({{\displaystyle \sum }}_{i=1}^J{C}_{tij}\right)-C_{tij}\right]/\left[{{\displaystyle \sum }}_{j=1}^J\left({{\displaystyle \sum }}_{i=1}^J{C}_{tij}\right)\right]}{Y_{t+1} -Y_t}}x100$$

(3)

$$\mathrm{U}={\displaystyle \frac{\left.{{\displaystyle \sum }}_{t=1}^{T-1}\left\{{{\displaystyle \sum }}_{j=1}^J\left[\left({{\displaystyle \sum }}_{j=1}^J{C}_{tij}\right)-C_{tij}\right]\right.\right\}/\left[{{\displaystyle \sum }}_{j=1}^J\left({{\displaystyle \sum }}_{i=1}^J{C}_{tij}\right)\right]}{Y_T-Y_1}}x100$$

(4)

$${\mathrm{G}}_{\mathrm{tj}}={\displaystyle \frac{\left[\left({{\displaystyle \sum }}_{i=1}^J{C}_{tij} \right)-C_{tij}\right]/\left(Y_{t+1}-Y_t\right)}{{{\displaystyle \sum }}_{i=1}^J{C}_{tij}}}\mathrm{x}100$$

(5)

$${\mathrm{L}}_{\mathrm{ti}}={\displaystyle \frac{\left[\left({{\displaystyle \sum }}_{j=1}^J{C}_{tij}\right)-C_{tij}\right]/\left(Y_{t+1}-Y_t\right)}{{{\displaystyle \sum }}_{j=1}^J{C}_{tij}}}\mathrm{x}100$$

(6)

$$\mathrm{tin}={\displaystyle \frac{\raisebox{1ex}{$C_{tin}$}\!\left/ \!\raisebox{-1ex}{$\left(Y_{t+1}-Y_t\right)$}\right.}{{{\displaystyle \sum }}_{j=1}^J{C}_{tij}}}\mathrm{x}100$$

(7)

$${\mathrm{W}}_{\mathrm{tn}}={\displaystyle \frac{\begin{array}{c}\\ \raisebox{1ex}{$\left[\left({{\displaystyle \sum }}_{i=1}^J{C}_{tin}\right)- C_{tnn}\right]$}\!\left/ \!\raisebox{-1ex}{$\left(Y_{t+1}-Y_t\right)$}\right.\end{array}}{{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{J}}\left[\left({{\displaystyle \sum }}_{i=1}^J{C}_{tij}\right)-C_{tnj}\right]}}\mathrm{x}10$$

(8)

$${\mathrm{Q}}_{\mathrm{tmj}}={\displaystyle \frac{\left[\raisebox{1ex}{$C_{tmj}$}\!\left/ \!\raisebox{-1ex}{$\left(Yt+1-Yt\right)$}\right.\right]}{{{\displaystyle \sum }}_{i=1}^J{C}_{tij}}}x100$$

(9)

$$Vtm={\displaystyle \frac{\raisebox{1ex}{$\left[\left({{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{J}}{\mathrm{C}}_{\mathrm{tmj}}\right)-{\mathrm{C}}_{\mathrm{tmm}}\right]$}\!\left/ \!\raisebox{-1ex}{$\left(Yt+1-Yt\right)$}\right.}{{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{J}}\left[\left({{\displaystyle \sum }}_{i=1}^J{C}_{tij}\right)-C_{tim}\right]}}x100$$

(10)

3. Results and Discussions

3.1. Magnitude and Trend of LUC

The study's results in (Figure 3and Table 2) reveal that although the values differed during the FTIS(1986 -2002) and STIS(2002 - 2020), all land use types experienced considerable changes in spatial extent and percent proportion as a result of the land transition processes covered above 2/3 of the study area from 1986 to 2020. Regarding area-shared proportion from 4140ha of NNP, in all reference time points(1986, 2002 & 2020) WA and FL were in the 2^nd and 5^th position, respectively. In 1^st and 2^nd time points, WL shared the largest proportion (39.24 and 38.25%) while CL had the smallest (1.17 and 2.00%), respectively. During the last time point, WL moved to the 3^rd level when BS took over the largest proportionate size (34.82%) mainly by displacing the land of WL (3426 & 6738 ha), FL (330 & 314ha) and GL (280 & 714ha) at the FTIS and STIS, respectively (Table 3).

Concerning the changes in spatial extent, WL exhibited an increased trend of decline and lost 7043ha (43.53%), of land over the last 34 years; of which 6632ha was in the STIS, while CL expanded by 615 ha(126.44%) with the larger expansion in FTIS (341ha). Similarly, GL downsized by 1965ha over the study period; of which 9204ha was in FTIS. Even though the contraction tended down significantly during the STIS, till, GL was under losing land types with 204ha net loss. Moreover, significant shrinkage(44.60%) with an increasing decline trend was also observed in FL; 658ha in FTIS and 768ha in STIS from its total area of 3198ha in 1986. On the contrary, the largest areal expansion was observed in BS by 8587ha, of which 6229 ha was during STIS. Likewise, WA enlarged by 14.84%; 132ha in FTIS, and 1101ha during STIS. Overall, the study indicates that NNP is under an extensive land transition process profoundly dominated by net gain and progression in bush/shrub, cultivated and water areas at the cost of negative net change in forest, grass and woodland. Researchers, like [9,50] pinpointed that such types of land transition processes not only in biodiversity hotspot sites but elsewhere have considerable implications for biodiversity degradation and disturbance of the natural land transition processes among land types and interventions are needed to stop or reverse these transition processes and ecosystem restoration. Moreover, similar to the findings of this study, sizable agricultural land expansion at the expense of more than 3 million km² of vegetation cover of PAs was reported in USA [16]. Moreover, a continued transition of PA’s forest and woodland to anthropogenic land use in developing countries has also been revealed by several researchers [17,37,39,51]. Similarly, previous studies in Ethiopia‘s PAs found rampant LUC dynamics with significant wildlife habitat and biodiversity loss because of increasing human land uses and degradation. For example, a recent study by [26] in Semien Mountain Park found about 366 and 159% expansion in cultivated and built-up areas at the cost of 31 and 16% forest and grassland decline, respectively from 1984 to 2020.

Likewise, the study by [31] in Gibe Sheleko National Park also reported an increment in bush and shrubland by 51.5% and a decline in forest coverage by 66.8%. Furthermore, consistent with the finding of this study, several recent studies conducted on non-protected areas in Ethiopia [45,52,53] have found a rapid expansion of bush and shrubs by displacing forests, grass, and woodlands. This indicates push and shrub encroachment is becoming a serious and common environmental problem, so urgent measures are imperative.

3.2. Gross Gain, Loss, and Persistence of Land Use Types

The total number of pixels for gross gain, loss and persistence, swap and net change of each land type are displayed in (Table 3). The change analysis reveals from the entire landscape of NNP, 53.65% and 64.79% were under LUC/land transitions/ during FTIS and STIS, respectively. Of the total LUC, the net change accounted for 13.67% in FTIS and 36.73% in STIS and the rest for swap. Similarly, from 1986 to 2020 the net change(50.41%) also was higher than the swap(28.41%); which signifies the prevalent of quantity/net gain and loss/ in land types instead of location (land exchange/ change which further escalate the tendency of advancement in gaining land types and degradation in losing types(forest, grass and woodland). Moreover, with different sizes; all land use types experienced gross gains, losses and persistency (Table 3). During the FTIS and STIS, the gross loss of FL, WL and GL were markedly higher than the loss; the reverse is true for other types. BS was first gaining land type with gross gain of 50.34% (4122.36ha) and 54.19%) (7812.72ha) and loss of 30.27% (1764.81ha) and19.34% (1583.46ha) during FTIS and STIS, respectively; while CL was the second by gross gain of 53.25% in FTIS and 42.21% in STIS. Over the study period, with 57.12% of gross loss, WL was the first losing land type followed by FL(55.78%) and GL(49.00%). Regarding persistency, the diagonal axis of (Table 3) shows that from the total land of the study area, 30312.18ha in FTIS and 27988.74ha in STIS were persistent. Among the land types, the highest percentage of persistent land was observed in WA(99.53%), followed by CL (85.54%) and BS (73.40%) with positive net change values of 1233, 615, and 8587ha, respectively from 1986 to 2020. On the contrary, the lowest persistent land was in WL(42.88%), followed by FL(44.22%) and GL(51%). Similarly, these land types experienced net loss in FTIS and STIS. Additionally, the absolute value of their net change was higher than swap change, except GL in STI (Table 3); demonstrating the exceeding of destructive interventions than the rehabilitation and protection measures in the study area.

Furthermore, the change vulnerability indices were computed and presented (Table S 5). According to [23,26] when the value of gain/persistent and loss/persistence ratios are greater than one, a given land type is more likely susceptible to gain from and loss to other types rather than to persist, respectively. The net change/persistence ratio signifies the gain and loss of a land type is how many times its persistent land [45]. The gain/persistent ratio for BS and CL was > 1 and highest among the land types between 1986 and 2020, as shown in Table S 5. At the same time, their loss/persistence ratio was < 1 and the lower, following WA(≈ 0). These results indicate that in the study area, these land types had the highest tendency to expand their territory by displacing other types instead of remaining as persistent with a very low(0.36 =BS and 0. 17= CL) likelihood of losing to others. On the contrary, the gain/persistence ratio of FL, WL and GL was < 1 while its loss/persistency was above and equal to 1; indicating its exposure to intensive degradation with a very low tendency for restoration/gaining/. Similarly, the net change/persistence ratio demonstrates the net gain of BS and CL was 201 and 148% times its persistence, respectively while the net loss of FL and WL was 101% times its persistence. Overall, the vulnerability indices show almost all the land types tended to enlarge/dwindle rather than remain persistent; which indicates the extensive and intensiveness of the LUC processes. Further, the pressures caused deterioration in grass, forest and woodland; provided the opportunities for advancement of bush and shrubs. This result is consistent with the findings of studies done by [25,26] elsewhere in protected areas. Regarding vulnerability status, similar findings were reported by studies conducted in non-protected areas adjacent to NNP and Northern Ethiopia [23,45]. This indicates protection given to NNP is poor or its exposure to anthropogenic pressures is very high. Thus, quantifying and mapping the human influence index is important for the future management plans of the park.

3.3. Intensities of LUC

3.3.1. Interval Level of Intensity Analysis (ILI)

The ILI results present in (Figure 4); the graphs on the left side show the percentage of total land change and on the right side the change intensity. The intensity bar ends to the right of the uniform line indicating the change is fast; whereas if it ends to the left the change is slow [35]. Accordingly, the LUC undertaken in the study area was relatively quicker during STIS(1.80 %) than FTIS(1.66%). Similarly, the total area subjected to land transitions/gain and loss/ was larger during the STIS(64.79%) than FTIS(53.56%). Although the STIS was longer than FTIS by 2 years, the ILI indicates that more intensive change was observed in this time interval. Thus, the plausible reason for accelerating LUC is the rapidly intensified and diversified socioeconomic pressures to realize the need for the increasing population in and around the park as suggested [6,30]. Supporting this argument, the study by [54] identified human interventions as the main driver for pervasive land change in Nijhum Dwip National Park. During the study time, if the LUC was uniformly distributed, . However, the graphs on the right side in Figure 4 show that the overall change intensities in both time intervals were not equal with the intensity of uniform change (1.73 %) indicating the observed LUC in the study area was not uniform and stationary for ILI. similar results were reported by [8,35,36]. The overall statistical output IA displayed in (Table S 6)

3.3.2. Category Level of Intensity Analysis (CLI)

The graphs on the left side of Figure 5A&B display the size of the annual gain and loss of each land type. In terms of size, BS had the largest annual gain while WL followed by GL had the largest loss in both time intervals. The smallest percentage of annual gain and loss was observed in WA and CL. However, the intensities of changes were not proportional to the size of the changes for most of the land types(see on the right side of Figure 5A&B). For example, the gain intensity of CL was the highest in FTIS and higher in STIS, although the size of its gain(≈ 4%) and total area proportion was the smallest. From this it can be noted that land types with the largest annual gain/loss may not have the most active intensity and land types with the lowest annual gain/loss also may not have the most dormant intensity because the transitions/gain and loss/ in a given land type could be caused by either the proportionate size at initial time or intensity of activities on it [8,35]. To understand the reason/s for the gain/loss/ of a given land type; quantification the intensity of observed changes and comparison with intensity of uniform change is, therefore, necessary and helpful [8].

On the right side of Figure 5A &B, the bar graphs show the intensity of annual gain and loss and the vertical broken lines indicate the intensity of uniform change. If the intensity’s bar of a land type stops to the left of the uniform line, the change is dormant; while it extends to the right, the change is active. As the study indicates WA was the only land type that had dormant gain and loss during both time intervals. However, this land type accounted for relatively the larger portion of the study area. It had also a 0.72% gross loss with 0.32% net gain in FTIS and 0.05% gross loss with a 2.63% net gain during STIS. Therefore, changes in WA were significantly associated with its proportionate size and partially for the occasional increase and decrease of water in Lake Abaya and Chamo rather than the intensiveness of the changes. This implies if the LUC was distributed uniformly among the land use types, WA would have less intensive gain and loss, proportional to its size at the initial points; which is consistent with the argument given by [8].

The results indicate WL was leading land types intensively subjected to the LUC process with the largest share of (21%, 8799.30ha in FTIS and 26%,10742.40ha in STIS) the total change. The annual gain of WL was dormant with a lower intensity value (1.60% in FTIS and 1.21% in STIS) in the study area. On the contrary, the annual loss intensity was intensively active and three times higher in STIS (3.05%) than the FTIS (1.77%). As a result, the land transited to other types was increased from 11.2% (287.80 ha yr^-1) with a net loss of 1.38% (411.84 ha) in FTIS to 21.2% (486.62 ha yr^-1) with net loss of 16.39% (6631.19 ha) in STIS. Throughout the study period, the size of the loss accounted for the highest percentage (56.88%) of the entire loss in the study area. These results signify that the WL’s gain was markedly contributed to its larger proportionate size, whereas, its loss was mainly caused by high-intensity driving pressures; i.e., overutilization by local people for livestock grazing and wood production. Similarly, the study done by [49] reported dormant gain and anthropogenic-driven active loss in woodland. On the contrary, due to distant location from human influences; dormant loss in woodland was reported by [55].

The gain intensity(3.20%, FTIS and 3.03%, STIS) of BS was intensively active and higher than its loss intensity(1.89%, FTIS and 1.07%, STIS). In terms of size, it had a gain of 10.03% (262.65 ha yr^-1) and 19.08% (434.07 ha yr^-1) with a net gain of 6.05% (235.55ha) and 15.42% (6229.26 ha) during the FTIS and STIS, respectively. However, the loss intensity of BS was less active in FTIS and dormant in STIS with 4% (110.30 ha yr^-1) and 3.67% (87.97 ha yr^-1) of land loss, respectively. These results indicate that the active gain of BS was primarily linked with the intensiveness of its gaining processes rather than its larger area proportion and the opposite is true for its loss. Researchers, [50] explained that when other vegetation ecosystems are under pressure from anthropogenic activities and changes in abiotic factors: climate change, and soil erosion, the bush and shrubs have advanced opportunities for reproduction and territorial expansion with less likelihood to be replaced by others. Regarding the less active loss in FTIS in the study area, the plausible reason could be the non-discriminable destruction of natural forests throughout the country during the period of political power transition from Dergue to the People Revolutionary Democratic Front; which is in line with the justification given by several researchers for massive deforestation in protected and other areas of Ethiopia [8,24,26,47,56].

Recent studies such as [25,32,38,39,54] found a substantial and continuous decline in forest and grassland in PAs linked with anthropogenic land expansion and overutilization. Similarly, our findings confirm that during the study time, FL and GL experienced intensively active losses. The net loss in FL was 1.62% (657.81 ha) and 1.86% (768.50 ha), while, in GL 4.22% (1761.30ha) and 0.47% (203.49ha) during FTIS and STIS, respectively. FL’s annual gain declined from 1.5% (39.35ha yr^-1) in FTIS to 0.88% (20.61ha yr^-1) in STIS, but GL’s gain increased from 3% (78.95ha yr^-1) to 3.78% (88.11ha yr^-1). However, the graphs on the right sides of Figure 5A&B show the gain intensity of both land types was dormant in both time intervals; which indicates the annual gain of FL and GL was mainly attributed to their relatively larger area proportion rather than the changes intensiveness. In line with this finding, forest’s dormant gain was reported by [49,57]. Concerning intensity of loss, during the FTIS, GL (2.58%) had the most active; followed by FL (2.52%). However, in the STIS, the loss intensity in GL significantly declined to1.79% and became dormant. The convincing reason is due to the high degradation in the grassland during FTIS, the local community shifted to woodland areas for their livestock grazing and fodder collection which resulted in a reduction of pressures and active regeneration in grassland during STIS [30]. The second important factor could be the resettlement of people outside the park during the early time of STIS. Overall, at CLI, the study displays that, the gain and loss in WA, CL, FL, and WL, and gain in BS and GL were stationary, although not perfectly stationary( Figure 5).

3.3.3. Transition Level of Intensity Analysis (TLI)

In Figure 6 the graphs under columns A and C show the annual area (count of pixels) transited to the gaining land type from other types and under columns B and D show the intensity of gain transition. Similarly, in Figure 7 the graphs under E and G display the annual area transition from losing land type, and under columns F and H show the intensity of loss transitions. Since 1986, the LUC observed in the study area is predominantly characterized by the intensive expansion in BS and LC at the cost of FL, GL, and WL, although vital for biodiversity conservation and ecosystem services including livelihood opportunities for local people. Therefore, our TLI focused on the transitions to BS and CL and from FL, WL and GL.

As displayed in Figure 6 and Table 3, the gain of BS was from all land types, but the largest portion was from WL, 83.95% (214.10ha yr^-1) in FTIS and 86.81% (374.30ha yr^-1) in STIS. Likewise, the transition intensity confirms that the gain of this land more intensively targeted WL and avoided others in both time intervals. Therefore, the transitions to BS were stationary, although not perfect. Moreover, the largest gain of BS from WL could be associated with its proportionate size or strength/intensity/of the transition. However, As shown in the graph almost half of the gain transition intensity bars from WL extends above the uniform lines; signifying the systematically targeted transition to BS from WL was significantly contributed to the intensiveness of the change. Similarly, CL’s gain was transferred from all others, except WA. However. the largest portion (above 80%) was intensively targeted from FL in both time intervals, slightly targeted from WL in STIS and avoided others. Thus, even though not perfect, all the transitions to this land were stationary, except from WL. The main logical reason for FL being intensively targeted by the expansion of CL was its closest location to cultivated and rural settlement areas, particularly in the Southeastern part of the park.

Regarding the loss transition, TLI shows from the total annual area transited from GL (189 and 99ha) to all other land types, the largest portion (79.29% in FTIS & 47.62% in STIS) intensively targeted towards WL. The all transitions from GL were stationary but not perfectly stationary. Similarly, When the area of WL annually declined by (287.84 and 482.62ha) the largest transition (74.39%, FTIS, and 77.50%, STIS) was intensively targeted by BS and avoided by others; which were also stationary but none of them perfectly stationary. Subsequently, when FL shirked down by (80.46 and 63.23ha yr^-1), a significant portion was taken by BS (27.53 and 25.77%, WL (19.68 and 34.83%) and CL (16.88 and 26.44%) during STIS and FTIS, respectively. However, intensively targeted by CL in both time intervals, less intensively targeted WA in STIS, and avoided others; therefore, its transition to WL, GL, and CL were stationary.

3.4. Implications of Land Use Change

A widespread and anthropogenic-oriented LUC occurring in biodiversity hotspot areas has significant adverse implications for natural ecosystems and biodiversity [12,56]. The natural forest and woodland are the major natural ecosystems in accommodating biodiversity and provision of multifold ecosystem services in NNP. However, the study’s results show during the last 34 years 1426ha forest and 7043ha woodland was neatly degraded from 3198ha and 16247ha in 1986; the largest portion contributed to the expansion of bush/shrubland (147.27%) and cultivated (126.44%); which profoundly impoverished the functionality of the forest and woodland ecosystems and resulting in substantial wild animals’ habitat destruction and biodiversity loss accompanied by species extinction ( See Figure S 3). In line with this, researchers, such as [12,39,54,58] found that human-induced LUC and degradation have been the major threat to the drastic loss of habitat and biodiversity in many protected areas across the world. Similar findings were also reported in Ethiopia by [25,26].

Moreover, In NNP, Savanna grassland is another important natural ecosystem, playing crucial in biodiversity conservation, and ecosystem services (economic and ecological). It has been serving as a home and vital ground for forage, drinking water, and reproduction for many birds, herbivorous and omnivorous animals [30]. Particularly, it is known for the highest number of plain Zebra; one of the most threatened wild animals in Ethiopia. However, our study reveals that due to the observed LUC, the grassland declined by 1965ha (27%) of net loss, of which more than 80% was attributed to invasive bush and shrub. Certainly, this vast land transition, particularly towards bush and shrubs considerably daunt the roles of grassland in NNP and comes up with multiple ecological consequences; scarcity of fodders and water for wild animals and decline of biodiversity(see Figure 4 and Figure 5). Supporting this finding, the experimental study carried out by [50] confirmed that the land transition process from grassland to bushland often ended with a wide range and prolonged environmental and economic problems.

Apart from biodiversity conservation, NNP is known for its multitudinous ecological services: regulation of precipitation and hydrological processes, carbon sequestration, and climate monitoring which have global and local importance. In addition, the natural springs and rivers that supply water for local people (urban and rural dwellers) and wild animals, emanate and flow inside the park, particularly in the forested areas. In flowless speaking, continued deforestation for human land uses and benefits has deleterious consequences on the ecological services of the park including disturbance of precipitation cycles and shortage of rainfall, prolonged drought, and scarcity of water supply in the park and surrounding areas. In addition, anthropogenic-driven land conversions and degradation have been considered the major causes of soil erosion and sedimentation [8,59]. Among the land types, grassland is more important for controlling soil erosion and transportation activities [59]. Therefore, the degradation of grassland coupled with forest destruction significantly escalates the rate of soil erosion from terrestrial parts and sedimentation in water areas of NNP, which have also negative consequences on the aquatic ecosystems. Furthermore, the expansion of cultivated, and settlement areas certainly amplifies the degradation of important ecosystems and ecological services in different ways. In line with this, [8] explained that forest and grassland degradation and conversion to cultivated land lead to disturbance in hydrological processes, high runoff, increased soil degradation, and water pollution.

Empirical studies [12,37] reported most of PAs are important sources of livelihood for people living in the adjacent areas, particularly, in developing countries. Likewise, for a long time, NNP has played a paramount role in supporting the socioeconomic life of many rural and urban dwellers in different ways: livestock rearing, crop cultivation, fishing, fuel wood production, and so on. Thus, the widespread degradation of forest, woodland and grassland resulting in long-term socioeconomic repercussions, particularly for farmers living in and around the park profundal depending on natural resources for their survival. Therefore, monitoring the ongoing serious LUC and directing the management practices and human interventions towards sustaining the natural resources not only for biodiversity conservation but also significantly important for harnessing socioeconomic opportunities of the park sustainably. In light of this further studies should be carried out on anthropogenic and natural drivers integrated with spatiotemporal LUC in the park.