Remote Sensing of Storage Fluctuations of Poorly Gauged Reservoirs and State Space Model (SSM)-Based Estimation

Abstract

To reduce hydrological uncertainties in the regular monitoring of poorly gauged lakes and reservoirs, multi-dimensional remote sensing data have emerged as an excellent alternative. In this paper, we propose three methods to delineate the volume of such equipotential water bodies through a combination of altimetry (1D), Landsat (2D) and bathymetry (2D) data, namely an altimetry-bathymetry-volume method (ABV), a Landsat-bathymetry-volume method (LBV) and an altimetry-Landsat-volume-variation method (ALVV). The first two data products are further merged by a Kalman-filter-based state space model (SSM) to obtain a combined estimate (CSSME) time series and near future prediction. To validate our methods, we tested them on the well-measured Lake Mead and further applied them on the poorly gauged Aral Sea, which has inaccurate bathymetry and very limited ground observation data. We updated the lake bathymetry of the Aral Sea, which was more than half a century old. The resultant remote sensing products have a very good long-term agreement among each other. The Lake Mead volume estimations are very highly coherent with the ground observations for all cases (R² > 0.96 and NRMSE < 2.1%), except for the forecast (R² = 0.75 and NRMSE = 3.7%). Due to lack of in situ data for the Aral Sea, the estimated volumes are compared, and the entire Aral Sea LBV and ABV have R² = 0.91 and NRMSE = 5.5%, and the forecast compared to CSSME has R² = 0.60 and NRMSE = 2.4%.

1. Introduction

Lakes and reservoirs are vital not only because they are a major source of water for domestic and industrial usage for human beings, but also because of their riparian zones, where some of the most bio-diverse ecosystems exist. Eventually the rapidly changing water volume of lakes and reservoirs disturbs not only human settlements that are dependent on them, but also whole ecosystems. Therefore regular monitoring and dynamic modeling of water volume variations are crucial for the proper management of this limited but highly essential resource. On a global scale, relatively few water bodies are regularly monitored by traditional in situ measurements. In the last few decades, satellite remote sensing has evolved as a promising alternative for regular global monitoring of water resources [1,2,3]. Satellite altimetry is now a well-established tool for inland water level estimation [4,5,6,7] and Landsat, with its long archive, free availability and relatively high-resolution database, delivers one of the most frequently used remote sensing data sets [8,9,10,11].

In this paper, we explored the potential of multi-dimensional satellite remote-sensing-based reservoir water volume estimation methods. We evaluated the monthly volume dynamics of the Aral Sea and Lake Mead using three independent methods. First, we used the altimetry-bathymetry-volume (ABV) method, which intersects the digital elevation model (DEM) of the reservoir bathymetry with the water level time series generated from satellite altimetry observations [12,13]. Secondly, we used the Landsat-bathymetry-volume (LBV) method, which is relatively less explored. In this method, the water height time series is generated by the intersection of the bathymetry of the reservoir with the temporal evolution of the shorelines generated from the Landsat images. Andreoli et al. [14] integrated a 90 m SRTM (Shuttle Radar Topography Mission) DEM with 75 m–1 km ENVISAT images and averaged the water levels at the boundary buffer to derive water height. In other previous studies, in situ water levels were combined with satellite-derived areas to obtain the volume [13,15,16]. Thirdly, we used the altimetry Landsat volume variation (ALVV) method by combining the water level from satellite altimetry with the respective surface area of the reservoir to obtain the volumetric variations [16,17,18,19]. Other groups [16,20] used a similar approach of estimating volumes by three methods; however, they either combined in situ observations or used existing elevation-volume relationships or storage capacities of the reservoirs from different databases. In this paper, we present a new LBV method for reservoir volume estimation and compare it to the results obtained from the ABV and ALVV methods. Furthermore, we design a Kalman-filter-based state space model to combine the estimated absolute volume from LBV and ABV. In our study, we did not use any ground observation information other than the bathymetry of the water body. If bathymetry is not available then one can extrapolate the best available DEM to fill the masked waterbody, considering gradual change in the Earth’s geometry. For our methods, even an approximate bathymetry can be sufficient, which is later updated (for the range of observed heights during the study).

These methods have inherent unknown uncertainties because of the inaccuracy of the DEM and errors in the measurement of satellite observations (discussed later). For the ungauged reservoirs where there is no ground information to validate, it is difficult to say which is the more accurate estimation of reservoir volume. Therefore, we propose combining the absolute volumes from ABV and LBV using a Kalman-filter-based state space model (SSM). The state space model or state space time series analysis (as applied in this work) has applications in many different fields [21,22,23,24]. Wallerman et al. [25] presented a Bayesian state space model of forest attributes using field measurements and remote sensing data. In this study, we present a Kalman-filter-based state space model for remote sensing data assimilation and time series analysis. SSM is robustly capable of dealing with the trends, seasonality, interventions and uncertainties in a time series. In this study, the SSM has also been used to fill the missing values in ABV and LBV and to predict two years’ water volume, based on the time series of combined estimates. Besides, as a by-product, we generated the latest bathymetry for the receded reservoirs, which could be used for estimating future refilling and variations. The methodology of this paper is described in Figure 1; within this paper, we use the word “reservoir” for all water bodies including both artificial and natural lakes.

To validate our algorithms, we estimated the absolute volume and volumetric variations of two drastically declining water bodies: the Aral Sea (an ungauged/poorly gauged lake) and Lake Mead (a very well monitored and recorded reservoir). The Aral Sea was the fourth largest lake in the world until the 1960s [26]; however, it has lost 85% of its area and 92% of its water volume in the past half century (1960s–2012) [27]. Moreover, the population dependent on the Aral Sea has increased by four times [28]. Lake Mead is the largest capacity reservoir in the United States and dropped to 39% of its capacity in July 2014 (Bureau of Reclamation). The topography of the two test sites is quite different; the Aral Sea is located in the central Asian plains whereas Lake Mead is in the foothills of the Rocky Mountains, enabling analysis of our methods in different scenarios.

Figure 1. Methodology of the paper.

Figure 1. Methodology of the paper.

2. Data and Pre-Processing

2.1. Changes in the Shoreline

The monthly time-series of the Landsat TM/ETM+/OLI images, from January 2002 to December 2014 are downloaded from www.earthexplorer.usgs.gov. To cover the entire Aral Sea, a minimum of four sets of Landsat images are required, while Lake Mead is covered by only one Landsat image for each observation. Furthermore, in the Aral Sea region almost no interpretable images are available for 3–4 months during winter (from December to March) because of ice, snow and storms, while Lake Mead does not have such climatic conditions. For the Aral Sea, only 84 mosaicked sets can be produced from a combination of the all Landsat satellites during the study period, while Lake Mead is over sampled by 250 snap shots for the same duration. To produce a monthly time series, data gaps are filled by a SSM-based missing data modification method (described in Section 3.2.3 and Section 4.3.1) and oversampled data within a month are averaged.

