Agroclimatic Zone-Based Analysis of Rainfall Variability and Trends in the Wabi Shebele River Basin, Ethiopia

: The amount and annual distribution of rainfall caused a major socioeconomic and environmental problem where rainfed agriculture is predominant. This study assessed the long-term variability and trends of rainfall in the Wabi Shebele River Basin (WSRB), Ethiopia. The basin was discretized into 9 local agroclimatic zones (ACZ) based on annual rainfall and elevation. The coefﬁ-cient of variation (CV) was used to check the variability of rainfall while modiﬁed Mann-Kendall (MK) and Innovative Trend Analysis (ITA) methods were used to detect rainfall trends. For each ACZ, stations with long-term records and less than 10% of missing data were selected for further analysis. The mean annual rainfall in the basin ranges from 227.2 mm to 1047.4 mm. The study revealed most of the ACZs showed a very high variation in Belg/Spring rainfall (CV% > 30) than Kiremt/Summer and annual rainfall. Seasonal and annual rainfall trend analysis revealed that no uniform trend was detected in all ACZs. However, most of ACZs in the arid and semi-arid areas showed a non-signiﬁcant decreasing trend in annual rainfall. From seasonal analysis, Belg and Kiremt rainfall showed relatively decreasing and increasing trends respectively. In comparison, a similar result was observed using MK and ITA methods.


Introduction
Climate change is already affecting every inhabited region across the globe with human influence contributing to many observed changes in weather and climate extremes [1]. In many regions, agricultural production is being adversely affected by rising and variability of temperatures, changes in amounts and frequency of precipitation, the increasing intensity of extreme weather events, rising sea levels, and the salinization of arable land and freshwater [2]. Africa is one of the most vulnerable regions highly affected and to be affected by the impacts of climate change [3]. Due to its low adaptive capacity and high sensitivity to socio-economic systems. Similarly, Ethiopia relies on subsistence rainfed agriculture as the main source of national income; consequently, the importance of the timing and amount of rainfall that occurs cannot be overstated [4].
Besides, Ethiopia generates most of the energy needed from hydropower plants which depend on the amount of water stored in the reservoirs. The government of Ethiopia is currently implementing large-scale investments in hydropower and irrigation projects along the major river basins [5]. Thus, the country is vulnerable to climate variability, and climate change is likely to increase the frequency and magnitude of climate disasters [6][7][8]. These have made the country most susceptible to famine usually caused by drought [7]. As a result, climate variability and trend analysis at appropriate spatial and temporal scales are crucial for understanding and mitigating problems in hydrology and water resources, such as water resource development, environmental protection, and ecological balance [9][10][11].
In the past three decades, climate variability and trend analysis in hydro-meteorological data has been conducted by numerous researchers using different methods, data sources, and spatial and temporal scales. Among the methods, the Mann-Kendall test is the most widely used test for detecting monotonic upward or downward trends in hydrometeorological and environmental data on the assumption that the observations in the time series are independent [12]. Accordingly, checking autocorrelation in the time series became a common practice before any analysis of trend detection is performed using the MK test. This test is mostly supported by Sen's slope estimator for estimating the magnitude of a trend in the time series [13]. Both tests are non-parametric tests that do not require the data set to be normally distributed. A score of studies used this method to estimate temporal trends for different climatic variables [6,7,11,[14][15][16][17][18][19][20][21][22].
Recently, Ref. [23] developed the innovative trend analysis (ITA) method to avoid all the restrictive requirements for the statistical test in the MK trend test. The basis of the approach rests on the fact that, if two-time series are identical to each other, their plot against each other shows scatters of points along 1:1 (45°) line on the Cartesian coordinate system [23,24]. National and international researchers have used the ITA method in different parts of the world [24][25][26][27][28]. The main limitation of the ITA method is that it does not allow the determination of the differences between each point and the 1:1 line is statistically significant [28].
Even though climate variability and change have location-specific impacts especially in Ethiopia due to high spatial and temporal variation of rainfall and topography effect, few studies have been carried out regarding rainfall variabilities and trends in the Wabi Shebele River basin. Ref. [26] assessed the spatiotemporal distribution and variability of rainfall in seasonal and annual time series in the Upper Wabi Shebele river basin using station data in MK and ITA tests. In the northern district of WSRB, Ref. [18] investigated the spatiotemporal variability and trends of rainfall and its association with Pacific Ocean SST using satellite-based rainfall data and the Mann-Kendall trend test. Because of the uneven distribution of rainfall recording stations in the basin, a study based on station data may not be sufficient for planning and decision-making. Therefore, this study investigated variations and trends in maximum daily, rainy seasons and annual rainfall in the WSRB, at the local agroclimatic zone, using modified MK and ITA methods. In general, a study based on local agroclimatic zonation overcomes the data scarcity challenge and simplifies sharing the findings with the local community.

