Operational reservoir inflow forecasting with radar altimetry: The Zambezi case study

. River basin management can greatly beneﬁt from short-term river discharge predictions. In order to improve model produced discharge forecasts, data assimilation allows for the integration of current observations of the hydrological system to produce improved forecasts and reduce prediction uncertainty. Data assimilation is widely used in operational applications to update hydrological models with in situ discharge or level measurements. In areas where timely access to in situ data is not possible, remote sensing data products can be used in assimilation schemes. While river discharge itself cannot be measured from space, radar altimetry can track surface water level variations at crossing locations between the satellite ground track and the river system called virtual stations (VS). Use of radar altimetry versus traditional monitoring in operational settings is complicated by the low temporal resolution of the data (be-tween 10 and 35 days revisit time at a VS depending on the satellite) as well as the fact that the location of the measurements is not necessarily at the point of interest. However, combining radar altimetry from multiple VS with hydrological models can help overcome these limitations. In study, a rainfall runoff model of the Zambezi River basin is built using remote sensing data sets and used to drive a routing scheme coupled to a simple ﬂoodplain model. The extended Kalman ﬁlter is used to update the states in the routing model with data from 9 Envisat VS. Model ﬁt was improved through assimilation with the Nash–Sutcliffe model efﬁciencies increasing from 0.19 to 0.62 and from 0.82 to 0.88 at the outlets of two distinct watersheds, the initial NSE (Nash–Sutcliffe efﬁciency) being low at one outlet due to large errors in the precipitation data set. However, model reliability was poor in one watershed with only 58 and 44 % of observations falling in the 90 % conﬁdence bounds, for the open loop and assimilation runs respectively, pointing to problems with the simple approach used to represent model error.


Introduction
Accurate short-term predictions of river flows are necessary for optimal river basin management, in particular for river systems with large reservoirs or in areas subject to flooding. The hydrological models used to generate river flow predictions are however subject to high uncertainties due to uncertain model structure, inputs, parameterization and initial conditions (e.g., Liu and Gupta, 2007).
In order to reduce prediction uncertainty, data assimilation can be used to combine the information from models and independent observations. Taking into account their respective uncertainties, models and data are combined to obtain the best possible estimate of the current state of the hydrological system. The improvements obtained from assimilation of in situ data to hydrological models, in particular water levels and discharge, have been successfully proven since the 1980s and data assimilation is commonly used in operational flood forecasting models (e.g., Kitanidis and Bras, 1980;Refsgaard, 1997;Madsen and Skotner, 2005). However, such applications require the availability of timely in situ data, which can be challenging in large remote river basins or in situations where riparian countries are unwilling to share their data. A solution to bypass such challenges is the use of remote sensing data. The direct measurement of river discharge from space is not possible with current technology, but radar altimetry can be used to track water level variations in surface water bodies. While initially designed for ocean monitoring, radar altimetry has successfully been used to measure river level variations in many areas of the globe (e.g., Koblinsky et al., 1993;Birkett, 1998;Berry et al., 2005;Frappart et al., 2006).
The two main challenges in using radar altimetry for hydrological models are the conversion of river level variations to discharge as well as the low temporal resolution, which has been of between 10 and 35 days for radar altimetry missions up to now.
Previous studies have focused on using radar altimetry in combination with hydrological models in order to overcome these limitations. Leon et al. (2006) and Getirana et al. (2009) obtained rating curves from discharge estimates of calibrated hydrological models. Getirana (2010) and  showed that altimetry could be used in the calibration of hydrological models with similar results to those obtained using in situ flow data.
Work preparing for the Surface Water Ocean Topography (SWOT) mission has shown that virtual wide-swath altimetry could be used to update hydrodynamic models and improve modeled depths and discharge (Andreadis et al., 2007;Biancamaria et al., 2011). The hydrodynamic models however rely on the availability of detailed bathymetric data that are not globally available. In terms of assimilation using nadir altimetry, Pereira-Cardenal et al. (2011) used altimetric measurements of reservoir levels to improve modeled reservoir levels and Paiva et al. (2013) showed that streamflow and water level forecasts in the Amazon could be improved through the assimilation of river altimetry.
The objective of this study is the assimilation of river altimetry into a routing model of the Zambezi River basin in order to improve inflow predictions for the Kariba and Itezhi-Tezhi reservoirs. The assimilation is carried out using the extended Kalman filter and model states are updated using altimetry data from the Envisat mission.

