Combining multisource satellite data to estimate storage variation of a lake in the Rift Valley Basin, Ethiopia

Integration of remote sensing data sets from multiple satellites is tested to simulate water storage variation of Lake Ziway, Ethiopia for the period 2009-2018. Sixty Landsat ETM+/OLI images served to trace temporal variation of lake surface area using a water extraction index. Time series of lake levels were acquired from two altimetry databases that were validated by in-situ lake level measurements. Coinciding pairs of optical satellite based lake surface area and radar altimetry based lake levels were related through regression and served for simulating lake storage variation. Indices for extracting lake surface area from images showed 91–99 % overall accuracy. Lake water levels from the altimetry products well agreed to in-situ lake level measurements with R= 0.92 and root mean square error of 11.9 cm. Based on this study we conclude that integrating satellite imagery and radar altimetry is a viable approach for frequent and accurate monitoring of lake water volume variation and for long-term change detection. Findings indicate water level reduction (4 cm/annum), surface area shrinkage (0.08km/annum) and water storage loss (20.4Mm/annum) of Lake Ziway (2009–2018).


Introduction
Worldwide many lakes are dramatically dwindling or completely drying up as shown in Pekel et al. (2016); Wurtsbaugh et al. (2017) and Tan et al. (2018). For instance, Lake Urmia in Iran lost 86 % of its surface area and reached only 0.83 % of its original storage within 15 years (Jeihouni et al., 2017). In Tanzania, the surface area of Lake Manyara shrank by 94 % in a period covering only eleven years (Deus and Gloaguen, 2013). Lake Mead in USA lost nearly 40 % of its surface area within a decade (Forsythe et al., 2012). Water level of Lake Qinghai in China gradually declined over the past half of a century at a rate of 8 cm per annum (Chang et al., 2017).
In Ethiopia, Lake Haramaya and Lake Adele entirely vanished a decade ago (Alemayehu et al., 2007) and Lake Abiyata is currently under a massive threat (Seyoum et al., 2015). Similarly, Lake Ziway which is the fourth largest natural water body in Ethiopia has been severely affected due to increasing water abstraction competition by various sectors and individuals (Getnet et al., 2014;Shumet and Mengistu, 2016;Desta and Lemma, 2017). There is urgent need to evaluate the change of lake water storage and to provide evidence for better management of the lake.
As shown in Bing et al. (2011); Chang et al. (2015) and Treichler et al. (2019), monitoring variation of lake water storage is an indispensable indicator to identify the long-term impact of climate and land use change and human activities on the inland water resources. Traditionally, water storage of lakes has been monitored using groundbased bathymetric survey data and lake level recordings (Rientjes et al., 2011;Rakhmatullaev et al., 2011;Han et al., 2016). In many developing countries, reliable data sets however are often not recorded or scarcely available (Coe and Birkett, 2004;Alsdorf et al., 2007;Munyaneza et al., 2009). The time delay to disseminate collected data is another major constraint to conduct scientific studies (Zhang et al., 2006). Moreover, data sets may be incomplete with missing observation records or of poor quality for long term change assessments (Biancamaria et al., 2010). Therefore, development of novel lake monitoring approaches with alternative data sources to overcome the limitations of in-situ data is essential and motivated this study.
Currently several types of remotely sensed data sets are available for studying lakes from different perspectives (Politi et al., 2016;Dörnhöfer and Oppelt, 2016). Since the launch of the first mission in 1970s, Landsat series are the most widely used optical remote sensing data for time series mapping of lakes surface area (Work and Gilmer, 1976;Sharma et al., 1989;Cui et al., 2017). Various water indices are developed by algebraically combining and taking a ratio of two or more bands of the imageries (e.g. McFeeters, 1996;Xu, 2006;Feyisa et al., 2014). These indices are widely applied to distinctively identify, extract and monitor the surface area of water bodies (Sun et al., 2017;Ogilvie et al., 2018;Schwatke et al., 2019). These indices mainly combine the green, near infrared and shortwave infrared bands to compose a single band gray-scale image. The bands are selected because water more strongly absorbs radiation within these bands (Li et al., 2013;Xie et al., 2016) than land features.
Selection of a specific water index (WI) depends upon the accuracy of the extracted lake surface area when compared to reference data sets. The accuracy of WI can be affected by location and time of image acquisition (Jiang et al., 2014), uncertainties from applied atmospheric correction method (Ayana et al., 2015), turbidity of water and class of the background land cover (Sharma et al., 2015) and type of sensor (Zhou et al., 2017). Due to these factors, there is consensus that there is no single one specific WI that is robust, stable and efficiently applicable for all water bodies under all circumstances. As a result, indices should be compared to select the one which most accurately extracts the surface area of any specific lake (Rokni et al., 2014;Taravat et al., 2016;Moknatian et al., 2017;Buma et al., 2018).
Radar altimetry is a remote sensing data source that is used for more than quarter-century for monitoring water level of lakes (Alsdorf et al., 2001;Ričko et al., 2011;Hwang et al., 2016;Wang et al., 2019). The data is collected from satellites that are equipped with sensors that continuously emit microwave signal towards the surface of the earth. The sensors register two-way traveling time of the emitted signal and of the reflected echo. The registered data helps to obtain satellite-to-water surface distance and later to estimate the lake water level. However, the data always requires specific correction for delay of the signal at different layers of the atmosphere and for crustal motions (Fernandes et al., 2014;Boergens et al., 2016;Quartly et al., 2019).
Currently, radar altimetry satellites are capable of effectively monitoring water levels of lakes larger than 50 km 2 in size (Crétaux et al., 2016). The reported measurement error of lakes water level from radar altimetry is highly variable (e.g. Jarihani et al., 2013;Okeowo et al., 2017). The error commonly ranges from less than 10 cm (e.g. Villadsen et al., 2016;Crétaux et al., 2018) to tens of centimeters (e.g. Kleinherenbrink et al., 2015;Nielsen et al., 2017) and beyond in rare cases (e.g. Liu et al., 2016). The major factors governing the accuracy of altimetry derived lake level data are (i) length of footprint over the lake that depends on shape and size of the lake and location the footprint passes over, (ii) complexity of terrain surrounding the lake that alters the magnitude of the reflected echo, (iii) height of wave on the surface of the lake that depends on wind speed over the area and (iv) method of processing the raw altimetry data called waveforms (Mercier et al., 2002;Birkett and Beckley, 2010;Ričko et al., 2012). According to Fekete and Vörösmarty (2007), gradually changing land cover, smooth transition of topography and wide water surface area are ideal conditions to obtain highly accurate altimetry measurements.
Satellite imagery provides information about estimation of water surface and lake shorelines whereas radar altimetry provides measurements of water levels. Information from each satellite data source is insufficient to comprehensively describe storage variations of lakes. As a result, the combined use of both satellite imagery and radar altimetry to study lake storage variation is advocated as shown in Singh et al. (2015) and Busker et al. (2019). However, the available studies mostly targeted large lakes (with surface area larger than 1000 km 2 ) that geographically are concentrated in the Middle East (e.g. Sima and Tajrishy, 2013;Muala et al., 2014;Nguy-Robertson et al., 2018) and in the Tibetan Plateau (e.g. Kropáček et al., 2012;Song et al., 2013;Li et al., 2019). In most of similar studies (e.g. Arsen et al., 2013;Avisse et al., 2017;Fang et al., 2019) in-situ based lake data is either not available or used to limited extent to reveal accuracy of the remotely sensed data sets for estimating respective lake variables and thus to track changes in lake water storage over time.
In Africa, studies mostly relied on altimetry data sets only and assessed water storage change of the African Great Lakes (e.g. Swenson and Wahr, 2009;Hassan and Jin, 2014;Moore and Williams, 2014). These lakes are Victoria, Tanganyika and Malawi that have a surface area greater than 150,000 km 2 . Similarly, other large lakes of the continent such as Lake Turkana (Velpuri et al., 2012) and Lake Chad (Buma et al., 2016) were also studied by applying satellite radar altimetry only. Duan and Bastiaanssen (2013); Tong et al. (2016) and Zhu et al. (2017) are among the very few studies that reported storage variation of largest lakes of Africa by integrating satellite imagery and altimetry. As such, there is a scope to evaluate integration of freely available satellite imageries and synergizing them with radar altimetry data to monitor storage variations of medium and smaller lakes of the continent. Lakes which are situated in the Eastern Africa Rift (EAR) system are not well monitored for lake water storage changes. Among the lakes in EAR basin, Lake Ziway is selected for this study by availability of in-situ based lake level observations time series. Temporal change of this medium sized lake is not fully reported due to lack of lake surface area data.
This study uses in-situ data as a reference for evaluation of accuracy of combined optical and radar altimetry data sets and for examining water storage variation of Lake Ziway. First, we systematically compared lake surface area extracted from Landsat image using various water indices to select the most appropriate index. Next, the accuracy of water level time series data which is obtained from different altimetry databases was compared against in-situ measured water levels. Through integrating both remotely sensed data sets, volumetric variation of the lake is estimated at time instants at which coinciding data pairs are available and validated using observed data collected from bathymetric survey. A temporal variation of the lake was detected from the chronologically assembled surface area, water level and storage variation time series data. The materials and approaches which were applied in this study can be adopted for lakes of similar size that suffer from data scarcity.