Study Area
Wabi Shebele River Basin (WSRB) is one of the largest river basins in Ethiopia in terms of area coverage. The basin, which occupies a total area of 189,655 km 2 comprising nearly 17% of the country's total land area, lies between 4.91°-9.59°N latitude and 38.69°-45.33°E longitude Figure 1. The elevations of the basin vary from 179 m.a.s.l along the Ethio-Somalia border in the south to more than 4188 m.a.s.l in the northwest part of the Bale mountains. More than 80% of the area has an elevation of less than 1500 m and was considered low land (arid and semi-arid). The upstream of the basin is characterized by steep valley slopes while the lowland has gentle slopes. Annual average rainfall over the basin varies from 227 mm in the southeast to 1240 mm in the northwest. Similarly, the minimum and maximum average temperatures are 6.29°C in the northwest and 35.98°C in the southeast respectively. Due to bimodal rainfall distribution in the basin, there are two rainy seasons, the main rainy season (Kiremt/Summer) spans from June to September and the short rainy season (Belg/Spring) extends from March to May. The dominant soil types are Gypsisols, Calcisols, and Leptosols covering 68% of the basin. A large percentage of the population in the highlands depends on crop cultivation while the lowlanders, in general, are pastoralists.

Data Source and Preparation
For this study, the observed daily rainfall data of 43 stations located in the basin and 21 of its surrounding for the period from 1987 to 2016 was obtained from the National Meteorology Agency (NMA).
For each observation station, high-resolution (4 km) satellite-gauge merged rainfall data was collected from NMA. The missing data were filled in using satellite-gauge merged rainfall data because of the unpredictability of the climate in the tropical zone and the lack of reference stations to use statistical methods. For the analysis, the daily rainfall records were added to seasonal and annual rainfall. Also, daily maximum rainfall was extracted from daily rainfall data.

Agroclimatic Zonation
There are 17 local agroclimatic zones (ACZs) in Ethiopia based on the annual average rainfall (AARF) amount and elevation of the area [40]. Similarly, DEM of 30 m grid size and annual average rainfall were used for agroclimatic zonation of WSRB. No data and sink values of the DEM were filled using the QGIS tool. The annual average station rainfall data was converted to a spatial rainfall map using an interpolation algorithm in the QGIS environments. The elevation (filled DEM) and rainfall spatial data were reclassified based on the value range given in Table 1. Then, the reclassified maps merged to create new classes representing the unique agroclimatic zone. The basin consists of 9 out of 17 ACZs in Ethiopia. The spatial coverage and characteristic of each ACZ are presented in Figure 2 and Table 1 respectively. Based on annual rainfall amount, ACZ1, ACZ2, ACZ3, and ACZ4 are categorized as dry regions covering 81.21% of the basin whereas ACZ5, ACZ6, and ACZ7 are categorized as humid regions covering only 18.79% of the basin. This makes the basin significantly vulnerable to climatic-related issues, especially drought. On other hand, ACZs are grouped as Low land (ACZ1, ACZ2, and ACZ5), Mid-Highland (ACZ3 and ACZ6), and Highland (ACZ4, ACZ8, ACZ7, ACZ9) based on the elevation.  After local agroclimatic zonation, one station from each ACZ was selected based on the long-term availability of records and less than 10% of missing records. Dry Wurch and humid Wurch agroclimatic zone were excluded from the analysis due to the unavailability of a rainfall recording station. Besides, these zones cover less than 0.6% of the basin. Thus, rainfall data representing only 7 ACZs were considered for further analysis. The location and spatial distribution of selected stations are shown in Figure 1. The details of the meteorological stations used in this study are presented in Table 1. The reliability and quality of the data to be used in the analysis should be checked statistically. In this study, Standard Normal Homogeneity Test (SNHT) [41] and Pettitt tests [42] were used for detecting the inhomogeneity of seasonal and annual rainfall data. Due to uneven distribution and scarcity of stations in the study area, only the station found inhomogeneous in each test was eliminated from further analysis. The result of the SNHT and Pettitt test is summarized in Table 2.

