# Climate variability and trends in the Endorheic Lake Hayk basin: implications for Lake Hayk water level changes in the lake basin, Ethiopia – Environmental Systems Research

May 12, 2022

### Description of the study area

Ethiopia is situated in Eastern Africa, between 3°  and 15°  latitude and 33°  and 48°  longitude (Horn of Africa). The Lake Hayk basin is a naturally closed (endorheic) drainage that belongs to one of Ethiopia’s most vulnerable zones to climate change and variability. Its areal extent is within 39.68°E to 39.81°E, 11.24°N to 11.39°N and its area coverage including the lake water surface area of 2156.76 ha is 8592.68 ha (Fig. 1).

The Lake Hayk basin is under a subhumid tropical climate with bimodal precipitation regimes (kiremt and belg). The lake basin received a mean annual precipitation of 1192.31 mm; the mean annual surface temperature was 17.58 °C.

### Data sources

The hydroclimate variability/trend in the Endorheic Lake Hayk basin was analyzed using mean monthly historical datasets of precipitation, mean temperature (Tmean) and Lake Hayk’s Water Level (LWL) from 1986 to 2015. Precipitation data from only one Hayk meteorological station (11.31°N, 39.68°E; 1984 m amsl) outside the study area (Fig. 1) has been used to observe the Lake Hayk water level response to climate change/variability. Evidently, the degree to which precipitation amounts vary across an area is an important characteristic of the climate of an area that affects hydrology of lakes. Keeping this in mind, we believe that the Hayk meteorological station is the sole relevant and appropriate station from which precipitation data is enough to evaluate the impacts of climate on Lake Hayk water levels. This is from two perspectives. The first is that, despite being outside the basin, it is very close to Lake Hayk (less than 6 km away), even closer than several points within the basin. Its proximity to Lake Hayk allows it to collect data that is nearly identical to what the lake and its surroundings receive. The other reason is that the lake basin is small (85.93 km2), resulting in a density of meteorological stations in the lake basin of 85.93 km2 per station, which adequately represents the Lake Hayk basin according to the World Meteorological Organization (WMO) recommendation of 300–1000 km2 per station in Temperate, Mediterranean and Tropical zones (Dingman 2002). In light of these considerations, the authors used only the Hayk meteorological station rather than interpolating data from other stations, which could compromise data quality. The Ethiopian National Meteorological Agency has provided us with data for the lake basin’s mean monthly precipitation (1986–2015) and temperature (1994–2015). In addition, due to a lack of station data, the reanalysis temperature products (RTPs) of the same station for the years 1986–2015 were retrieved from the climate explorer (https://climexp.knmi.nl/) portal.

Furthermore, the LWL data of Lake Hayk are measured using the water level measuring gauge situated on the southwest shore of Lake Hayk (Fig. 1). The lake average daily water level time series from 1986 to 2015 were provided by Ethiopia’s Ministry of Water Resources. However, the LWL time series was riddled with severely missed data. Daily data with more than 10% missing values must be excluded from the analysis (Seleshi and Zanke 2004). Full daily LWL data were available for the years 1999–2005 and 2011–2015. Therefore, we fused the water level observations from these periods with remote sensed water extents to bridge the gap between 2005 and 2011. Cloud free (clouds cover ≤ 10%) Landsat 5 Thematic Mapper (TM) images for 2009–2011 and Landsat 7 Enhanced Thematic Mapper Plus (ETM +) images for 2005 and 2008 years were retrieved from the Earth Explorer (http://earthexplorer.usgs.gov/) archiving system. To ensure greater accuracy of interpretation, all Landsat images were downloaded for months of the dry (bega) season of the year.

### Data analysis

This study combined hydroclimate data from gauging station, gridded data (reanalysis) and remotely sensed satellite data to analyze climate variability/change and its implications on changes in the water level of endorheic Hayk Lake at the local level, using statistical approaches with the integration of remote sensing and geographic information system. The general methodology of the study is depicted schematically in Fig. 2.

### Evaluating the reanalysis temperature data

Reanalysis products are thought to be useful in situations when meteorological stations are insufficient and unevenly dispersed, as well as in cases where missing records and short period observations exist (Dee et al. 2011). Climate reanalyses such as the European Center for Medium-Range Weather Forecasts (ECMWF) ReAnalysis 5th generation (ERA5) (Hersbach et al. 2020), the Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA-2) (Reinecker et al. 2011) and the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) (Kalnay et al. 1996) are currently in use. Due to the scarcity of historic station temperature records in the Endorheic Lake Hayk basin, we relied on reanalysis products (historic estimates produced by combining a numerical weather prediction model with observational data from satellites and ground observations) as the best alternative solutions, but their performance evaluation should no longer be overlooked. Therefore, the ERA5, MERRA-2 and NCEP/NCAR reanalysis temperature products (RTPs) in the Lake Hayk basin were quantitatively evaluated against ground station temperature data for the 1994–2015 time series on annual and seasonal scales using coefficient of determination (R2), root mean square error (RMSE) and relative bias (Alemseged and Tom 2015; Nkiaka et al. 2017).