Additional pre-processing is required for Landsat ETM+ images from 31 May 2003 because of a scan line corrector (SCL) failure in the sensor, which caused a 22% loss in any given scene. Scaramuzza et al. [29] demonstrated that for a single image, local linear histogram matching provides a good result for SLC-off gap filling, which employs a triangulation technique in a moving window of each missing pixel. Therefore, for Lake Mead, the Delaunay triangulation method is used to fill the data gap by calculating triangles from surrounding good pixels. In the case of the Aral Sea, the Delaunay triangulation filling method can introduce classification errors because of the adjoining very bright pixels from the salt deposits of the shrinking lake. Therefore, two extra Landsat scenes are mosaicked to fill the SLC-off stripes from adjacent paths. Hence, a total of six scenes for each timeframe is mosaicked for the Aral Sea, which were acquired on different dates within a month. The overlapping areas are color matched by histogram matching to obtain a seamless mosaic. The images are geo-rectified by the nearest-neighbor resampling method to co-register them during the mosaic operation. The resultant image still had some residual gaps. Therefore, the mosaicked scenes with an SLC-off gap are first classified, then filled by focal filling (3–5 times). Consequently, the mask size increased; therefore, after complete filling, one border pixel is deleted per focal filling instance to reduce the mask to approximately its original size and shape.

The mosaicked near-infrared band of Landsat images are classified by an unsupervised isodata classifier to obtain the water mask at a 30 m spatial resolution. The Aral Sea has a complex spectral signature due to salt deposits and shallow water; therefore it is first classified into 10 spectral clusters, which are later merged into two, as water and non-water. First, two spectral clusters are mostly identified as water bodies by visual interpretation; the remaining classes are set as land. In the shallow East Aral Sea, the water boundary is not very clear, so that even the visual interpretation is uncertain. Unfortunately, we have no ground data from the Aral Sea to verify the classification results. For Lake Mead, direct two-class isodata classification was sufficient because of its deep structure and clear water boundary. The shallow East Aral Sea disappeared completely during the study period, while relatively little change is observed in the size and shape of Lake Mead (Figure 2).

Figure 2. Receding shorelines observed from Landsat. Left: the Aral Sea, and right: Lake Mead.

Figure 2. Receding shorelines observed from Landsat. Left: the Aral Sea, and right: Lake Mead.

2.2. Changes in Water Level

In this study, observations of satellite altimetry are taken from high frequency data of Envisat, Jason1, Jason2, Jason1 extended and Saral/Altika mission, using all passes and missions available at coordinates shown in Table 1. The time series of water levels have been derived from DGFI’s data-base for inland waters (DAHITI) at http://dahiti.dgfi.tum.de/en/ [5].

Table 1. Details of used altimetry missions.

Reservoir	Longitude (Degrees)	Latitude (Degrees)	Mission	Pass
Mead	−114.49	36.137	Envisat	0811
			Jason1 extended	180
			Saral/Altika	0811
North Aral	60.7489	46.5211	Envisat	0126, 0167, 0625
			Jason1	218, 107
			Jason2	218, 107
			Saral/Altika	0126, 0167, 0625
East Aral	59.7146	44.9415	Jason1	142, 107
			Jason2	142, 107
			Saral/Altika	0670
West Aral	58.5626	45.1947	Envisat	0797, 0212
			Jason1 extended	107
			Jason2	142
			Saral/Altika	0797, 0212

Table 1. Details of used altimetry missions.

Reservoir	Longitude (Degrees)	Latitude (Degrees)	Mission	Pass
Mead	−114.49	36.137	Envisat	0811
			Jason1 extended	180
			Saral/Altika	0811
North Aral	60.7489	46.5211	Envisat	0126, 0167, 0625
			Jason1	218, 107
			Jason2	218, 107
			Saral/Altika	0126, 0167, 0625
East Aral	59.7146	44.9415	Jason1	142, 107
			Jason2	142, 107
			Saral/Altika	0670
West Aral	58.5626	45.1947	Envisat	0797, 0212
			Jason1 extended	107
			Jason2	142
			Saral/Altika	0797, 0212

Lake Mead in situ elevation data are taken from a local power-house datum at coordinates of 36.02°N, 114.74°W. The power-house datum is 0.762 m above the North American Vertical Datum of 1988 (NAVD-88). Therefore, we added the difference to bring the altimetry observations to the in situ power-house reference system.

2.3. Bathymetry

In this study, we collected DEMs from different sources and modified them to obtain the best possible bathymetry. For the Aral Sea, the digital contour map at 1 m contour spacing, derived from a 1960s 1/500,000 scale map, is provided by Prof. Renard [30] by personal communication. The given map is re-projected to UTM Zone 40 over WGS84 North datum. The poly-lines are converted to a point feature class for interpolation. Universal kriging of the point data, modeled by spherical semi-variogram and re-sampled by bilinear interpolation generates the raster digital bathymetry model at 30 m resolution. The resultant digital bathymetry model showed better agreement with the shoreline time-series generated from Landsat images than the model used in our earlier publication [31].

The 1-m spatial resolution DEM for Lake Mead is generated from the 5 ft. resolution triangulated irregular network (TIN) data provided by the Bureau of Reclamation, projected to UTM Zone 11, NAD83 datum. For security reasons, the Hoover dam storage region was masked out in the data provided. To replace the missing data, another 10 m resolution DEM of Lake Mead was downloaded from USGS database [32] and interpolated to 1 m intervals. The combined 1 m spatial resolution bathymetry is used in the study for Lake Mead.

3. Methodology

3.1. Methods for Reservoir Volume Estimation

3.1.1. Altimetry-Bathymetry-Volume (ABV) Method

Water level observations from satellite altimetry are intersected with the bathymetry of the water body to produce the absolute water volume (ABV) by integrating all water height columns within a reservoir (discussed in Section 4.3.1 for Lake Mead and Section 4.3.2 for the Aral Sea). The resultant absolute volumes are further reduced by the first observation of the time series to obtain volume variation during the study period (discussed in Section 4.2).

3.1.2. Landsat-Bathymetry-Volume (LBV) Method

The Landsat-bathymetry-volume method is based on the assumption that the lake and reservoirs can be considered flat equipotential surfaces. Theoretically, the shorelines of the masked bathymetry image can be considered as one contour line, and all border pixels should have the same height. Extracting the water height from the masked bathymetry is not trivial because in reality, the extracted border pixels have a large range and elevation variance. This can be due to uncertainty in the water mask, coarse bathymetry and high gradients. Taking into account all these issues, we obtain a range of water heights from the shoreline and bathymetry intersection. The algorithm developed to obtain the final height from the Landsat mask is described in Figure 3.

We traced the border of the masked bathymetry clockwise and collected the boundary pixel values to obtain the entire boundary vector (EBV), starting from the first column of the mask from the left (Figure 4, red asterisks). At the same time, we selected a subsection of the water body (Figure 4, red rectangles) based on the knowledge of the study site to obtain selected region boundary (SRB) values. We proposed two criteria for the selection of the subsection of the basin. (1) The subsection should have a comparatively lower gradient, so that there is less height variation within each 30 m pixel; (2) It should be comparatively stable, i.e., not under direct flow of water (Figure 4: red magenta arrow shows the river mouth) to uphold the concept of an equipotential surface.

Figure 3. Algorithm of water level estimation in the LBV method: rectangular boxes indicate processes, trapezoids indicate the manual operation and multiple document boxes indicate the time series of the images.

Figure 3. Algorithm of water level estimation in the LBV method: rectangular boxes indicate processes, trapezoids indicate the manual operation and multiple document boxes indicate the time series of the images.