Rainfall Seasonal and Annual Variability
For each ACZ, descriptive statistics for the long-term rainfall data were computed using standard statistical procedures using Microsoft Excel, including mean, minimum, maximum, variance, standard deviation (SD), and Coefficient of Variation (CV). The CV is a commonly used statistical measure to examine the inter-annual variation of rainfall. In this study also, the inter-annual variability of daily maximum, rainy seasons, and annual rainfall was evaluated using CV. Greater values of CV indicate larger variability, and vice versa. The CV (%) value for each series can be computed as follows: where σ is the standard deviation and µ is the mean of time series data. In the literature, the degree of variability of rainfall based on CV classified as low (CV < 20), moderate (20 < CV < 30), and high (CV > 30) [6,15,43].

Rainfall Trend Analysis
Two non-parametric (Mann-Kendall and Sen's slope estimator) and the ITA methods were used to assess the rainfall temporal trend in the Wabi Shebele River basin. The detail of these methods presented in consecutive subsections.

Modified Mann-Kendall (MK) test
The Mann-Kendall test [12] is perhaps the most widely used non-parametric test for detecting trends in hydro-meteorological and environmental data. It detects monotonic trends where data is either consistently increasing or decreasing over time. MK test does not assume the data to be normally distributed and it is flexible to outliers. While the test requires time series data to be serially independent. The test assumes a null hypothesis, H O , of no trend and an alternate hypothesis, H A , of increasing or decreasing monotonic trend. The mathematical equations for Mann-Kendall statistics S, the variance associated with S V(S), and standardized test statistics Z are as follows [12,14,44].
In these equations, x i and x j are the time-series observations in chronological order, n is the length of the time series, t p is the number of ties for the pth value, and q is the number of tied values. A positive (negative) value of Z indicates that the data tend to increase (decrease) with time. The strength of the trend is proportional to the magnitude of the Mann-Kendall Statistic. The null hypothesis of no trend is rejected if the absolute value of Z is greater than 1.96, and 2.58 at 5%, and 1% significance levels respectively.

Sen's Slope Estimates
Sen developed a more robust and non-parametric procedure for estimating the magnitude of a trend in the time series [13]. The slope of data value pairs was calculated to get Sen's slope Q i by the relationship: The N values of Q i are ranked from smallest to largest and the median of slope or Sen's slope estimator (Q med ) is computed as: if N is odd.
, if N is even.
Positive/negative values of Q med indicate an increasing/decreasing trend, respectively, while its value indicates the magnitude of the trend.
Simulation experiments demonstrated that the existence of serial correlation alters the variance of the estimate of the Mann-Kendall (MK) statistic; and the presence of a trend alters the estimate of the magnitude of serial correlation [45]. Several pre-whitening procedures have been proposed by different authors to remove the influence of the autocorrelation including variance correction [46], effective sample size [47], and trend-free pre-whitening (TFPW) [44,45]. In this study, Ref. [46] procedure was applied to check and eliminate the autocorrelation before performing trend analysis. Serial autocorrelation was checked with an autocorrelation function in the modifiedmk R package [48]. Numerous researchers have also used this approach to eliminate serial correlation in time series data [7,14,17].

Innovative Trend Analysis (ITA)
Recently, Ref. [23] developed the ITA method to avoid all the restrictive requirements for a statistical trend test such as the MK test. In this trend analysis method, graphs are used as templates. The basis of the approach rests on the fact that, if two time series are identical to each other (no trend), their plot against each other shows scatters of points along 1:1 (45 • ) line on the Cartesian coordinate system [23,24]. In this method, the time series is divided into two equal parts, which are separately sorted in ascending order. Then, the first and the second half of the time series are plotted on the X-axis and the Y-axis, respectively. If data points are located in the lower triangular area of the no-trend line, there is a decreasing trend in the time series, and vice versa [23,27]. In this study, the ITA method was applied to the daily maximum, rainy seasons, and annual rainfall from 1987 to 2016. The rainfall data were divided into two sub-series 1987-2001 (the first half) and 2002-2016 (the second half). R 'Trendchange' package was used to investigate the trend developed by [49]. The rainfall data points are further divided into "low", "medium", and "high" categories based on the position in the ITA graph. Strong trends of extreme ("high" and "low") rainfall are associated with floods and drought (Alifujiang et al., 2020). Besides, to illustrate graphically the distance of the data points from the 1:1 line, in this paper, a confidence band (0.05 and −0.05) has been plotted. This is to help the reader to better realize the differences between the data points and the trendless line, without any statistical implications [24]. A flowchart was presented showing major activities in the methodology Figure 3.

Rainfall Distribution and Variability
In WSRB, the annual rainfall distribution is bimodal with two rainy seasons Figure 4. The first one is the short rainy season (Spring) which starts from March to May months is locally called 'Belg'. The second one is the long rainy season (Summer) which spans from June to September locally called 'Kiremt'. However, the amount of Belg rainfall was found higher than Kiremt rainfall in ACZ1-3. In this document Spring or Belg and Summer or Kiremt are used interchangeably. The temporal annual rainfall pattern at each agroclimatic zone is shown in Figure 5. A clear decrease in rainfall was observed in ACZ3 in the analysis period.  For each ACZ, the rainfall statistical analysis results are presented in Table 3. The minimum mean annual rainfall recorded was 227.2 mm in ACZ1 whereas the maximum mean annual rainfall recorded was 1047.6 mm in ACZ 6. Seasonally, a high amount of rainfall was observed in Summer except in ACZ1 and ACZ2. Based on the CV result, the annual rainfall CV ranges between 16.7% and 46.7%, these shows high inter-annual variability of annual rainfall in the basin. Seasonally, Belg/Spring rainfall showed the highest inter-annual variability in all ACZs. But, the maximum inter-annual variability was detected in Kiremt rainfall in ACZ1 (87.4%) and ACZ2 (68.9%). In comparison, annual rainfall showed a lower value of CV in all ACZs.

Trend Analysis
To detect the trend in a time series, the statistical tests assume the subsequent data in the series to be independent. Serial autocorrelation was checked for the data series of daily maximum, Belg (short rainy season), Kiremt (heavy rain season), and annual rainfall. For a series with 30 observations, the critical value is ±0.358 at the 95% confidence interval. Autocorrelation analysis indicated only three series were affected by autocorrelation two in daily maximum rainfall and one in annual rainfall Figure 6. Both Belg and Kiremt rainfall were found serially independent. Thus, pre-whitening of the data series was conducted before using it for trend analysis.

MK and Sen's Slope Test
Results of MK and Sen's slope analysis on the daily maximum, spring, summer, and annual rainfall for each ACZ in the basin are shown in Table 4. From trend analysis result, daily maximum rainfall has shown statistically significant increase in ACZ1, ACZ5, and ACZ6, an insignificant increase in ACZ4 and ACZ7, and an insignificant decrease in ACZ2 and ACZ3.
As depicted in Table 4, the declining trend of Belg/Spring rainfall is not statistically significant in all ACZs except ACZ4 where rainfall is statistically significant decreasing by −4.11 mm/year. In Kiremt/Summer rainfall, only one ACZ showed a significantly decreasing trend at a 5% significance level, on the contrary, the remaining ACZs showed a statistically insignificant increase. The analysis of annual rainfall showed a statistically insignificant decrease (ACZ2 and ACZ4) or increase (ACZ1, ACZ5, ACZ6, and ACZ7) in the trend. ACZ3 which is a dry middle elevation zone showed a nonsignificant decrease in daily maximum and Belg/Spring rainfall, a significant decreasing trend in Kiremt/Summer, and annual average rainfall.  Note: *, ** statistically significant at 0.05 and 0.01 significance level.

ITA Method
Figures 7-10 shows the results of the ITA method applied to the seven ACZs of the research area on daily maximum, spring, summer, and annual rainfall. ITA result of daily maximum rainfall Figure 7, the rainfall indicated an increasing trend in all rainfall categories in ACZ1 and ACZ5. Contrary, maximum daily rainfall showed a strong decreasing trend in ACZ2. All of the other ACZs showed a quite distinct pattern in each rainfall category. All humid ACZs showed that most rainfall data points fall above the trendless line, implying an overall increasing trend. However, dry ACZs indicated an increase in ACZ1 and decreasing in ACZ2 and ACZ4 trends. Whereas, ACZ3 showed a decreasing trend in low and an increasing trend in high rainfall categories.  At the seasonal scale Figures 8 and 9, a decreasing trend was observed in ACZ1 and ACZ4 in all categories of Belg rainfall Figure 8. No clear trend was detected in ACZ2, ACZ5, and ACZ7. Negative patterns were detected in ACZ3, for the medium and the high values. In ACZ6, negative patterns were detected, but only for the low and the high rainfall values. In Summer, a strong negative trend was observed in ACZ3 with all the values positioned below the 0.05 relative band limit. In ACZ1 the data points show an increasing trend in medium and high rainfall categories whereas decreasing in low rainfall categories Figure 9. A contrasting pattern was found in ACZ2 and ACZ5. In ACZ2 medium rainfall categories showed an increasing trend whereas low and high rainfall categories showed a decreasing pattern. The opposite pattern was detected in ACZ5. In ACZ6 and ACZ7, an increasing trend was observed in the medium and high rainfall categories. Overall, an increasing trend was detected in dominant ACZs WSRB in Summer rainfall. However, a negative pattern was observed in Belg rainfall as a result of the majority of the data points located below no trend line.  On an annual scale Figure 10, a prevailing negative trend for all categories of rainfall (low, medium, and high) has been detected in ACZ3 and ACZ4 with almost all values exceeding the 0.05 relative band limit in ACZ3. In ACZ2 and ACZ6, medium rainfall categories showed decreasing trend whereas low and high rainfall categories showed an increasing trend. In ACZ5, the data points show an increasing trend in medium and high rainfall categories whereas decreasing in low rainfall categories. Relatively, ACZ7 showed unclear trend with the value concentrated within the 0.05 band limit.

Rainfall Variability in the Basin
From the monthly distribution of annual rainfall, Kiremt/Summer rainfall is the main rainfall season for the area with an elevation greater than 1500 m.a.s.l. On the contrary, Belg/Spring rainfall is the main rainfall season in the area less than 1500 m.a.s.l. Comparing the value of CV, ACZ1 showed a higher seasonal and annual variability of rainfall. Relatively, most of the ACZs in higher elevations showed lower inter-annual variation (CV < 20%) whereas ACZs in lower elevations showed moderate to high variation (CV > 20%) of annual rainfall. Higher inter-annual variability of maximum daily and Belg rainfall was observed in all ACZs. The results of this study are generally consistent with others research findings that were conducted in different parts of the country. These research documents reported that rainfall in Belg/Spring exhibits higher inter-annual variability than Kiremt/Summer rainfall [25,29,36,43]. Similarly, previous studies found that the areas with high rainfall amounts exhibit a low coefficient of variations at different time series [4,15,26,31,36,43].

Rainfall Trend
From MK analysis, in daily maximum rainfall a statistically significant increasing trend was detected in ACZ1, ACZ5, and ACZ6 with increases of 1.02, 0.72, and 0.52 mm/year, respectively. In addition, seasonal and annual trend tests showed a mix of increasing and decreasing trends in different ACZs. Analysis of Belg rainfall indicated that all of the ACZ in the study basin showed a decreasing trend over the last 30 years (1987-2016). The maximum rate of change in Belg/Spring rainfall was −2.22, −2.38, and −4.11 mm/year in ACZ7, ACZ3, and ACZ4, respectively. The Kiremt rainfall , however, showed a non-significant increasing trend for all ACZs except ACZ3 showed a significantly decreasing trend at a 0.05 significance level. The analysis also revealed that three ACZs (ACZ2-4) depicted decreasing trends of annual rainfall, besides ACZ3, which is statistically significant. The maximum annual decreasing rate was observed in ACZ2, ACZ4, and ACZ3 at −0.59, −1.96, and −12.3 mm respectively between 1987 and 2016. Among the highland ACZs, ACZ6 showed the maximum increase in annual rainfall with 3.73 mm/year.
In the WSRB, non-uniform trend was detected among ACZs that were considered at different temporal scales using the MK method. The finding agrees with previous studies conducted in the different river basins of the country [25,26,30,[36][37][38]43]. In Lake Tana Sub-basin, Ref. [6] revealed a general decreasing trend in annual rainfall. Ref. [15] reported dominantly decreasing trends in small and main rain seasons and annual rainfall in the Choke Mountain Watersheds (1981-2016). Ref. [36] found an increasing trend in the Kiremt season and annual rainfall in most of the stations investigated, however, the Belg season rainfall showed a non-significant declining trend in four of the five representative stations in Modjo River Watershed, Awash River Basin of Ethiopia. Similarly, Ref. [19] reported a decreasing trend in all stations located in the semi-arid areas of Western Tigray, Ethiopia. Ref. [26] found mixed results of increasing and decreasing seasonal and annual rainfall in the Upper Wabe Shebelle River Basin, Ethiopia. On the contrary, Ref. [9] found annual and seasonal rainfall totals showed increasing trends in the Alwero watershed, in western Ethiopia.
In this study, the ITA method was also applied to investigate the existence of a statistical trend in maximum daily, rainy season, and annual rainfall. The advantage of ITA over the MK method is that detailed trend evaluation can be made by dividing the rainfall into three categories that are "low", "medium", and "high". Uniform trends were not detected in rainfall in all ACZs. Except for ACZ2, all ACZs showed more points above the no-trend line implying an overall increasing trend of daily maximum rainfall. On seasonal rainfall, a decreasing trend was detected in ACZs based on the percentage of data points located below the no-trend line. On other hand, Kiremt rainfall showed a generally increasing trend in all ACZs except for ACZ3. Similarly, ACZ3 showed decreasing trend in annual rainfall in low, middle, and high categories. In the annual time series, the data points are concentrated within the ±0.05 band line. Even though a limited documents were found on trend analysis using the ITA method, the finding of this study is consistence with previous studies that showed a non-uniform trend result [24,26,28].

MK and ITA Method Comparison
The results of the ITA have been compared with the ones obtained through the MK method. In this study, all statistically significant increasing or decreasing trends at 0.05% of the level of significance using the MK method showed a similar trend in all categories of rainfall using the ITA method. On the contrary, even though most of the data points fall above the no-trend line in all categories, there is no significant increasing or decreasing trend detected in the time series data using the MK method. In general, a similar result was observed using the MK and ITA method. But, an in-depth assessment of rainfall trends in the low, middle, and high categories was only possible using the ITA method. Because of the higher variability of rainfall in the basin, the ITA method is more complex to decide the existence of a trend at a different period than the statistical method. The data points are very scattered in the plot. Even though ITA has the advantage of estimating hidden trends at three categories of the data point, in this study, it was found that defining the boundary of low, middle, and high data values is not as simple as stated in the initial work [23].

Conclusions
This study investigated variability and trend in maximum daily, rainy seasons, and annual rainfall in the Wabi Shebele River basin over 30 years (1987-2016) using MK and ITA methods. From the agroclimatic zonation of the WSRB, arid and semi-arid climates cover 81% of the area whereas only 19% are considered humid climates. The distribution of annual rainfall in the basin is bimodal. In this study, lowland and low rainfall areas tend to experience high rainfall variability than highland and higher annual rainfall areas. The findings indicate that, according to the dataset and methods used in this study, annual rainfall in WSRB shows a mixed upward or downward trend with respect to time. The analysis of trends in seasonal rainfall for different parts of the basin found that while the highly variable Belg rain displayed any insignificant decrease over the study period, the Kiremt rain was insignificant increasing. In line with this, agricultural practices in Belg/Spring season will be affected if a similar trend continues particularly in lowland areas where Belg is the main rainy season. In general, the amount of annual rainfall shows a decreasing trend with higher interannual variability. Coupled with climate change, a decreasing trend of rainfall severely impacts the livelihoods of the rural communities relying on subsistence rainfed agriculture.
Finally, the findings of this study provide valuable information on the distribution and variability of rainfall in the WSRB essential for planning and designing sustainable water management strategies and minimizing the impact of droughts and floods. Further study should be conducted to predict the future pattern of the rainfall in the basin using an appropriate climate model under different projection scenarios. Data scarcity was the major limitation of this study in that each agroclimatic zone is represented by one meteorological station. Thus, maximum effort should be made to establish new meteorological stations, especially in the middle and lower part of the basin.
Author Contributions: A.T.T. collected and analyzed data, and wrote the original manuscript; A.M. major supervisor and A.K.K. co-supervisor supervised and instructed the analyses of the results and modified the manuscript. This research article is part of PhD research for the first author. All authors have read and agreed to the published version of the manuscript.
Funding: This PhD study was funded by the German Academic Exchange Service (DAAD) (EECBP HomeGrown PhD Scholarship Programme) and the Universität der Bundeswehr München (Scholarship and Support Program for Foreign Students and Doctoral Candidates (STIBETIII)).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Data will be available upon request to the authors.