$$R^{2} , = ,,,left[ {{{left( {sumlimits_{t = 1}^{n} {left( {T_{r} – overline{{T_{r} }} } right),,left( {T_{s} – overline{{T_{s} }} } right)} } right)} mathord{left/ {vphantom {{left( {sumlimits_{t = 1}^{n} {left( {T_{r} – overline{{T_{r} }} } right),,left( {T_{s} – overline{{T_{s} }} } right)} } right)} {left( {sqrt {sumlimits_{t = 1}^{n} {left( {T_{r} – overline{{T_{r} }} } right),,sumlimits_{t = 1}^{n} {left( {T_{s} – overline{{T_{s} }} } right)} } } } right)}}} right. kern-nulldelimiterspace} {left( {sqrt {sumlimits_{t = 1}^{n} {left( {T_{r} – overline{{T_{r} }} } right),,sumlimits_{t = 1}^{n} {left( {T_{s} – overline{{T_{s} }} } right)} } } } right)}}} right]^{2}$$

(1)

$$RMSE,, = ,,,sqrt {{{sumlimits_{t = 1}^{n} {left( {T_{r} – T_{s} } right)^{2} } } mathord{left/ {vphantom {{sumlimits_{t = 1}^{n} {left( {T_{r} – T_{s} } right)^{2} } } n}} right. kern-nulldelimiterspace} n}}$$

(2)

$$Bias,, = ,,,left( {{{sumlimits_{t = 1}^{n} {left( {T_{r} – T_{s} } right)} } mathord{left/ {vphantom {{sumlimits_{t = 1}^{n} {left( {T_{r} – T_{s} } right)} } {sumlimits_{t = 1}^{n} {T_{s} } }}} right. kern-nulldelimiterspace} {sumlimits_{t = 1}^{n} {T_{s} } }}} right),,times,100%$$

(3)

where Tr and Ts denote reanalysis and ground station temperature records respectively and n is the length of data. R2 varies within 0 ≤ R2 ≤ 1; R2 = 0 reveals no correlation and R2 = 1 indicates perfect correlation between the reanalysis product and station temperature record. Bias detects a systematic error in temperature values. Zero bias indicates absence of systematic error, whereas negative/positive biases reveal respectively underestimation and overestimation of values (Alemseged and Tom 2015). The RMSE measures residual dispersion (estimation errors) around the best fitting line. RMSE near zero would be a better fit to the data.

### Variability and trends analysis of hydroclimate time series

Various statistical approaches were used to examine the variability/trend in the hydroclimate time series of the Endorheic Lake Hayk basin from 1986 to 2015. The coefficient of variability (CV), the standardized rainfall anomaly (SRA) and the precipitation concentration index (PCI) were employed to study variability of the data. The Modified Mann Kendall (MK) trend test method and the Sen Slope estimator were applied to analyze the significance and magnitude of trend respectively using XLSTAT software. The CV value represents the level of variability in the dataset and is defined as the standard deviation (SD) to mean value (μ) ratio (Hare 2003).

$$CV, = ,,,left( {{{SD} mathord{left/ {vphantom {{SD} mu }} right. kern-nulldelimiterspace} mu }} right),times100$$

(4)

Hare (2003) characterizes variability as being less for CV values less than 20, moderate for CV values between 20 and 30 and high for CV greater than 30. PCI examines the heterogeneity of mean monthly precipitation data. For Pi is the ith month precipitation magnitude, Oliver (1980) defines PCI as follows:

$$PCI,, = ,,,,,left[ {{{left( {sumlimits_{i, = 1}^{12} {P_{i}^{2} } } right)} mathord{left/ {vphantom {{left( {sumlimits_{i, = 1}^{12} {P_{i}^{2} } } right)} {left( {sumlimits_{i, = 1}^{12} {P_{i} } } right)^{2} }}} right. kern-nulldelimiterspace} {left( {sumlimits_{i, = 1}^{12} {P_{i} } } right)^{2} }}} right],,times,100$$

(5)

Precipitation concentration can be identified as low concentration (uniform distribution of precipitation) for PCI values lower than 10, high for values from 11 to 20 and very high for values above 21 (Oliver 1980). SRA offers insights on the occurrence and severity of drought periods. Standardized rainfall and temperature anomalies have no units. They are dimensionless. To get more information about the magnitude of the anomalies, standardized anomalies are calculated by dividing anomalies (deviations of mean value from each observation) by the standard deviation to remove influences of dispersion. Therefore, for Pt is annual precipitation at a year of interest t and Pm is the mean annual precipitation value during the study period, SRA can be estimated according to Agnew and Chappell (1999):

$$SRA, = ,left( {P_{t} – P_{m} } right),/,SD,,,,$$

(6)

Then, severity of drought can be categorized as extreme drought (SRA <  − 1.65), severe drought (− 1.28 > SRA >  − 1.65), moderate drought (− 0.84 > SRA >  − 1.28) and no drought (SRA >  − 0.84) (Agnew and Chappell 1999).

The modified Mann Kendall (MK) trend test was used to examine the monotonic trends of hydro climatic time series in the endorheic Lake Hayk basin from 1986 to 2015 at a significance level of 5% on a monthly, annual and seasonal basis. It was chosen because it is a rank-based (i.e., less affected by low-quality hydro climatic data-data with missing values and/or outliers) nonparametric (i.e., less sensitive to skewed datasets-applies for all distributions) method (Hirsch and Slack 1984). The MK tests the null hypothesis (H0) assuming no trend against the alternative hypothesis of monotonic trend (Ha) using either the S statistics (n < 10) or the standardized normal Z statistics (n ≥ 10) (Hirsch and Slack 1984; Yue et al. 2002). The MK test S statistic is calculated using the following equations (Eqs. 7 and 8) as:

$$S, = ,sumlimits_{i = 1}^{n = 1} {sumlimits_{j = i + 1}^{n} {{text{sgn}} left( {x_{j} – x_{i} } right)} } ,,,,,,,$$

(7)