Study area: the Zambezi River basin
The Zambezi River basin is the largest of southern Africa, and the fourth largest in Africa. The basin covers over 1 390 000 km 2 and eight countries have land areas within its boundaries. Precipitation in the basin is highly seasonal with almost all of the rainfall occurring in the rainy season between the months of October and March.
The study focuses on two distinct watersheds, both located in the western part of the Zambezi River basin: (1) the Zambezi River upstream of Lake Kariba (draining approximately 5.2 × 10 5 km 2 ), and (2) the Kafue River upstream of Lake Itezhi-Tezhi (draining approximately 9.5 × 10 4 km 2 ) ( Fig. 1). Both Lake Kariba and Lake Itezhi-Tezhi are used for hydropower generation and their operation could benefit from improved predictions of inflow.
The two watersheds were divided into subbasins based on the availability of in situ data as well as altimetric virtual stations (VS), which are the locations where the satellite track and river network intersect. The outlet of watershed 1 is located approximately 160 km upstream of Lake Kariba and the outlet of watershed 2 is located approximately 90 km upstream of Lake Itezhi-Tezhi in order to coincide with in situ gauging stations.
The major feature located in the study area is the Barotse floodplain, which is located in watershed 1 (see Fig. 1) and has a storage capacity of 8.5 km 3 and extends over 7700 km 2 . The floodplain has a damping effect on flow through storage and evaporation and was taken into account in the routing scheme.
While the water resources in the Zambezi River basin are not currently subject to major stress, there are large variations in the temporal and spatial distribution of water availability across the basin. Water demand is also expected to increase rapidly with population and economic growth and there is a large potential for further development including plans for further hydropower installations and expansion of irrigated areas.
In this context, recent studies have focused on water management issues in the basin. Among these, Tilmant et al. (2010) focused on the optimization of reservoir operation in order to take into account tradeoffs between ecological conditions downstream and hydropower generation. Beck and Bernauer (2011) studied the combined effects of increased water demand and climate change. They predicted that water shortages were likely to occur in the basin, stressing the need to further develop water management strategies in the basin.
In order to achieve the objectives of efficient water management, reservoir inflow predictions are highly valuable. Assimilation of remotely sensed data to hydrological models can be used in order to improve discharge forecasts. One example of such a study is the work by Meier et al. (2011), who used remotely sensed soil moisture in a data assimilation framework to improve discharge forecasts with the objective of improving reservoir management. The aim of the current study is similarly to use data assimilation to improve inflow forecasts using radar altimetry.

Altimetry data
The altimetry data used in this study was the River AlTimetry (RAT) product developed at the Earth and Planetary Remote Sensing Lab (E.A.P.R.S.) . The RAT data product was obtained by retracking the 18 Hz Envisat waveforms. The data points have a 369 m along-track spacing and the return period for one virtual station is of 35 days. The data extraction procedure for the Zambezi River basin is detailed in Michailovsky et al. (2012).
The location of the 9 VS -6 in watershed 1 and 3 in watershed 2 -used in the study is shown in Fig. 1. The VS used are classified as "good" with standard errors (std) of less than 40 cm or "moderate" with expected standard errors of less than 70 cm. The values reported in Table 1 include amplitude adjustments for VS which had to be evaluated against in situ gauges at a different location along the same reach to account for cross section variability (see Michailovsky et al., 2012, for details on the uncertainty estimation and classification procedure).

Modeling
The altimetry data was assimilated to a routing model of the Upper Zambezi and Kafue rivers. The Barotse floodplain was modeled as interacting with the adjacent reaches through a first order exchange driven by head differences between the reach and floodplain. Inflows to the routing model were generated using a rainfall-runoff (RR) model of the study area.