Geographic settings and location
The endorheic Central Rift Valley (CRV) basin of Ethiopia is part of the Great East African Rift Valley system. The basin encompasses four medium sized lakes named Ziway, Langano, Shala and Abyata (Fig. 1). Unlike other CRV lakes, Lake Ziway is the largest in size and the only freshwater lake in the CRV basin. The mean and maximum depths of the lake are about 2.5 m and 9 m respectively (Belete et al., 2016).
The lake is geographically the most upstream in the CRV basin at about 08̊00′ N and 39̊50′ E. This water body is located at 160 km south of the capital, Addis Ababa. It is situated on the left side of the major road heading to Nairobi and in vicinity to Ziway Town. Five islands, which are known as Tulugudo, Gelila, Funduro, Tedecha and Debre Sina are present on the lake and mostly serve as inhabitation for the Zay people.

Physical and hydro-climatic characteristics
Lake Ziway is fed by two perennial rivers named Katar and Meki with annual average inflow of about 464 and 276 million cubic meters, respectively (Gadissa et al., 2018). These two rivers originate from Arsi and Gurghe highlands, in the western and eastern escarpments of the lake catchment respectively. Annually around 170 million cubic meters of water is discharged from the lake (Getnet et al., 2014;Goshime et al., 2020) into Lake Abiyata via Bulbula River.
According to the Köppen-Geiger climate classification (Peel et al., 2007), the lake catchment area is categorized under dry-temperate climate. The lake is receiving mean annual rainfall of 750 mm and has average daily temperature of 21C. Radiation and wind blown over the surface of the lake results in 1660 mm of evaporation annually (Melesse et al., 2009).