Figure 4. Contour of reservoirs; * starting point of boundary pixel collection; the red rectangle is the subsection of the reservoirs for selected area boundary pixels and the large magenta dashed arrow shows the mouth of the river.

Figure 4. Contour of reservoirs; * starting point of boundary pixel collection; the red rectangle is the subsection of the reservoirs for selected area boundary pixels and the large magenta dashed arrow shows the mouth of the river.

As examples, Landsat image shorelines for every fifth available epoch from the East Aral Sea are plotted in Figure 5. The plot shows the elevation variation along the shoreline; the blue lines on the left axis are EBV, the red lines are the pixels from SRB; and the green plot on the right axis shows the standard deviation of the 50 adjacent pixels of EBV. The 50-pixel moving window standard deviation indicates that there are regions in the reservoir that always have relatively low elevation variance. For example, in the East Aral Sea, less than 10 cm of dispersion is observed from the eastern and southern ends (Figure 4). The beginning and end of EBV have relatively high standard deviations (Figure 5) because of close contours in the western and northern sides of the East Aral Sea (Figure 4). This demonstrates that the area selected for SRB is relatively stable.

Figure 5. (Left) Boundary pixels of the East Aral Sea obtained from Landsat mask at different time points: blue line is the entire boundary vector (EBV), red portions are pixels from the selected region boundary (SRB); (Right) 50-pixel moving standard deviation of EBV.

Figure 5. (Left) Boundary pixels of the East Aral Sea obtained from Landsat mask at different time points: blue line is the entire boundary vector (EBV), red portions are pixels from the selected region boundary (SRB); (Right) 50-pixel moving standard deviation of EBV.

In addition, we rejected pixels from the southern part of the East Aral Sea, because it is at the river mouth, and during flood events, it would not be an equipotential surface. We can observe this phenomenon in the subplots of Figure 5, e.g., floods of August 2010 and April 2013. A wide range of water heights are observed during the flood events from the southern part (Figure 5), although a 50-pixel standard deviation is still low. Therefore, EBV estimates higher water height during the flood events than SRB, because of the violation of equipotential surface assumption (discussed in Section 4.1). The total number of pixels along the shoreline are also significantly reduced during the study period because of the shrinking of the East Aral Sea (Figure 5). A similar phenomenon is observed from all other reservoirs; therefore, pixels from the SRB regions shown as red rectangles in Figure 4 are processed to obtain the water height from all reservoirs.

The outliers from the SRB vectors (Figure 5: red line) are rejected on the box plot distributional assessment basis [33,34]. To find the outliers, a lower quantile (Q1), upper quantile (Q3) and inter-quantile range (H) are calculated for each SRB vector. Values lower than (Q1 − Q3) × H and higher than (Q3 + Q3) × H are considered outliers and hence rejected. The final water height for each shoreline is assigned by taking the mean of the remaining elevation vector.

Finally, the absolute water volume time series (LBV) is generated for each reservoir from the SRB elevations by integrating the entire water height column within the shoreline (discussed in Section 4.3.1 for Lake Mead and Section 4.3.2 for the Aral Sea). Similar to ABV, LBV time series are also reduced by the first observation to obtain volume variations (discussed in Section 4.2).

3.1.3. Altimetry-Landsat Volume Variation (ALVV) Method

In many cases, bathymetry of a water body is not available. Here a completely remote sensing based approach can be the best alternative for volume variation estimation. Unlike the first two methods, instead of using a DEM, we combined the area calculated from Landsat and the elevation calculated from altimetry to produce the volumetric changes in the water body using a pyramid frustum volume estimation method. This method assumes the water body to be an inverted truncated pyramid [18]. It computes the change in volume at each time step with the change in the elevation and area of the truncated inverted pyramid with respect to the previous time step. The cumulative summation of each time step volume change derived from the trapezoidal volume formula is given in Equation (1).

$$ALVV={\displaystyle \sum }_{t=1}^n\frac{1}{3}\times \left(H_t-H_{t-1}\right)\times \left(A_t+A_{t-1}+\sqrt{\left(A_t\times A_{t-1}\right)}\right)$$

(1)

where

$ALVV$ = Total volume change with respect to the initial state (t₀) at the nth month

$A_t$ = Area of the water extent at month t and $A_{t-1}$ = Area of the water extent at the previous month

$H_t$ = Elevation of the water body at month t and $H_{t-1}$ = Elevation of the water body at the previous month. n = Number of months.

The two observation time series of the water surface area and water level are often difficult to synchronize because the dates of acquisition of data from Landsat and different altimetry missions are often not the same. For Lake Mead, area and elevation data points that fall within the same 10-day period are used. For the Aral Sea, a monthly window is used because the whole lake is covered by combining three Landsat paths; therefore, a single date cannot be assigned per observation. The outliers are rejected on the basis of the quantile quantile plot [35] and then from the extracted area and water level samples, a linear relationship is derived (Figure 6). The coefficient of determination (R²) between the area and elevation shows a strong linear correlation (>0.92) with NRMSE <2% for all of the reservoirs, except the East Aral Sea (NRMSE = 6%), because of the large area-level ratio due to its almost flat structure. Therefore, the derived relationship can be used to fill the missing values in any of the two sets of data. In the case of unavailability of both datasets, the time series is interpolated based on nearest neighbors. The resultant continuous monthly area level data are used in Equation (1) to estimate storage variations (discussed in Section 4.2).

Figure 6. Surface area and water height relation.

Figure 6. Surface area and water height relation.

3.2. State Space Model for the Estimation and Prediction of Absolute Volume Time Series

3.2.1. Structure of a State Space Model (SSM)

State space modeling provides a unified methodology for treating a wide range of problems in time series analysis such as forecasting, signal extraction and estimation of parameters. It assumes that an unobserved series of vectors, in our case true volume $x$, is associated with a series of observations (ABV and LBV) y₁ and y₂, which determine the development of the system over time. The relation between $x$_t and y_t is specified by the state space model. The state-space model or dynamic linear model (DLM), in its basic form, employs an order one vector auto regression as the state equation (Equation (2)) [36].

$$x_t=\mathsf{\varphi }{x}_{t-1}+w_t$$

(2)

where the state equation determines the rule for the generation of the $p\times 1$ system state vector $x_t$ from the past $p\times 1$ state $x_{t-1}$ for time points $t=1,\mathrm{...},n$. We assume that the, process noise $w_t$ is $p\times 1$ independent and identically distributed, zero-mean normal vectors with covariance matrix $Q$ (i.i.d. N (0, $Q$)) and $\mathsf{\varphi }$ is a state transition matrix. The process noise represents some sudden change in the trend, e.g., a flood in our case. In the DLM, we assume the process starts with a normal vector $x_0$ that has mean ${\mathsf{\mu }}_0$ and $p\times p$ error covariance matrix${\mathsf{\sigma }}_0$.

Further, the SSM assumes that we observe the state vector $x_t$ indirectly by its linear transformed version with added noise (Equation (3)) i.e., unknown observational error associated with ABV and LBV in this case.

$$y_t=A_t{x}_t+v_t$$

(3)

where $A_t$ is a $q\times p$ measurement matrix and $v_t$ is a measurement noise vector ~i.i.d. N (0, R) with measurement noise covariance matrix R. Equation (4) is called the observation or measurement equation.