Rainfall-runoff model
The rainfall-runoff model was built using the Soil and Water Assessment Tool (SWAT), which is a widely used semidistributed, semiphysically based hydrological modeling tool that operates on a daily time step (Neitsch et al., 2005). The SWAT model was chosen because it is well suited for largescale applications and is easily applicable in data sparse areas (Gassman et al., 2007).
The model was set up using remote sensing data only. The general set up as well as soil and vegetation parameters were taken from Schuol et al. (2008). Other data sets were taken from freely available data sources: the digital elevation model used was from the Shuttle Radar Topography Mission (SRTM, Farr et al., 2007), precipitation forcing was the Famine Early Warning Systems Network (FEWS-Net) rainfall estimate product (RFE) and temperature forcing was the European Centre for Medium Weather Forecast (ECMWF) ERA-Interim product.
The SWAT model was run using the Hargreaves method for the calculation of evapotranspiration, which requires only temperature data as input. The calibration of the RR model focused mostly on the groundwater parameters that were found to be the most sensitive parameters not related to soil and land cover data. Such parameters were not calibrated in order to preserve the physical representation and reduce the number of calibration parameters. Table 2 presents the main calibration parameters and their values.

Reach routing
Channel routing was modeled using a Muskingum routing scheme expressed in terms of water storage. The propagation model for stored water volume in the Nth reach downstream can be written as where A N is defined as and C 1N and C 3N are defined as in Chow et al. (1988): where s N,j is the storage in reach N at time step j (m 3 ) and M N,j is the rainfall-runoff model generated inflow to reach N at time step j (m 3 s −1 ) and t is the model time step (days). The two Muskingum parameters, X N a weighing factor and K N , the travel time of the flood wave through the reach (days), were assumed constant for each reach segment (e.g., Chow et al., 1988, p. 258).

Floodplain model
Because storages in Eq. (1) are expressed only as a function of the states of the previous time step, the routing through the reaches and the floodplain processes can be carried out sequentially. A simple floodplain model was built following an approach similar to that used by Dincer et al. (1987) to model the Okavango swamp. Two processes were modeled in the floodplain: water transfers with the main reach and evapotranspiration from the floodplain. Direct rainfall onto the floodplain was not considered here as it is already taken into account in the rainfall-runoff model. The floodplain-reach interaction was modeled as a first order exchange driven by the difference in water levels between the floodplain and reach. The open water evaporation rate was assumed equal to the potential evaporation (PET) rate from the subbasin in which the floodplain is located. The PET value was obtained from the RR model. The basic equations for the floodplain-reach interaction for a floodplain located in a reach "rc" are where V fp is the floodplain volume (m 3 ), coeff is the transfer coefficient (m 2 s −1 ), h rc and h fp are the water levels in the reach and the floodplain, A fp is the floodplain area (m 2 ), ET 0 is the potential evaporation (ms −1 ), s rc is the water stored in the reach (m 3 ), msk is the Muskingum routing operator as presented in Eq.
(1), s is the state vector of volumes in all reaches and M the inputs from the RR model (m 3 s −1 ). The one-day time step used for the modeling was assumed small relative to the timescale of the floodplain processes and the floodplain equations were therefore solved assuming mean daily volume in the floodplain equal to the volume at the end of the previous day, minus the evaporation that was assumed to be removed before any transfers take place. The explicit solution is then where h rc,k * is the level in the reach after the addition/subtraction of volume from the Muskingum routing but before any transfers with the floodplain have taken place.
The geometry of the reach and floodplain need to be known in order to obtain levels and floodplain areas from the modeled water volumes. Figure 2 presents the geometry that was assumed for the reach and floodplain on one side of the reach. The floodplains extend on both sides of the reach, and are assumed symmetrical with respect to the reach.
Hydrol. Earth Syst. Sci., 18, 997-1007 The reaches were assumed to have trapezoidal cross sections with constant bank slope, α b , and the elevation of the bottom of the floodplain was assumed to rise with distance from the reach following where β and m are shape parameters and x is the distance from the side of the floodplain closest to the reach (Fig. 1). The shape parameters were determined by extracting SRTM heights along one representative cross section per reach where the floodplain is located and fitting the β and m parameters.
Reach width and bank slope were determined based on Landsat imagery. Bank slope was estimated by measuring low and high flow widths as well as high and low flow altimetric heights from the same location. Bank slope was then calculated as In order to obtain bottom widths, a low flow depth was computed based on historical measured low flows (or modeled low flows where unavailable) using a method detailed in Michailovsky et al. (2012). The calibration of the RR model was carried out manually using in situ discharge data. The Muskingum K parameter as well as the floodplain exchange coefficient, coeff, were calibrated using both in situ flows and altimetric levels.