$$,,{text{sgn}} left( {x_{j} – x_{i} } right),, = ,,left{ begin{gathered} + 1,,,,if,,theta ,,, > ,0 hfill \ 0,,,,,,,if,theta , = ,0 hfill \ – 1,,,,,if,theta , < ,0 hfill \ end{gathered} right.,,,,,,,,$$

(8)

where n is the data size and xi and xj are the data values at times i and j respectively, i = 1, 2,…, n−1 and j = i + 1, i + 2…, n. Every value in the chronologically ordered time series is compared to every value preceding it, yielding a total of n (n – 1)/2 pairs of data. The total of all rises and falls result in the ultimate value of S (Yue et al. 2002). S values can be positive to show rising trends or negative to indicate falling trends.

When n ≥ 10, the S statistic is assumed to have a normal distribution, with the mean becoming zero and the variance computed using the following equation (Eq. 9) (Kendall 1975):

$$Vleft( S right), = ,frac{1}{18}left[ {nleft( {n – 1} right),left( {2n + 5} right), – ,sumlimits_{t = 1}^{m} {t_{i} left( {t_{i} – 1} right),left( {2t_{i} – 5} right)} } right],,$$

(9)

where V (S) is the variance of S statistics, m denotes the size of tied groups (groups with similar values) and ti represents the size of data points in the ith tied group. Then, the Z test statistics can be calculated from the known values of S and V(S) using the following equation (Eq. 10):

$$Z, = ,,left{ begin{gathered} {{left( {s – 1} right)} mathord{left/ {vphantom {{left( {s – 1} right)} {sqrt {v,,left( s right)} ,,}}} right. kern-nulldelimiterspace} {sqrt {v,,left( s right)} ,,}},,,,,,if,,S > 0 hfill \ 0,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,if,,,S = 0 hfill \ {{left( {s + 1} right)} mathord{left/ {vphantom {{left( {s + 1} right)} {sqrt {v,,left( s right)} }}} right. kern-nulldelimiterspace} {sqrt {v,,left( s right)} }},,,,,,,,,if,,S < 0 hfill \ end{gathered} right.,,,,,,,,,,,,$$

(10)

The resulting Z values indicate the direction of a trend (values can be positive to show rising trends or negative to indicate falling trends). Furthermore, the Z statistic is used to measure significance of a trend. When testing for a trend (2 tailed) at significance α, H0 is rejected if (left| Z right|,) equals or exceeds its critical value (,left( {left| Z right|,,, ge ,,Z_{alpha /2} } right)). For instance, if the 5% significance level is used, H0 is rejected when (,left( {left| Z right|,,, ge ,,1.96} right)) or P ≤ 0.05, indicating that a trend exists (A time series has a trend when it is significantly correlated with time). The letter P symbolizes the probability of risk to reject/accept H0 while it is true.

Prior to trend testing, it is critical in time series analysis to examine autocorrelation or serial correlation, which is frequently overlooked in many trend detection studies. To account for the effect of autocorrelation, Hamed and Rao (1998) propose a modified Mann–Kendall test rather than the original MK test, which should be used only on datasets with no seasonality or significant autocorrelations. This is because the presence of significant autocorrelation in a dataset can alter the variance of the original MK test (the existence of positive autocorrelation will lower the actual value of V (S) and vice versa). Hence, when data exhibit autocorrelation, the modified MK test calculates the modified variance using the following equations (Eqs. 11 and 12):

$$V^{*} left( S right),,, = ,,,frac{1}{18}left[ {n,,left( {n – 1} right),,left( {2n + 5} right)} right]frac{n}{ns*},,,,,,,,,,$$

(11)

$$frac{n}{ns*},, = ,,,1 + ,frac{2}{{n,,left( {n – 1} right),left( {n – 2} right)}},,,sumlimits_{i = 1}^{p} {left( {n – i} right),} ,left( {n – i – 1} right),,left( {n – i – 2} right),,,p_{s} ,left( i right),,$$

(12)

1998) modified MK trend test facility (at a significance level of 10%), to account for the autocorrelation effect.

Sen’s slope estimator computes the linear annual rate and direction of change (Sen 1968). It is a nonparametric approach for dealing with skewed datasets and outlier effects. The linear model f (t) is defined by the equations (Eqs. 13 and 14) (Sen 1968) as follows:

$$fleft( t right),, = ,,Qt, + ,beta ,,,,,,,,,,,$$

(13)

$$Q,,, = ,,Median,,,,,,{{left( {X_{i} , – ,X_{j} } right)} mathord{left/ {vphantom {{left( {X_{i} , – ,X_{j} } right)} {left( {i – j} right)}}} right. kern-nulldelimiterspace} {left( {i – j} right)}},,,,,,forall j,, < ,,i,,,,,$$

(14)

### Lake Hayk water level response to climate change/variability

Due to the endorheic (closed) nature of the Lake Hayk basin, the main underlying hydrological processes are surface runoff and evapotranspiration, with precipitation and temperature being the most prominent climatic factors. Under such conditions, water level is the primary response variable that serves as an indicator to better reflect the climate change/variability effects on lake storage. In addition, it can easily be measured at observation stations and the changes in lake water levels can be monitored easily, accurately and continuously (Tan et al. 2017). However, in situations where lake level data is patchy, as it is in Lake Hayk, remotely sensed water extents derived from Landsat images would allow us to bridge the data gap of the water level time series (McFeeters 1996; Xu 2006). This is achieved by developing spectral water indices to extract water bodies from remotely sensed Landsat images, which is typically accomplished by computing the normalized difference between two image bands and then applying an appropriate threshold to segment the results into two categories (water and nonwater features). The Modified Normalized Difference Water Index (MNDWI) can efficiently extract lake waters from Landsat images by easily suppressing signals from various environmental noises (such as vegetation, built-up areas and shadow noises) compared to its predecessor, the Normalized Difference Water Index (NDWI) using Shortwave Infrared (SWIR) rather than Near Infrared (NIR) used in the NDWI (Xu 2006). The formula used for the MNDWI calculation is:

$$MNDWI = ,,{{left( {Green, – ,SWIR} right)} mathord{left/ {vphantom {{left( {Green, – ,SWIR} right)} {left( {Green,, + ,SWIR} right)}}} right. kern-nulldelimiterspace} {left( {Green,, + ,SWIR} right)}}$$

(15)

For Landsat 5 Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper plus (ETM +), MNDWI becomes:

$$MNDWI = ,,{{left( {band2 – ,band5} right)} mathord{left/ {vphantom {{left( {band2 – ,band5} right)} {left( {band2, + ,band5} right)}}} right. kern-nulldelimiterspace} {left( {band2, + ,band5} right)}}$$

(16)

MNDWI values for water and nonwater features are determined by the reflectance ability of the features to band 2 and band 5 of Landsat 5 TM and Landsat 7 ETM + satellite images (Xu 2006). Water features have greater reflectance in band 2 than in band 5, resulting in positive MNDWI values, whereas nonwater features have negative MNDWI values (i.e., index values range from − 1 to + 1). A standard threshold value of zero is used to signify that a feature is water if MNDWI > 0 and nonwater if MNDWI ≤ 0 (Xu 2006). When the standard threshold value of zero is used, the MNDWI can accurately determine the spatial position of shorelines at the land–water boundary and successfully extract them from the multi-temporal Landsat TM and ETM + images. Thus, MNDWI values are still positive for the shallowest parts of water bodies and waterlogged areas (areas inundated with water). A change in MNDWI values at a specific time occurs when a sensor detects a spatio-temporal change in nonwater features (change in land use and land cover types within the lake basin) and/or change in depth and quality of water features.

Therefore, to bridge the time gap between 2005 and 2011, water level observations from 1999 to 2005 and 2011 to 2015 were fused with the MNDWI index extracted remotely sensed water areas from Landsat images in the ArcGIS 10.1 software environment. The analysis was conducted using cloud free (clouds cover ≤ 10%) Landsat 7 Enhanced Thematic Mapper plus (ETM +) satellite images of 2005 and 2008 and Landsat 5 Thematic Mapper (TM) satellite images from 2009 to 2011 obtained from the http://earthexplorer.usgs.gov/ portal. The obtained Landsat images (Level 1 Terrain Corrected (L1T) product) were pre-geo-referenced to the Universal Transverse Mercator (UTM) zone 37°N projection system using a World Geodetic System 84 (WGS84) datum. A single Landsat image is sufficient to encompass the entire Lake Hayk basin, which has an area of 8592.68 ha. Landsat 5 TM and Landsat 7 ETM + images specifications are shown in Table 1.

The MNDWI index was validated by correlating remotely sensed water areas to water level observations both obtained on the same date in each year by using the Pearson’s (parametric) and Kendall’s tau (nonparametric) correlations at a 0.01 significance level with the Statistical Package for the Social Sciences (SPSS) version 20 software. For n sample sizes of X and Y variables, the Pearson’s coefficient (r) can be computed as:

$$r,, = ,,{{sumlimits_{i = 1}^{n} {left[ {left( {X_{i} , – ,overline{X} } right),left( {Y_{i} , – ,overline{Y} } right)} right],} } mathord{left/ {vphantom {{sumlimits_{i = 1}^{n} {left[ {left( {X_{i} , – ,overline{X} } right),left( {Y_{i} , – ,overline{Y} } right)} right],} } {,left( {sqrt {sumlimits_{i = 1}^{n} {left( {X_{i} , – ,overline{X} } right)^{2} ,sumlimits_{i = 1}^{n} {left( {Y_{i} , – ,overline{Y} } right)^{2} } } } } right)}}} right. kern-nulldelimiterspace} {,left( {sqrt {sumlimits_{i = 1}^{n} {left( {X_{i} , – ,overline{X} } right)^{2} ,sumlimits_{i = 1}^{n} {left( {Y_{i} , – ,overline{Y} } right)^{2} } } } } right)}}$$

(17)

This can be confirmed by the nonparametric Kendall’s correlation. The Kendall’s correlation helps to minimize the effects of extreme values and/or the effects of violations of the normality and linearity assumptions (Kendall 1938). The Kendall’s tau (τ) is calculated based on signs as:

$$tau = ,,frac{2}{{n,left( {n – 1} right)}},sum {_{i,, < ,,j} ,sign,left[ {left( {x_{i} , – ,x_{j} } right),,left( {y_{i} , – ,y_{j} } right)} right]}$$

(18)

In both cases, the correlation coefficients are within − 1 and + 1. Correlation values close to ± 1 indicate strong relationships. For result interpretation, the hypothesis for a 2 tailed test of the correlation at a given significant level is defined as: H0: r, tau = 0 versus Ha: r, tau ≠ 0.