The primary aim of any analysis involving the state-space model would be to produce estimators for the underlying unobserved signal $x_t$ given the data $Y_s=\left\{y_1,y_2\right\}$, to time $s$ by minimizing error covariance. When $s<t$, the problem is called forecasting or prediction; when $s=t$, the problem is called filtering; and when$s>t$, the problem is called smoothing, which is used for filling the data gap. The solution to these problems is accomplished via the Kalman filter and smoother [36]. Throughout the discussion, we will use the following definitions:

$$\text{Estimation }{x}_t^s=E\left(x_t|Y_s\right)$$

(4)

and

$$\text{Error covariance }{P}_{t_1,t_2}^s=E\{(x_{t_1}-x_{t_1}^s){(x_{t_2}-x_{t_2}^s)}^{\prime }\}$$

(5)

The Kalman filter minimizes the state error covariance matrix P by weighted least squares using Q (Equation (8)) and R (Equation (12)). Therefore, the Kalman filter can be considered as a predictor-corrector system.

3.2.2. The Kalman Filter

For the state space model specified in Equations (2) and (3), state propagation and state error propagation is defined by Equations (6) and (7) for $t=1,\dots , n$

$$x_t^{t-1}=\mathsf{\varphi }{x}_{t-1}^{t-1}+w_t$$

(6)

$$P_t^{t-1}=\mathsf{\varphi }{P}_{t-1}^{t-1}{\mathsf{\varphi }}^{\prime }+Q$$

(7)

The prior state estimate (Equation (6)) and the prior error covariance (Equation (7)) are updated based on observations $y_t$ to obtain the posterior state estimate (Equation (8)) and the posterior error covariance (Equation (9)). The prediction for $t>n$ is accomplished via Equations (6) and (7) with initial conditions$x_n^n|P_n^n$.

$$x_t^t=x_t^{t-1}+K_t\left(y_t-A_t{x}_t^{t-1}\right)$$

(8)

$$P_t^t=\left[I-K_t{A}_t\right]P_t^{t-1}$$

(9)

where

$$K_t=P_t^{t-1}{A}_t^{\prime }{\left[A_t{P}_t^{t-1}{A}_t^{\prime }+R\right]}^{-1}$$

(10)

The Kalman gain $K_t$ is computed by a recursive least squares solution, which trades off state uncertainty and measurement noise. Equations (8) and (10) suggest that when the magnitude of R is small, the measurements are accurate, and then the state estimate depends mostly on the measurements and vice versa. The new measurements add to knowledge of the system; therefore, the estimated error decreases in Equation (9). The posterior uncertainty or covariance is always less than the prior as it is reduced by the multiple of the Kalman gain (Equation (10)). Important byproducts of the filter are the innovations (prediction errors)

$${ϵ}_t=y_t-E\left(y_t\text{|}{Y}_{t-1}\right)=y_t-A_t{x}_t^{t-1}$$

(11)

and the corresponding variance-covariance matrices

$$var\left({ϵ}_t\right)=var[A_t\left(x_t-x_t^{t-1})+v_t\right]=A_t{P}_t^{t-1}{A}_t^{\prime }+R$$

(12)

The proof of these Kalman filter equations can be found in Shumway and Stoffer, (2011) [37].

3.2.3. Missing Data Modifications

An attractive feature available within the state-space framework is its ability to treat time series that have been observed irregularly over time. The definition of smoothing described above is used to obtain the estimate for missing data [37].

3.2.4. Estimates of Underlying Signal from Two Observation Time Series

In the present study, the absolute water volume estimates have been obtained from two observation time series LBV and ABV. We suppose that both series observe the same signal of actual water volume of the reservoir $x_t$ with different error in the measurements ($v_{t1}\text{ and} v_{t2}$). In our case the state transition model $\mathsf{\varphi }$ and measurement model $A_t$ are the identity matrix.

$$y_{t1}=x_t+v_{t1}\text{ and }{y}_{t2}=x_t+v_{t2}$$

(13)

more compactly,

$$\left(\begin{array}{c}{y}_{t1}\\ y_{t2}\end{array}\right)=\left(\begin{array}{c}1\\ 1\end{array}\right)x_t+\left(\begin{array}{c}{v}_{t1}\\ v_{t2}\end{array}\right)$$

(14)

where

$$R=var\left(\begin{array}{c}{v}_{t1}\\ v_{t2}\end{array}\right)=\left(\begin{array}{cc}{r}_{11}& 0\\ 0& r_{22}\end{array}\right)$$

(15)

Here r₁₁ and r₂₂ are the variance of the observation error $v_{t1}$ and $v_{t2}$ respectively, which is equal to one at t = 1. In our case, we do not know the noise in ABV and LBV. Therefore at t = 1 we randomly put some value for $v_{t1}$ and $v_{t2}$, which is recursively optimized to estimate error by Equations (11) and (12). The Combined SSM estimate (CSSME) takes the mean of LBV and ABV reduced by their respective error (Equation (16)).

$$x_t=1/2\left\{\left(y_{t1}-v_{t1}\right)+\left(y_{t2}-v_{t2}\right)\right\}$$

(16)

It is reasonable to suppose that the unknown common signal $x_t$, can be modeled as a random walk with drift (Equation (17)) of the form

$$x_t=\mathsf{\delta }+x_{t-1}+w_t$$

(17)

with $\mathsf{\delta }$ = stochastic trend component of the signal, computed on every step by maximum likelihood estimator. The estimation of unknown parameters contained in $Q \text{and }R$, which define this model, is carried out using the Kalman recursion as described above.

3.2.5. Signal Extraction and Forecasting

Once the CSSME is obtained from ABV and LBV time series data by applying the aforementioned procedure, the recursive six monthly forecast for the next two years has been carried out based on CSSME time series. Since CSSME exhibits trend as well as seasonality, we consider the series to be the sum of a trend component, a seasonal component, and white noise. The CSSME series can be expressed as the Equation (18) [37]

$$y_t=T_t+S_t+v_t$$

(18)

where $T_t$ is trend and $S_t$ is the seasonal component. Suppose we allow the trend to increase exponentially; that is,

$$T_t=\mathsf{\varphi }{T}_{t-1}+w_{t1}$$

(19)

where the coefficient $\mathsf{\varphi }>1$ characterizes the increase and $\mathsf{\varphi }<1$ charaterizes the decrease. Let the seasonal component be modeled as

$$S_t+S_{t-1}+S_{t-2}+S_{t-3}=w_{t2}$$

(20)

which corresponds to assuming the seasonal component is expected to sum to zero over a complete period or four quarters. $w_{t1} \text{and} w_{t2}$ are the trend and seasonal noise respectively.

The model can be expressed in state-space form as below. Let $x_t=\left(T_t,S_t,S_{t-1},S_{t-2}\right)\prime$ be the state vector, then the state equation (Equation (2)) can be written as (Equation (21))

$$\left[\begin{array}{c}{T}_t\\ S_t\\ S_{t-1}\\ S_{t-2}\end{array}\right]=\left[\begin{array}{cccc}\mathsf{\varphi }& 0& 0& 0\\ 0& -1& -1& -1\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{T}_{t-1}\\ S_{t-1}\\ S_{t-2}\\ S_{t-3}\end{array}\right]+\left[\begin{array}{c}{w}_{t1}\\ w_{t2}\\ 0\\ 0\end{array}\right]$$