The extended Kalman filter
The Kalman filter is a sequential data assimilation scheme that can be split in a propagation and an analysis phase. When the measurement and model operators are nonlinear, the Extended Kalman filter (EKf), which is the nonlinear extension to the Kalman filter, can be used.
In the propagation phase of the standard Kalman filter, a forecasted state and covariance are calculated using where s f and P f are the forecasted state vector and state covariance matrix, u is the model forcing, w is a sequence of white Gaussian noise with covariance Q, F is the state transition matrix, G the control input matrix and the noise input matrix.
In the analysis phase, the analysis or updated state vector and covariance at a time m when a measurement is acquired are obtained through the following equations: where s a and P a are the analysis state vector and covariance matrix, and H is the measurement operator, which is defined as where y m is the measurement at time m and v is a sequence of white Gaussian noise with covariance R m . The quantity (y m − H m · s f m ) is called the innovation or measurement residual.
In the EKf, the nonlinear model and measurement operators are used directly in Eqs. (10) and (14). The H and F matrices needed in Eqs. (11)-(13) are obtained by linearizing the measurement and model operators around the forecasted state. So if h denotes the nonlinear model operator and f the nonlinear model operator H = ∂h ∂s s=s f and F = ∂f ∂s s=s f .

Measurement operator
The measurement operator, h, maps the model state in the observed space; i.e., in this study, h is used to convert the stored volume in a reach to an altimetry reading. Reaches were assumed to have a trapezoidal cross section with bottom width w, bank slope α b and length L. The storage in the reach, s, can then be expressed as a function of the water depth d as Solving for depth yields Finally, a common reference was needed between modeled depth and measured altimetry. This was done by running the routing model and shifting the altimetric heights by the difference in mean between coincident modeled depths and measurements over the calibration period leading to the definition of h: where (alti) are the altimetric height measurements.

Error model
The measurement error on the altimetry values was assumed normally distributed with zero mean. Standard error estimates were based on the values reported in Michailovsky et al. (2012) (see Table 1). Specification of model uncertainty is one of the major tasks in data assimilation because of the many different sources of error, which are poorly known and typically extremely difficult to separate from one another (Liu and Gupta, 2007).
The approach chosen for this study was to assume that the rainfall-runoff forcing was the main source of model error. As the magnitude of the error on model-generated runoff is typically proportional to the magnitude of the runoff, the error was applied as a multiplicative term on the RR forcing and the model error representation was therefore determined by analyzing the normalized runoff residuals during the calibration period.
In order to obtain in situ measurements of runoff, gauged flow was assumed equal to runoff for upstream catchments. For catchments located further downstream, gauged runoff from a given area was assumed equal to the difference between downstream and upstream gauged runoff.
While model error is assumed white Gaussian (see Eq. 11), the runoff residuals were found to be highly autocorrelated. This was taken into account by assuming a first-order autoregressive (AR1) model: where (w) are the runoff residuals, α is the AR1 parameter and ε is a sequence of Gaussian white noise, which has a covariance Q . This type of AR1 error model can easily be implemented in the EKf by augmenting the state vector with the correlated noise term. By setting (e.g., Jazwinski, 1970) where I is the identity matrix and all other terms have been defined previously, Eq. (10) can then be written as and all other equations for the EKf can then be directly applied using the matrices and vectors defined in Eq. (19). Because of the low dimensionality of the problem, the added states will not cause any significant computational burden for the application: in this case study it will lead to having two states per reach plus one state per floodplain, which is a total of 32 state variables.