Optical satellite imagery
In this work, multi-temporal Level 1 T (terrain corrected) Landsat ETM+/OLI images of Lake Ziway were obtained from USGS-GloVis (http://glovis.usgs.gov) in GeoTIFF file format. Selected images had minimal cloud cover to capture the surface area of the lake. The area of the lake is covered by adjoining areas of two Landsat scenes: Path168/ Row054 and Path168/Row055 (Fig. 1). To reduce radiometric error during mosaicking, images that were acquired at the same date were downloaded and projected to UTM Zone 37 N using WGS84 datum.
The selected images were mainly captured in January, March, June and October, which represent the major inflectional points on the normal annual water level pattern of the lake. Images were identified by tolerating only minor difference (less than 20 days) between the image acquisition dates for different years. However, images were missing occasionally due to either unavailability or were excluded because of poor visual quality. Cloud cover was the other considered factor while identifying candidate images for analysis. Due to low reflectance, clouds are commonly miss-classified as a water feature (Ji et al., 2015) and subsequently deteriorate the accuracy of lake surface extraction from satellite imageries. This selection procedure resulted in, a time series of sixty cloud-free or nearly cloud free scenes that were further processed for analysis in this study.

Radar altimetry
Water level time series data of Lake Ziway was obtained from two satellite radar altimetry data sources and from a lake level gauging station. The altimetry sources were Global Reservoir and Lake Monitor (GRLM) (Birkett et al., 2011, https Both organizations avail the products after processing radar signals which are collected from different altimetry missions. In case of Lake Ziway, the waterbody is overpassed by Topex/Poseidon (TP) and Jason family radar altimetry missions that have 10 days repeat cycle and 315 km ground track spacing. Ground footprint with pass number 094 of these missions crosses a stretch of 11 km over the surface of Lake Ziway from north to northeastern direction ( Fig. 1). Daily recorded insitu lake level (2009)(2010)(2011)(2012)(2013)(2014)(2015) and bathymetric survey data (2006) for validation purpose were retrieved from the Ministry of Water, Irrigation and Electricity (MoWIE), Ethiopia. The in-situ water level data were recorded at Ziway Station which is located at the southern tip of the lake (Fig. 1).

Reference data sets
In this study, the surface area of Lake Ziway is extracted from Landsat imagery using three frequently applied water indices. In our approach, the most accurate index was identified to delineate the lake extent for lake surface estimation. The reference data set for evaluation comprised five hundred ground control points (GCP) that were collected from the southwestern portion of the lake (Fig. 1). The GCPs were collected during a field campaign using GPS device and digitized from high resolution image (HRI).
Four hundred of the GCPs were gathered from both water body and land features surrounding shore of the lake. We covered representative features including water surface from middle lake area and mixed water surface around edges of the lake, built-up areas, bare land and forest. The HRI was used to collect the remaining hundred GCPs from unreachable features in the lake surrounding because of access restrictions. These features include wetlands, floriculture shade nets and mechanized farms. Therefore, GeoEye-1 image with resolution of about 2 m 2 was freely obtained from GoogleEarth™ platform. The image was acquired and made available through ©DigitalGlobe.

Image pre-processing
Satellite images are always exposed to different type of distortions (Young et al., 2017). These distortions can be caused by malfunctioning of sensors, solar and atmospheric interferences, topographic errors and geometric inaccuracies. Therefore, preprocessing was done prior to application to correct and improve representation of the targeted feature.
First, we applied radiometric calibration as the initial stage in the series of image preprocessing stages. This calibration was done to extract more meaningful information regarding the magnitude of object reflectance arrived at the satellite sensor. Based on the input requirement of subsequent processes, it involves a conversion of digital number (DN) into one of either output: at the sensor radiance or top of atmosphere reflectance. In this study, the DNs were converted into radiance for feeding into the succeeding image correction stage. The required input variables for the conversion were obtained from the metadata (MTL.txt) file which is made available within the decompressed Landsat imagery folder.
Next, Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes, FLAASH (Atmospheric Correction Module, 2009) was applied on the resulting image to obtain surface reflectance. This correction helps to reduce the effect of scattered and absorbed light due to airborne matters in the atmosphere. We selected FLAASH because of its wide application (e.g. Feyisa et al., 2014;Yang et al., 2015;Xiong et al., 2018) and its user defined input inquiry that considers multiple parameters and a user defined range of wavelengths. The tool offers the ability to correct from visible (0.4 μm) up to shortwave infrared (3 μm) range of wavelength which is the common range for water body extraction.
Wedge shaped data gaps occur on Landsat 7 ETM + images since May 2003 due to failure of Scan Line Corrector (SLC). As a result, null data strips covering up to 22 % of the entire scene are formed on the images of the product. We filled these missing pixels using single file gap fill extension toolbox in ENVI which applies a triangulation interpolation method. Mosaicking of scenes was also done using ENVI Classic mosaicking tool which is based on geo-referenced images to compose an image that covers the lake area.
We subsequently resized all images to our area of interest that is Lake Ziway. A base image was selected and geometrically corrected using twenty spatially distributed ground control points (GCPs) to coregister multiple time series images into nearly exact geographic location. These GCPs were collected at distinguishable features including island, bridges, main road junctions, intersection and curves using hand-held GPS. The accuracy of the correction was assessed using RMSE by keeping the margin of an error below half a pixel size (15 m) and by removing four outlier points. Lastly, the entire images time series were co-registered using the geometrically corrected base image with an error of less than one-fourth of a pixel.

Water extraction indices and thresholding
We applied three water extraction indices on the pre-processed Landsat images to the extract surface area of Lake Ziway (Table 1). Water extraction indices combine selective bands of the satellite image (Table 2) which are available within the visible to infrared wavelength  range. The bands were used because of a more pronounced absorption and less reflectance property of water in the mentioned wavelength range. This resulted in a gray scale image having land and water feature representing bimodal type of pixel value distribution histogram. Here, an optimal threshold value demarcating the periphery of water surface from the adjacent land features is crucial to extract the lake surface area. Conventionally, a null value is taken as a default threshold for quick and ease separation of the two classes. However, the applicability of such thresholding approach could not be always effective and consistent. The location of the study area and time of image acquisition are perpetually causing variation of the threshold value and hence adjustment is required for each local condition.
Currently, various thresholding methods are available to convert the gray scale images into a binary (black and white) image that represents the two classes distinctively. In this study, one of the most commonly used methods called Otsu's Thresholding Method (Otsu, 1979) was applied for Normalized Difference Water Index (NDWI) and Modified NDWI (MNDWI) using MATLAB R2014b (MathWorks, 2014). The method was selected due to its time effective, automated and nonsubjective nature. In Otsu's method, the value in the histogram that reduces variability within a class and maximizes the variance between classes is selected as an optimal threshold. However, the suitability and effective applicability of Otsu's method is restricted to the case of seeking a single optimal threshold in the valley of a bimodal histogram having equivalent peaks (Fan and Lei, 2012) which is not a case for Automated Water Extraction Index (AWEI) histograms.
As a result, studies applied AWEI used the zero threshold value (e.g. Zhai et al., 2015;Sunder et al., 2017) because of its stability at different locations (see Feyisa et al., 2014). In contrary, other studies (Satgé et al., 2017;Acharya et al., 2018) revealed the deviation of AWEI threshold values from the recommended. Therefore, the Modified Histogram Bimodal Method (MHBM) was applied in this study for thresholding AWEI gray scale image following the procedure mentioned in Zhang et al. (2018). The method defined an optimal threshold as a minimum median value in the valley of a bimodal histogram formed from an area one and half times greater than the water body it encompasses.

Lake surface area accuracy evaluation
In this paper, accuracy of the extracted lake surface area from the three indices was evaluated using GCPs collected from the field and from a high-resolution image. The evaluation process included the following tasks: (i) selection of reference image, (ii) sample point generation from the reference data sets and (iii) preparation of a confusion matrix.
The selection of reference image was conducted by considering a minimal gap of acquisition date of two days with time reference to the Landsat image. This helps to ensure temporal consistency between the different data sources and to reduce the influence of time-born effects. The selected reference image encompasses an area of 120 km 2 and covered the southwestern portion of Lake Ziway with diversified features including water body, built-up areas, irrigated area and wetland (Fig. 1).
Subsequently, five hundred sample points were generated from both reference data sets. The samples were collected equally from both water and land cover features. Stratified random sampling technique was used to collect the samples from the diversified land covers based on their extent weightage on the reference image. The type of class for each sample point was assigned through on-screen digitization process based on expert knowledge by a site visit. Next, a comparison was conducted between the sample points and their respective point on the binary images of the three indices. After comparison, each point was categorized in one of the four bins based on the two classes (water and land) in the classified and reference data sets. These bins are Water classified as Water (WW), Water classified as Land (WL), Land classified as Water (LW) and Land classified as Land (LL).
Lastly, a confusion matrix, which has elements representing the number of counts that fall under each of the four bins, was developed for accuracy assessment. As shown in Table 3, we applied the overall accuracy (OA), producer accuracy (PA) and user's accuracy (UA) to quantitatively evaluate the accuracy of the classified images (Story and Congalton, 1986). The non-dimensional kappa coefficient (κ) was also computed (Eq. 1) to measure the extent of agreement between the classified image and reference data set (Congalton et al., 1983). TW TW TL TL  T  TW TW TL  (1) where: Ts stands for total number of samples; Tc and Tr represent total number of classified and reference data sets respectively; and W and L stand for water and land pixels respectively.
Among the indices, the best performed index having the least error was selected and adopted to extract multi-temporal surface area of Lake Ziway. These surface areas were later chronologically arranged to quantify the rate of variation from the general slope of the trend. Vector files of the surface areas were superimposed to visually detect changes on water shorelines of the lake.

Validation of radar altimetry
The lake level measurements which were obtained from in-situ and radar altimetry (GRLM and DAHITI) databases were measured with respect to different references. The in-situ and GRLM lake water level data sets were available in height variation from a fixed reference whereas DAHITI database provided an ellipsoidal height that is measured with respect to the WGS84 datum. The height variation obtained from GRLM database was measured with respect to Jason-2 mission reference pass level.
In this study, we selected WGS84 datum as a transformation frame to conduct a consistent and uniform comparison between the water level measurements from the three data sources. The datum was selected because of its wide use. Therefore, vertical shifting constants were added to in-situ and GRLM measurements that resulted in lower difference with DAHITI time series.
The accuracy of each altimetry water level product was evaluated using the in-situ lake level measurements as reference data. Mean Absolute Error (MAE), Root Mean Square Error (RMSE) and Coefficient of Determination (R 2 ) were used here to quantify the accuracy of these remotely sensed water levels (Willmott, 1982;Ji and Gallo, 2006;Chai and Draxler, 2014). The equations read as follows: where WLalt and WLgrd are water level measured by satellite altimetry and ground gauging station and i is time index for water level observation and n is total number of observations. Over-bars indicate the mean of each water level data sets. Lastly, the rate of lake level variation was quantified from the best fitting linear line slope of the time series measurements.

Estimation and validation of water volume variation
Water storage in the lake can be regarded by a reference water level that separates static water storage (i.e., lower storage) from dynamic or temporally varying water storage (i.e., upper storage). The identified reference is defined as the lowest lake level (LLL) that was recorded during the study period. Thus, the lower and the upper storage correspond to the volume of water that is stored below and above the LLL, respectively.
Since the data from satellites that are used in this study are restricted to the water surface, the applicability of the data is limited to estimating the upper storage of the lake only. Estimation of the lower storage is challenging because of (i) uneven nature and rapid changes of the lakes bottom shape in response to sedimentation and (ii) incapability of optical and radar altimetry satellite sensors to capture the lower storage. In this study, we estimated the upper storage variation by developing a second order polynomial equation (Eq. 5). The equation correlates the multi-temporal lake surface area (LSA) of the lake which is extracted from Landsat with the corresponding date radar altimetry relative lake level (RLL) that is obtained after subtracting LLL from the lake water level (LWL) time series.
Here, exact coinciding dates are not always attained for the two satellite data sets because of differences in revisit time interval between Landsat (16 days) and radar altimetry missions (10 days). Therefore, studies commonly allow 5-10 days difference to develop a remote sensing based water level-surface area relationship (Baup et al., 2014;Zhu et al., 2014;Tong et al., 2016). In this study, we selected paired images and altimetry measurements having less than eight days of acquisition differences. We reached at this threshold after setting three criteria and analyzing the optimal values that satisfied the considered factors. The considered criteria are (i) a minimal date difference, (ii) highest number of paired datasets and (iii) maximum possible coefficient of determination (Fig. 2). The polynomial equation that fitted the paired data points was used to estimate the lake surface area (LSA) and subsequently, integrated to estimate lake volume variation (LVV) (Duan and Bastiaanssen, 2013). The equations read as follows: Acquisition date difference (DD) between paired lake surface area and water level measurements with respective coefficient of determination which resulting from regression analysis for developing equations relating water level of the lake to surface area and to water volume variation.
W. Asfaw, et al. Int J Appl Earth Obs Geoinformation 89 (2020) 102095 where: a, b and c are regression constants of the second order polynomial that correlates relative lake level (RLL in m) that is measured above the lowest lake level (LLL in m) and lake surface area (LSA in km 2 ) of the same date to obtain lake volume variation (LVV in Mm3). Note: the constant is specified as zero with the condition of LVV = 0 at RLL = 0, but LSA ≠ 0. In this study, bathymetric survey data of the lake by MoWIE in 2006 is used as a reference to validate accuracy of the estimated LVV. We quantified accuracy of our estimation with three statistical indicators (RMSE, MAE and R 2 ) as previously described in altimetry data validation section. Subsequently, we generated the volumetric variations corresponding to different periods using the integrated equation and chronologically assembled to represent temporal storage variation of the lake.

Accuracy of water extraction indices
The accuracy of the three water indices (NDWI, MNDWI and AWEI) was evaluated by use of five hundred ground control points and by inter-comparison of the lake extent extracted from the indices (Fig. 3). The reference dataset was obtained through a field survey and from high resolution image (HRI). Visual interpretation of lake surface area using the HRI indicated that MNDWI and AWEI best distinguished the water surface from the adjacent land features. In contrast, NDWI showed limitation in uniquely detecting the water body from some of the land features, since some urban area with higher reflectance (e.g. concrete surface on top of buildings and new metal roofing materials) were misclassified as water body.
The classified images of MNDWI and AWEI showed high coincidence and similarity. Hence, visual comparison is not adequate to compare the performance of the two indices. Therefore, confusion matrix based accuracy indicators were applied to compare results of the two indices (Table 4). The comparison shows that both indices performed well with higher than 98 % for all accuracy indicators and 0.98 for kappa coefficient. However, AWEI performed slightly better (by +0.2 % for overall accuracy and +0.4 % for κ) than MNDWI in distinguishing water feature. As a result, AWEI was preferred for this study to extract multi-temporal surface area of Lake Ziway from satellite images. We note that, MNDWI showed adequate accuracy and thus we consider it appropriate for Lake Ziway. NDWI showed to have lowest value that we considered too low to accept related outcomes. Fig. 3. The reference and classified images covering the southwestern part of Lake Ziway and nearby land feature.

Validation of radar altimetry data
In this study, two radar altimetry products obtained from GRLM and DAHITI databases were validated using in-situ lake level measurements. For comparison purposes, vertical shifting constants of 1645.03 m and 1646.15 m were added to in-situ and GRLM water levels, respectively, for measuring the lake height from WGS84 datum as DAHITI. These constants are obtained by defining an added value that results in zero difference when mean value of DAHITI measurements is subtracted from mean values of in-situ and GRLM datasets, respectively.
The result showed that DAHITI and GRLM altimetry products well agreed to the measured water levels of Lake Ziway. In the validation process, these products strongly agreed with in-situ measurements and resulted in similar accuracy but with slight superiority of GRLM. The combined mean of both data sets yielded improved results as compared to individual data sets (Fig. 4). The use of the combined mean value resulted in an improved coefficient of determination (by +8.2 % for R 2 DAHITI and +7 % for R 2 GRLM ) and reduced errors (by -29.2 % for RMSE DAHITI , -24.5 % for RMSE GRLM ; -30.3 % for MAE DAHITI and -25.2 % for MAE GRLM ). Villadsen et al. (2016) describe that small errors can result from three main sources of error that are (i) the geophysical corrections applied for radar altimetry, (ii) the difference in time of water level recording and (iii) distance between location of in-situ gauging station and ground track of radar altimetry missions. For this study, effects of the third cause can be considered as minor due to relatively flat water surface of Lake Ziway. In general, lake level measurements of radar altimetry databases well matched to counterparts by in-situ measurements. As a result, the water level of Lake Ziway can be monitored at high accuracy by combining the radar altimetry satellite products.

Estimation and validation of lake volume variation
In this study, regression constants which are required for estimation of the time dependent LVV, were obtained from twenty five paired data points. As shown in Fig. 2, these pairs of surface areas (LSA) from Landsat ETM+/OLI and relative lake levels (RLL) from radar altimetry have a maximum of eight days difference in acquisition dates. Coefficient of determination for a second order polynomial equation for remotely sensed LSA and RLL data sets resulted in a value of 0.91 (Fig. 5).
This suggests strong relation and hence, the equation is used to establish lake surface area data from available water level data. The equation can also be used in case the lake surface area extracted from satellite images is known but related lake water level is unknown. This allows to reconstructing water storage changes for past dates when lake water surface area is extracted by available Landsat ETM+/OLI images. Finally, we develop the Eq. 10 to estimate LVV of Lake Ziway at different RLL measurements above the lowest lake level, LLL (i.e. 1645.43 Validation of estimated LVV was conducted by use of bathymetric survey data developed by MoWIE in 2006. Validation resulted in very high coefficient of determination (R 2 ) of 0.9956. In terms of error indicators, MAE of 47.2 Mm 3 and RMSE of 56.8 Mm 3 account only for 2.9 % and 3.5 % of the lake total storage (as reported on WWDSE and CECE, 2008), respectively. The error potentially is caused by (i) errors that must be associated with the applied satellite data sets and regression technique (Van Den Hoek et al., 2019), (ii) inaccuracy of the bathymetric data (Ayana et al., 2015), and (ii) the historic time period of the bathymetry survey (i.e. 2006).

Lake storage variation
In this study, we established relations between water level (LWL), surface area (LSA) and volume variation (LVV) to unveil decadal variation of Lake Ziway. Findings show that, the lake water level is fluctuating annually in response to the seasonal rainfall pattern. Based on the long-term data , oscillation of LWL showed highest correlation with rainfall for a lag time of three months (Fig. 6).
For the period 2009-2018, fluctuation of the lake can be characterized by three major phases (Fig. 7). The first phase was the period of highest lake levels in October 2010. In the second phase, from mid of 2011 up to mid of 2015, the lake has consistent water level fluctuations with gradual changes that reflect on ordinary lake inflow and outflow conditions. However, in the third time span that extends from June 2015 to the end of 2018, the lake did not attain highest lake water storage at the end of the wet season. In this period, LWL, LSA and LVV W. Asfaw, et al. Int J Appl Earth Obs Geoinformation 89 (2020) 102095 consistently remained at a relatively lower level than the previous years.
In terms of anomalies, the lake attains its highest and lowest level in October 2010 and February 2016 respectively. Highest lake level was 1647.43 m with a surface area and volume variation of 417.74 km 2 and 791.97 Mm 3 respectively. The high annual rainfall was a major factor that caused the high lake level. In February 2016, Lake Ziway indicated the LLL of 1645.33 m, the smallest surface area of 414.35 km 2 and LVV reached the null level. This situation potentially results from the 2015−2016 El Niño Southern Oscillation (ENSO) event which affected water resources and caused devastating drought in East African region (Siderius et al., 2018;Kolusu et al., 2019) and in Ethiopia particularly.
From the analysis for the period of study (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018), an average annual lake level decline of 0.04 m and surface area decrease of 0.08 km 2 was shown. This converts to an average annual lake water storage reduction of 20.4 Mm 3 . These findings indicate that over a period of 10 years, Lake Ziway lake level has dropped 0.4 m and the lake area reduced by 0.8 km 2 to cause a water storage loss of 204 Mm 3 . Based on the reported 1.6 Bm 3 water storage of the lake in 2008 (WWDSE and CECE, 2008); the lake lost some 12.75 % of its storage volume over the past decade.
As shown in Fig. 8, the northern shore of the lake has retreated from north to the south direction. However, the magnitude of change in terms of surface area was not significant. Therefore, the volume loss can be attributed to lowered lake water levels. This shows that, the water volume variation more relates to lake level change rather than the changes in lake surface area. Similar finding is also shown in Gao et al. (2012); Ye et al. (2017) and Keys and Scott (2018).

Impacts of anthropogenic activities and climatic variables
Some water storage loss is noticed in Lake Ziway over the past 10 years. The loss is potentially driven by increased lake water abstractions, by climatic changes and the El Niño event in 2015-2016 which caused substantial water stress and drought in the Eastern Africa region (Anyamba et al., 2018;Qu et al., 2019). Rapid expansion of irrigation based local and commercial farming is observed on the northern and western adjacent areas of the Lake Ziway. Water for these command areas is abstracted from the lake and contributed to reduction of the lake water storage; but, quantitatively not further assessed in this study. Effects of climate change and the El Niño event on lake water storage also are not further evaluated in this study. However, to better understand long term lake water balance implications we recommend assessment on rainfall and temperature following (Haile and Rientjes, 2015), and to perform hydrological impact assessments on runoff following (Haile et al., 2017).  W. Asfaw, et al. Int J Appl Earth Obs Geoinformation 89 (2020) 102095 In addition, sedimentation has a substantial impact on water storage capacity of the lake. In case of Lake Ziway, Meki and Katar rivers are main sources of sedimentation that contribute 57.2 % and 42.8 % of sediment inflow, respectively (Gadissa et al., 2018). Annually, 2.04 Mton of sediment having a depth of 4 mm is deposited in the lake resulting in annual storage loss of 0.11 % (Aga et al., 2019). However, unlike the aforementioned factors water loss due to sediment accumulation is not easy to identify. Its effect is significant on the lower storage of the lake rather than on time dependent storage variation of the lake. Based on sedimentation accumulation rate mentioned in Aga et al. (2019), water storage reduction of Lake Ziway is higher than shown by findings in this study.
Reduction of lake water storage may directly affect two salient components of the ecosystem. First, Lake Ziway drains to Lake Abiyata through Bulbula River that hydraulically connects both lakes. As a consequence, lower storage of Lake Ziway affects inflow of Lake Abiyata that is the CRV terminal lake. Secondly, the lowering of Lake Ziway water level may reduce fish diversity and bird population of the  W. Asfaw, et al. Int J Appl Earth Obs Geoinformation 89 (2020) 102095 lake because of changes in physical properties and chemical characteristics of the water. Abera et al. (2018) showed that fish production reduced from to 3180 to 1157 tons in a period between 1997 and 2014 although effects of reduced water storage remain uncertain. In general, reduction of the lake storage impacts the environment and socioeconomic activity of various stakeholder groups. We would like to note that the storage loss would have been more than reported in this study if it was not for the new barrage which is constructed at the outlet of Lake Ziway in 2016.

Summary and conclusions
This study evaluated variation of the lake storage using data sets from multiple satellites and ground based data. Multi-temporal Landsat ETM+/OLI optical satellite imageries were processed to extract lake surface areas using three indices. Accuracy assessments indicated that the performance of AWEI and MNDWI indices was high whereas that of NDWI was relatively low. Lake level data from two satellite altimetry databases was compared to ground measured water level data. The obtained result of R 2 = 0.92 and root mean square error of 11.9 cm was within the globally reported range of high accuracy (Sulistioadi et al., 2015;Schwatke et al., 2015).
Paired radar altimetry water level and Landsat ETM+/OLI based lake surface area images were related using regression and resulted in an equation that can be used to estimate water volume variations. Availability of a bathymetric survey allowed validation of the estimated lake storage variation and the result portrays a higher accuracy. Approaches used in this study can be applied to other similar sized lakes which have not been monitored due to scarcity of in-situ lake water level and bathymetric survey data.
Water levels, surface areas and volume variations were chronologically assembled to reproduce time series to indicate water storage variation of the lake. From the time series data, we showed that the lake reaches its highest and lowest lake water storage values by wet and dry climatic anomalies which occurred in 2010 and 2015 respectively. This study concludes that lake levels, lake surface area and lake volume gradually decrease at an annual rate of 0.04 m, 0.08 km 2 and 20.4 Mm 3 respectively. As a result, since 2009 the lake has lost 12.75 % of its water storage volume.
The decrease of lake water storage can be caused by increased lake water abstractions, climatic change and effects of the El Niño event in 2015-2016 which caused a substantial water stress and drought in the Eastern Africa region. This study did not investigate impacts by each of these causes; but findings indicate that further loss of storage or reduction in outflow by the operation of the newly constructed barrage can be anticipated. In consequence, socioeconomic activity in surrounding areas will be interrupted and aquatic life and the downstream ecosystem will be endangered. Findings of the study serve decision makers to establish water resource management intervention strategy to assure sustainability of the lake.