(21)

and the measurement equation (Equation (3)) as

$$y_t=\left[1 1 0 0\right]\left[\begin{array}{c}{T}_t\\ S_t\\ S_{t-1}\\ S_{t-2}\end{array}\right]+v_t$$

(22)

with $R=r_{11}$ and

$$Q=\left[\begin{array}{cccc}{q}_{11}& 0& 0& 0\\ 0& q_{22}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$$

(23)

The model reduces to state space form, (Equations (2) and (3)), with $p=4$ and$q=1$. The parameters to be estimated by maximum likelihood for the model are$$the noise variance in the measurement equations$(r_{11}$), the model variances corresponding to the trend and seasonal components ($q_{11}$ and$q_{22}$, respectively), and$$the transition parameter $\mathsf{\varphi }$ that models the increasing or decreasing trend. This model is used for forecasting purposes. Most of the equations presented in the SSM section are directly adapted from Shumway and Stoffer (2011) [37].

4. Results

4.1. Water Height

Lake Mead is often observed twice per year by Landsat and therefore the oversample height measurements were averaged per month. On the other hand, the lake has no altimetry observations for one year in 2012 due to an intermission gap. There, the Landsat-SRB height can be a good substitute. Furthermore, the altimetry measurements obtained from the East Aral Sea for 2008 are contaminated by the reflection from land due to shrinking of the reservoir. Later available observations for 2010 and 2011 from altimetry have overestimated the water height due to flood and land boundary contamination as altimetry tracks cross the East Aral Sea from its southern end. These effects can be noticed in Figure 7 (top right), the Landsat-EBV (green line) and altimetry (blue line) observed false high, because during the flood event, the assumption of an equipotential surface get abandoned. Consequently, water heights observed from the southern part of East Aral Sea are higher than the relatively stable part selected in the Landsat-SRB (red line). Therefore, for the East Aral Sea from 2008, only Landsat-SRB water heights are used, and meaningful comparisons of the models for the reservoir are restricted until December 2007.

Figure 7. (Left) Water height observed by altimetry, Landsat: entire boundary vector (EBV), selected region boundary (SEB) and in situ; (Right) error bar for Landsat-SRB water height.

Figure 7. (Left) Water height observed by altimetry, Landsat: entire boundary vector (EBV), selected region boundary (SEB) and in situ; (Right) error bar for Landsat-SRB water height.

The water heights derived by the Landsat-SRB (Figure 7: red dashed line) from 2D Landsat masks are in very good agreement with the observed heights from satellite altimetry for all of the reservoirs. Landsat-EBV (Figure 7: green line) has measured higher elevations particularly during extreme events; for example, during the extreme droughts of the West Aral Sea (2014) and Lake Mead (2010 and 2014), and during floods of the East Aral Sea (2010 and 2013).

The approximate decrease in water height observed in Lake Mead, the East Aral Sea and the West Aral Sea during 2002–2014 are 28 m, 6 m and 6.5 m respectively, while a 2 m gain is observed for the North Aral Sea (Figure 7).

4.2. Volumetric Variation

To compare volumetric variations from ALVV, absolute volumes from ABV and LBV are subtracted by the first observation of each time series (Figure 8), which changes with respect to the January 2002 observation. The ALVV method can be applied to the reservoirs where no ground observations of the structure of the lakebed are available; thus, the volume cannot be calculated from ABV and LBV. These methods are preliminarily tested for Lake Mead, which has daily in situ data to validate the results, and the algorithms are then applied to the three sub-parts of the Aral Sea. The reduced ABV and LBV has a very good agreement with the in situ observations of Lake Mead and between each other in all of the sub-parts of the Aral Sea (discussed in Section 4.4.2). The North Aral Sea results are statistically poor but qualitatively the results are within the acceptable range (discussed in Section 4.4.1).

Figure 8. Volumetric variation computed by ABV, LBV and ALVV methods for the three sub-parts of the Aral Sea and for Lake Mead.

Figure 8. Volumetric variation computed by ABV, LBV and ALVV methods for the three sub-parts of the Aral Sea and for Lake Mead.

4.3. SSM Estimation and Prediction

4.3.1. Standard Case of Lake Mead

Lake Mead has been considered as a standard case to demonstrate the effectiveness of SSM because in situ data are available for the reservoir to verify the results. The SSM estimations for Lake Mead are plotted in Figure 9 and Figure 10; a similar process is applied for the other reservoirs of the Aral Sea and their results are shown in Figure 11.

Figure 9. Missing data modification by SSM for the Lake Mead. (Top left) Landsat-bathymetry-volume (LBV) estimation and data gap filled by SSM; (Bottom left) difference between LBV-SSM estimates and in situ observations; (Top right) altimetry-bathymetry-volume (ABV) estimation and data gap filled by SSM; (Bottom right) difference between ABV-SSM estimates and in situ observations.

Figure 9. Missing data modification by SSM for the Lake Mead. (Top left) Landsat-bathymetry-volume (LBV) estimation and data gap filled by SSM; (Bottom left) difference between LBV-SSM estimates and in situ observations; (Top right) altimetry-bathymetry-volume (ABV) estimation and data gap filled by SSM; (Bottom right) difference between ABV-SSM estimates and in situ observations.

Figure 10. Lake Mead SSM analysis. (Top left) The combined SSM estimate (CSSME) (magenta line) and the forecast (green line) for 2013 and 2014; (Bottom left) difference between CSSME and in situ observations; (Top right) estimated seasonal component; (Bottom right) estimated trend component. The dashed cyan lines indicate the upper and lower 95% confidence limit.

Figure 10. Lake Mead SSM analysis. (Top left) The combined SSM estimate (CSSME) (magenta line) and the forecast (green line) for 2013 and 2014; (Bottom left) difference between CSSME and in situ observations; (Top right) estimated seasonal component; (Bottom right) estimated trend component. The dashed cyan lines indicate the upper and lower 95% confidence limit.

Figure 11. Aral Sea absolute volume (ABV and LBV), Combined SSM estimates, and SSM forecast.

Figure 11. Aral Sea absolute volume (ABV and LBV), Combined SSM estimates, and SSM forecast.

Filling up the Missing Data with SSM

Since our datasets (ABV and LBV) have many missing values from all the reservoirs, these gaps are filled dynamically by respective ABV-SSM and LBV-SSM estimates (Figure 9). Consistent positive offset is observed in ABV and also frequently in LBV (Figure 9: bottom plot), which might be due to overestimation of height by altimetry and Landsat-SRB because of the difference in the geoid of the sites of measurement. The in situ measurements are taken from the western side of Lake Mead near Hoover dam, and altimetry and Landsat-SRB observations are from Overton Arm of Lake Mead at north. However, LBV offset has a slightly lesser magnitude than ABV. Therefore, the ABV-SSM has a relatively high RMSE with respect to in situ observations, relative to the LBV-SSM (Table 2).

Table 2. Comparison of results from all the reservoirs by root mean square error (RMSE) and coefficient of determination (R²).

Reservoir	Compared Results		RMSE (NRMSE)	R2
Lake Mead	Landsat-SRB water height	In situ water height	0.59 m (0.2%)	0.99
	Altimetry water height	In situ water height	0.67 m (0.2%)	0.99
	Landsat bathymetry volume (LBV)	In situ absolute volume	0.32 km3 (1.6%)	0.98
	Altimetry bathymetry volume (ABV)	In situ absolute volume	0.41 km3 (2.1%)	0.97
	Altimetry Landsat volume variation (ALVV)	In situ reduced by first obs. (volume variation)	0.53 km3	0.96
	ALVV	LBV reduced by first obs.	0.31 km3	0.98
	ALVV	ABV reduced by first obs.	0.06 km3	0.99
	Combined SSM estimate (CSSME)	In situ absolute volume	0.35 km3 (1.8%)	0.97
	CSSME	LBV	0.32 km3 (1.6%)	0.98
	CSSME	ABV	0.41 km3 (2.1%)	0.97
	CSSME Forecast	CSSME	0.53 km3 (3.0%)	0.80
	CSSME Forecast	In situ absolute volume	0.66 km3 (3.7%)	0.75
West Aral Sea	Landsat-SRB water height	Altimetry water height	0.44 m (1.5%)	0.94
	LBV	ABV	1.58 km3 (4.5%)	0.94
	ALVV	LBV reduced by first obs.	1.73 km3	0.94
	ALVV	ABV reduced by first obs.	1.09 km3	0.97
	CSSME	LBV	1.00 km3 (2.7%)	0.97
	CSSME	ABV	0.50 km3 (1.6%)	0.99
	CSSME Forecast	CSSME	0.52 km3 (1.9%)	0.76
North Aral Sea	Landsat-SRB water height	Altimetry water height	0.50 m (1.2%)	0.32
	LBV	ABV	1.59 km3 (6.0%)	0.23
	ALVV	LBV reduced by first obs.	1.46 km3	0.56
	ALVV	ABV reduced by first obs.	0.15 km3	0.99
	CSSME	LBV	0.82 km3 (3.2%)	0.83
	CSSME	ABV	0.82 km3 (3.1%)	0.83
	CSSME Forecast	CSSME	1.09 km3 (3.9%)	0.23
East Aral Sea	Landsat-SRB water height (January 2002–December 2007)	Altimetry water height	0.33 m (1.0%)	0.82
	LBV (01.2002–12.2007)	ABV (January 2002–December 2007)	2.40 km3 (13%)	0.85
	ALVV	LBV reduced by first obs.	4.26 km3	0.51
	ALVV	ABV reduced by first obs.	0.67 km3	0.99
	CSSME (January 2002–December 2007)	LBV (January 2002–December 2007)	1.92 km3 (12.0%)	0.89
	CSSME (January 2002–December 2007)	ABV (January 2002–December 2007)	0.49 km3 (2.7%)	0.99
Entire Aral Sea	LBV	ABV	4.10 km3 (5.5%)	0.91
	CSSME Forecast	CSSME	1.40 km3 (2.4%)	0.60

Table 2. Comparison of results from all the reservoirs by root mean square error (RMSE) and coefficient of determination (R²).

Reservoir	Compared Results		RMSE (NRMSE)	R2
Lake Mead	Landsat-SRB water height	In situ water height	0.59 m (0.2%)	0.99
	Altimetry water height	In situ water height	0.67 m (0.2%)	0.99
	Landsat bathymetry volume (LBV)	In situ absolute volume	0.32 km3 (1.6%)	0.98
	Altimetry bathymetry volume (ABV)	In situ absolute volume	0.41 km3 (2.1%)	0.97
	Altimetry Landsat volume variation (ALVV)	In situ reduced by first obs. (volume variation)	0.53 km3	0.96
	ALVV	LBV reduced by first obs.	0.31 km3	0.98
	ALVV	ABV reduced by first obs.	0.06 km3	0.99
	Combined SSM estimate (CSSME)	In situ absolute volume	0.35 km3 (1.8%)	0.97
	CSSME	LBV	0.32 km3 (1.6%)	0.98
	CSSME	ABV	0.41 km3 (2.1%)	0.97
	CSSME Forecast	CSSME	0.53 km3 (3.0%)	0.80
	CSSME Forecast	In situ absolute volume	0.66 km3 (3.7%)	0.75
West Aral Sea	Landsat-SRB water height	Altimetry water height	0.44 m (1.5%)	0.94
	LBV	ABV	1.58 km3 (4.5%)	0.94
	ALVV	LBV reduced by first obs.	1.73 km3	0.94
	ALVV	ABV reduced by first obs.	1.09 km3	0.97
	CSSME	LBV	1.00 km3 (2.7%)	0.97
	CSSME	ABV	0.50 km3 (1.6%)	0.99
	CSSME Forecast	CSSME	0.52 km3 (1.9%)	0.76
North Aral Sea	Landsat-SRB water height	Altimetry water height	0.50 m (1.2%)	0.32
	LBV	ABV	1.59 km3 (6.0%)	0.23
	ALVV	LBV reduced by first obs.	1.46 km3	0.56
	ALVV	ABV reduced by first obs.	0.15 km3	0.99
	CSSME	LBV	0.82 km3 (3.2%)	0.83
	CSSME	ABV	0.82 km3 (3.1%)	0.83
	CSSME Forecast	CSSME	1.09 km3 (3.9%)	0.23
East Aral Sea	Landsat-SRB water height (January 2002–December 2007)	Altimetry water height	0.33 m (1.0%)	0.82
	LBV (01.2002–12.2007)	ABV (January 2002–December 2007)	2.40 km3 (13%)	0.85
	ALVV	LBV reduced by first obs.	4.26 km3	0.51
	ALVV	ABV reduced by first obs.	0.67 km3	0.99
	CSSME (January 2002–December 2007)	LBV (January 2002–December 2007)	1.92 km3 (12.0%)	0.89
	CSSME (January 2002–December 2007)	ABV (January 2002–December 2007)	0.49 km3 (2.7%)	0.99
Entire Aral Sea	LBV	ABV	4.10 km3 (5.5%)	0.91
	CSSME Forecast	CSSME	1.40 km3 (2.4%)	0.60

Combined SSM Estimate (CSSME) and Forecast

We generated SSM estimates for all the reservoirs to assimilate the calculated remote sensing product from Landsat and altimetry. All the series have a clear seasonal cycle of six months and a trend signal, which gradually decreases over time, except for the North Aral Sea. The trend and seasonal components of the time series are modeled according to Equations (19) and (20), and the results are shown in the right panel of Figure 10. The model CSSME has also been used to make six month recursive forecasts for the next two years. The CSSME is extrapolated first for the 6 months and then again for the next 6 months and so on, to get forecasts from January 2013 until December 2014 for all of the reservoirs, except the East Aral Sea.

The Lake Mead trend shows a drop of approximately 1 km³ per year with a maximum seasonality of 1 km³ except during the 2011 flood, when the lake gained almost 8 km³ water. The CSSME for 2012 is substituted by the LBV-SSM because of the lack of altimetry observations for this year; nevertheless, from March 2013 the CSSME series (Figure 10: black line) is again updated by the combined ABV and LBV observations. The forecast (Figure 10: green line) fits very well with the CSSME estimate (3% NRMSE) and with the in situ (3.7% NRMSE).

4.3.2. CSSME and Forecast for Aral Sea

A similar process of filling in missing data, assimilating ABV and LBV to obtain CSSME and model forecast is applied in all the sub-parts of the Aral Sea. The West Aral Sea decreased in volume by approximately 20 km³, showing a dominant trend component. The CSSME for the West Aral Sea showed good agreement of estimate with its both input ABV and LBV (R² > 0.97 and NRMSE < 2.7%) and with its forecast (0.76 R² and 1.9% NRMSE). On the other hand, the North Aral Sea has many jumps in the data and a relatively strong seasonal signal.

Because of shrinking of the East Aral Sea, the ABV series has many missing observations after December 2008. Therefore, CSSME is estimated until December 2008 and forecast is done only for one year (2009) because there was no CSSME to do recursive forecast. The LBV-SSM estimate is plotted for the January 2008–December 2015 duration in the Figure 11 for the East Aral Sea (the lake ultimately dried up in 2014). The entire Aral Sea plot in Figure 11 is an integration of its three parts. The largest difference between ABV and LBV is observed between 2010 and 2012 because of the East Aral Sea shrinking and overestimation by altimetry due to flooding (discussed in Section 4.1).

4.4. Validation

To validate our results R-squared values (which indicate a relative measure of fit) and root mean square error (RMSE) (which gives the absolute fit of one time series with respect to other time series) are estimated (Table 2). Additionally, NRMSE is computed by normalizing the RMSE against the mean of the measurements. However, NRMSE can be misleading if the variable contains both negative and positive values; therefore, volumetric variations (ALVV) are not compared in %.

4.4.1. Validation of Water Height

The water height obtained from the developed Landsat-SRB method fits perfectly well with the in situ and altimetry observations for all the study sites. Lake Mead has a coefficient of determination of 0.99, which means 99% variability of the in situ data can be explained by variability in the Landsat-SRB and altimetry observations (Table 2). The obtained 59 cm RMSE or 0.17% NRMSE show the relative closeness of the Landsat-SRB observed water height to the ground observations. Similarly, the water height observed by altimetry for Lake Mead is also highly accurate, with marginal (6 cm) difference in RMSE from the Landsat-SRB estimates.

Unfortunately, for the Aral Sea reservoirs very few in situ observations are available. They are mostly collected from publications [27,38], which are given without height reference information. Therefore we cross-validated the observations due to lack of possibility of validating with respective ground observations. Altimetry (Figure 7: blue line) and Landsat-SRB (Figure 7: red dashed line) observations follow the same trend and are quite close to the limited available in situ points. The R² values for the two observations for the West Aral Sea and the East Aral Sea (January 2002–December 2007) are 0.94 and 0.82 respectively (Table 2), which shows that the two data have a very strong linear relationship for these reservoirs.

The North Aral Sea has more dispersion between the two missions’ observations; nevertheless, they have good covariance of approximately 0.94. The maximum difference between the water height observed by Landsat and altimetry is less than 90 cm and RMSE is 50 cm, which is well within the error range of a typical Landsat height estimate (Figure 7). However, because of the low variation in the North Aral Sea (approximately 2 m) within the study period, the 90 cm difference has a significant impact on the statistical evaluation (Table 2). The F-distribution test is applied to evaluate variations among the methods, which shows no significant difference in the mean height of all the water bodies except for the North Aral Sea. The F-test rejected the null hypothesis for the North Aral Sea, failing the assumption that the two datasets (i.e., Landsat-SRB and altimetry water heights, and eventually their volumes LBV and ABV respectively) came from normal distributions with the same variance.

4.4.2. Validation of the Estimated Water Volume (LBV, ABV and ALVV)

The Table 2 shows that absolute water volumes estimated by LBV and ABV have a very high correlation with the ground observations for Lake Mead (R²_LBV = 0.98 with RMSE_LBV = 0.32 km³, R²_ABV = 0.97 with RMSE_ABV = 0.41 km³). For an ungauged basin like Aral Sea, cross comparison of different model outputs can be considered as a validation method. The entire Aral Sea and its sub-parts showed a good LBV-ABV agreement (R²_{entire Aral Sea} = 0.91 with 4.1 km³ RMSE). The West Aral Sea has a better agreement between the estimated volumes (0.94 R²) relative to the East Aral Sea (0.85 R²), similar to what we have observed for the estimated water heights. However, despite of the smaller water height RMSE for the East Aral Sea compared to the West Aral Sea, the RMSE of computed volume for the East Aral Sea increased (RMSE_west = 1.5 km³ and RMSE_east2007 = 2.4 km³), which is because of its substantially greater size until 2007.

The volumetric variations estimated by ALVV showed a better agreement with the ABV than LBV for all of the reservoirs (Table 2), simply because in Equation (1) height variations are taken from altimetry and not from the Landsat-SRB. According to Birkett et al. (2011) [39], the accuracy of the estimate should be better than 10% of the expected seasonal amplitude and our all results are well within this range.

4.4.3. Validation of CSSME Volume

For Lake Mead, CSSME shows an extremely high agreement with LBV and ABV (R²_CSSME-LBV = 0.98 and R²_CSSME-ABV = 0.97 with a marginal difference of 0.5% in NRMSE) as well as with the ground observation (R² = 0.97 and RMSE = 0.35 km³) (Table 2). The West Aral Sea CSSME fits slightly better to ABV (R²_CSSME-LBV = 0.97 and R²_CSSME-ABV = 0.99 with 1.1% difference in NRMSE). This shows a very good assimilation of the two input time series in both of the cases. The East Aral Sea CSSME has a higher tilt towards ABV (R²_{CSSME-LBV_2007} = 0.89 and R²_{CSSME-ABV_2007} = 0.99 with 9.3% difference in NRMSE), which implies that the ABV time series is smoother than LBV for the East Aral Sea. The North Aral Sea CSSME agrees almost equally with the LBV and the ABV (R²_CSSME-LBV = R²_CSSME-ABV = 0.83), showing their comparable contribution. The CSSME forecast (for 2012–2014) for all of the reservoirs also showed a very good agreement with the CSSME estimates (Table 2), except for the North Aral Sea, which is because of the discrepancies between the ABV and LBV observations.

4.5. Improved Bathymetry

The Aral Sea is located in an arid region with much sand, and salt deposits in the surrounding barren open land. This gives rise to the possibility of significant displacements into the lakebed of sediments by the action of wind, water, gravity, or snow. The bathymetry of the lake was last observed in the 1960s; therefore, within the last half century, the structure of the lakebed may have changed. As a by-product of our shoreline and water level time series, an updated bathymetric map is generated for the region that fell dry. Figure 12 shows the change in the structure of the lakebed of the East Aral Sea. The river mouth made its impression in the present East Aral Sea, because its shrinking led to the inflow of much sediment into the lake’s interior. In addition, salt deposits along the shorelines increased the corresponding elevation, a phenomenon that is more pronounced in the reduced lakebed because of increased salinity. Note, that height estimates from satellite altimetry were only possible down to a minimum elevation of 27 m above sea level due to a lack of reliable altimetry data over the largely reduced water surface. In order to construct a complete bathymetric chart of the lakebed of the East Aral Sea heights below 27 m were obtained solely from Landsat optical images following our SRB approach. The updated bathymetry map is available as supplement to this paper at http://dx.doi.org/10.1594/PANGAEA.855779 [40].

Figure 12. East Aral Sea. (Left) subset of the 1960s Aral Sea bathymetry; (Right) updated bathymetry using remote sensing data.

Figure 12. East Aral Sea. (Left) subset of the 1960s Aral Sea bathymetry; (Right) updated bathymetry using remote sensing data.

5. Discussion

The combination of multiple remote sensing data is a key for continuous monitoring and measurement of water resources. The size of reservoirs is the biggest limiting factor for the terrestrial application of satellite altimetry for water storage estimation, which can be compensated for by the remote sensing image based approach discussed above. Lake Mead has no altimetry observations for nearly one year (from February 2012 until March 2013), because of an intermission gap between the Jason1 extended and Saral/Altika missions. Similarly, the East Aral Sea has missing altimetry observations for nearly one year (from January 2009 until March 2010) as well as for the remaining study period from August 2011, because of its decreased size and the fact that no mission passed over the reduced Sea in this period [31]. In both cases, the water height estimated from the Landsat observations served as a good alternative. However, there are a few months in winter when neither of the two missions had observations of the Aral Sea due to snow, ice and fog, there SSM based gap filling is useful.

Our estimated water height from Landsat and altimetry have a very high agreement with the ground observations (0.99 R²) for Lake Mead and with the ICESat-GLAS and Hydroweb database discussed in Duan et al. (2013) [17] and Gao et al. (2012) [19]. For the Aral Sea we compare the results from the two height estimation methods, because of the lack of in situ data for verification. The West Aral Sea has a better agreement between them, compared to the East Aral Sea (R²_west = 0.94 and R²_east2007 = 0.82) because of the inconsistencies in water mask due to classification errors of the shallow saline East Aral Sea; therefore, the linear relationship between Landsat and altimetry observations was disturbed at some time points. However, RMSE of the West Aral Sea is 11 cm more than the East Aral Sea (January 2002–December 2007), because a few pixels of classification error will have a greater impact on the height estimation of a deeper reservoir. Depending on the accessibility of the datasets, a combination of the given methods can be used, as shown in Figure 1. If all three of the aforementioned parameters, i.e., DEM, water level and areal extent, are available, their combinations can produce two independent absolute water volume time series (LBV and ABV) and one relative water volume time series (ALVV). Extrapolated digital surface models from different sources like TanDEM-X, ASTER DEM, and SRTM etc. can be used as an approximate bathymetry. However, they do not allow for the computation of absolute volumes. Nevertheless, the methods (LBV and ABV) can be applied on the extrapolated global DEM to obtain volumetric variations until one reaches close to the lake bed.

The estimated absolute volumes (ABV and LBV) and volumetric variations (ALVV) agree very well with each other (discussed in Section 4.4.2). Therefore, a combination of different approaches and satellites can be helpful to obtain a continuous time series of observations. Further, we assimilated the two absolute water volume time series by means of a Kalman-filter-based SSM to produce one water storage estimate (CSSME) and a six-month recursive forecast for two years for each reservoir, which also produced very good results (discussed in Section 4.4.3). The statistical parameters highlight the quality of regression of all of the volume estimation methods and the CSSME in Table 2.

The observed inconsistencies can be explained by the possibility of error in data acquisition [38,41] and its processing. In addition, for the large water bodies there is the possibility of a geoid anomaly. The presence of the Rocky Mountains could have result in a positive gravity anomaly in the Lake Mead region, and might have caused some constant offset (Figure 9) because of an uneven geoid. The observations from altimetry and Landsat-SRB are from its northern arm while in situ observations were collected near the Hoover dam located at the western end of the lake. The computational error for altimetry is mostly less than 5 cm with a maximum of up to 15 cm [17]. The possible range of uncertainty in Landsat observations from our algorithm for the Aral Sea is less than 1.5 m and for Lake Mead is mostly less than 4 m. Compared to altimetry, Landsat water heights have a larger error bar because the latter is highly dependent on the accuracy and spatial resolution (especially in steep terrain) of the mask and bathymetry. The spatial resolution of the mask and the bathymetry is often not the same; we therefore lose some existing information. With regard to the Aral Sea, the bathymetry had a very poor spatial resolution (though interpolated to 30 m) and may not be accurate. For example, changes in the structure of the lakebed from 1960s to the present East Aral Sea can be seen in Figure 12. Furthermore, Lake Mead and the West Aral Sea have very close contour lines, and there is a possibility of significant changes in the water height within a 30 m pixel from the Landsat mask. This is the reason why the error bar in water height extraction from the Landsat in Figure 7 is relatively high for these two water bodies. The other possible source of error is the uncertainty in the quality of the water masks derived from Landsat, which can have a few pixel registration errors or classification errors. The offset of a few pixels could induce severe problems in the height calculation, specifically for steep-sloped reservoirs.

Given the above stated many possible uncertainties, the remote sensing based estimates for reservoir height and volume from all reservoirs are promising and similar to ground observations. The results show that Lake Mead, the East Aral Sea and the West Aral Sea lost nearly 11.5, 27, 21.7 km³, respectively, of water volume over 12 years, while in the same period the North Aral Sea gained 6–7 km³ because of the construction of the Dike Kokaral dam in 2005. Lake Mead’s downward trend was interrupted once by a flood event between May 2011 and January 2012. This was caused by heavy snowmelt in the upper Colorado and Rocky mountains, prompting the release of an extra 4.1 km³ of water from Glen Canyon into Lake Mead, and an increase of nearly 12 m in water elevation. Continuous drought for more than a decade and increasing water and electricity demand from Las Vegas (one of the fastest-growing cities in the United States) have put Lake Mead under great strain.

The two test sites suffered significant water losses, but the effect is more pronounced on the water height observed from satellite altimetry from the hilly terrain of the Lake Mead (~28 m loss). However, the Aral Sea showed a dramatic loss in volume (~43 km³) by shrinking a few meters in height, because of its relatively flat topography. Therefore, a three-dimensional assessment of the changes is essential and can be achieved by combining observations from different remote sensing data from different perspectives to estimate surface water hydrology at on a global scale.

6. Conclusions

The main objective of this paper was to estimate reservoir water volume from remote-sensing-based methods. We proposed a new method of absolute volume computation from Landsat and bathymetry intersection (LBV), and compared the results with absolute volume time series derived from the interaction of satellite altimetry and bathymetry (ABV) as well as with the volumetric variations from an inverted pyramid model (ALVV). Very good long-term agreement was obtained between water heights observed by satellite altimetry and the Landsat-SRB method (Figure 7) as well as between all volumetric variation estimation methods (ABV, LBV and ALVV: Figure 8) for all reservoirs. Furthermore, we demonstrated the application of a state space framework as a plausible tool to obtain one time series estimate from multiple remote sensing data products because of its attractive features of data assimilation, missing data modifications and near-future forecast. Our results are in very good agreement with the ground observations and with one another. We tested the proposed methods in both gauged and ungauged basins; the results from all are very promising.

The results show that lakes and reservoirs can be effectively monitored at a global scale using remote sensing data and even an inaccurate bathymetry. With the combined use of the proposed methods, we have a good alternative to ground-based observations, particularly for small or remotely located, ungauged or poorly gauged reservoirs. Depending on the availability of data, combinations of the given methods can be used in global hydrological models to make updated dynamic estimates of surface water storage.