Model evaluation
In any model prediction, the performance can be evaluated using a number of different criteria in order to characterize its performance in terms of accuracy (i.e., how close the value of the model estimate is to the observations) and precision (i.e., how uncertain the model prediction is).
In order to fully assess the performance of the open loop and assimilation model runs, the following measures were used.
-Coverage: the percentage of observations that fall within the predicted nominal confidence interval.
-Nash-Sutcliffe efficiency (NSE) -Root mean square error (RMSE) -Sharpness: the width of the predicted nominal confidence interval.
Because a tradeoff must usually be made between sharpness and coverage, a measure combining both criteria, the interval skill score (ISS) was also used. The ISS is defined as follows (Gneiting and Raftery, 2007): where where u and l are the upper and lower confidence bounds at the significance level α for the model estimate and x is the observed value. The ISS should therefore be minimized as a lower ISS value will indicate sharper confidence intervals and fewer observations located outside of the confidence bounds.

Results and discussion
Model calibration was carried out over the years 2001-2004 (calibration of Muskingum's K using altimetry levels could only be carried out from October 2002, which is the date from which Envisat data is available). The model was run up to the end of 2008 as this was the most recent year for which we were able to obtain in situ data for model validation.
Model calibration yielded NSE values between 0.55 and 0.90, which decreased to values between 0.07 and 0.82 over the validation period (Table 3). The low values at some subbasin outlets can be explained by the large errors in the precipitation data set. In particular, in 2005 large differences were observed between the different remote sensing precipitation data sets available. By excluding the year 2005, the NSE values over the validation period increased significantly and were found to be between 0.46 and 0.87.
Analysis of the RR residuals yielded AR1 parameters between 0.9918 and 0.9978 and the standard deviations on the white Gaussian error were found to be between 0.035 and 0.49. For some subbasins, this produced unrealistically high estimates of model error. In particular for subbasin 17 as well as for subbasins located downstream of the Barotse floodplain.
For subbasin 17, further investigation showed that applying the method detailed in Sect. 4.3.3 led to using only residuals from the year 2004 because it is the only year for which data was available upstream. However 2004 is a poor year in terms of model performance for subbasin 17 (Fig. 3) and the runoff residuals were therefore analyzed over the whole of watershed 2 to obtain the error parameters at subbasin 17.
For the Barotse floodplain, the assumption that the error is mainly attributable to the RR forcing breaks down as the floodplain model uncertainty is not taken into account though high uncertainties are expected on the volume-area relationship, the transfer coefficient and the ET rate from the floodplain. This led to the attribution of unrealistically high errors to the RR forcing in the residual analysis. For subbasins downstream of the floodplain (24, 29, 32 and 34), the RR error parameterization from the nearest upstream subbasin (sub. 14) was therefore used. Table 4 presents the results for all subbasins where in situ data from the validation period (post-2005) was available and Fig. 3 presents the results graphically at the outlets of the two watersheds.
Major improvements in RMSE and NSE were observed in all subbasins with the range of NSE values at the watershed outlets going from 0.19 to 0.62 for watershed 1 and from 0.82 to 0.88 for watershed 2.
All subbasins showed improvements in all measures except for coverage, which only improved for subbasin 12. The loss of coverage was particularly significant for subbasin number 34 where approximately 14 % fewer observations fell within the confidence bounds after assimilation. As a consequence, subbasin 34 showed an increase in the ISS through assimilation meaning that the loss of coverage outweighed the gains in sharpness at this outlet. The loss of coverage between open loop and assimilation runs can indicate an underestimation of the measurement error that is to be expected in this case as only measurement error itself was taken into account while the error term should also include measurement operator uncertainty. Inspection of the residuals showed that 13 % of residuals fell outside of the 90 % confidence bounds predicted by the filter equations for the calibration period, and 17 % for the validation period. The assimilation was run again adding a constant to the estimated standard deviation of the altimetry measurements. Adding 0.2 m led to a good match between the innovations and their predicted statistics with 10.4 and 11 % of innovations falling out of the 90 % confidence bounds for the calibration and validation time periods respectively.
The results using the new measurement error are presented in Table 5 and show an improvement in the coverage for the assimilation run compared to the previous run. However, the coverage was still slightly degraded at most outlets relative to the open loop run. For outlets 12 and 17, where the coverage was already good, almost no change was observed in the results suggesting that the initial measurement error was adequate. Further improvements could potentially be obtained by analyzing the virtual stations and measurement operators at each location in order to adapt the measurement error term.
The results in both cases show that the main weakness in the assimilation scheme is the representation of model errors. This can be observed in the coverage values, in particular for reach 34 where in the open loop run only 58 % of observations fall within the 90 % confidence bounds. This issue particularly affects subbasins downstream of the Barotse floodplain where the open loop model performance is poor, which suggests both that the modeling in the area needs to be improved, and that in complex river environments including floodplains, a separate error term representing the floodplain processes is needed.
A first step towards better error representation would be to take evaporation uncertainty into account by adding, for example, a multiplicative error term on the ET forcing. A separate additive error term could also be added to floodplain subbasins. However, as the model error representation becomes more complex, it may be necessary to use a different assimilation scheme such as the ensemble Kalman filter in which nonlinear processes can more easily be taken into account. Figure 3a presents the results for the outlets of watershed 1 and shows that improvements from the assimilation are not equally distributed over the simulation time period. For example, in 2007 the assimilation performs poorly for reach number 34. Figure 4b shows the timing of the altimetric measurements between October 2006 and October 2007 and it can be observed that only one altimetric measurement is available over a period of 67 days, between 9 February and 17 April 2007, and that the update that is carried out in this period degrades model performance. The satellite passes occur on different days at the different VS over the 35 day repeat period and the maximum delay between satellite passes at any VS within watershed 1 is of 16 days. However many Hydrol. Earth Syst. Sci., 18, 997-1007 of the VS are located on narrow rivers (< 200 m wide) and this increases the risk of no data points being acquired or of incorrect values/outliers because only one data point will be available per satellite pass.
In contrast, Fig. 4a shows that when no such data gaps exist the assimilation performs well; the use of multiple VS over the watershed compensating the low temporal resolution at each individual VS.
For watershed 2, the problem is magnified by the fact that only 3 VS are used for the update and that these 3 VS are all visited by the satellite within 6 days of the 35 day repeat period. Therefore, while the benefits from assimilation are clear in a year where the model consistently over-or underpredicts discharge such as in 2005 (Fig. 3), the altimetry data set will be unable to capture the shorter-term flow variability as is shown in Fig. 5.

Conclusions
In this study, radar altimetry from the Envisat satellite was used to update the reach storages in a Muskingum routing scheme coupled to a simple floodplain model and driven by the output of a rainfall-runoff model. Assimilation improved Nash-Sutcliffe model efficiencies from 0.19 to 0.62 and from 0.82 to 0.88 at the outlets of two distinct watersheds located upstream of Lake Kariba and Lake Itezhi-Tezhi. Model reliability was good for the outlet of watershed 2 but was found to be low at the outlet of watershed 1. This was due to lower model quality and error representation in watershed 1, which is more complex hydrologically, including the large Barotse floodplain. The study also highlighted the limitations of using the current altimetric data set in areas where only few VS are available due to its low temporal resolution. It was found that when the initial model run is of high quality, the assimilation only marginally improves the performance. However when the initial model does not perform well, in particular due to large errors in input data, it was shown that the assimilation or radar altimetry had the potential to significantly improve results. While this study has only used altimetry over rivers, floodplain levels can also be tracked through altimetry and further work including altimetric data over the floodplain could potentially improve the results obtained in and downstream of the Barotse floodplain.
Nonetheless, the high potential for the use of radar altimetry in assimilation has been demonstrated as the use of multiple VS was able to compensate for the low repeat period of the satellite where sufficient VS were available. These results should be greatly improved in the future, with higher spatial resolution altimeters or swath altimetry as planned in the Surface Water Ocean Topography (SWOT) mission, which will allow for more, narrower, rivers to be monitored through altimetry.