Water Resources Research

RESEARCH ARTICLE Discharge Estimation in Ungauged Basins Through Variational

10.1002/2017WR021735
Data Assimilation: The Potential of the SWOT Mission
Key Points: H. Oubanas1 , I. Gejadze1, P.-O. Malaterre1, M. Durand2,3 , R. Wei3, R. P. M. Frasson3 , and
 River discharge is estimated using
A. Domeneghetti4
synthetic remote sensing
measurements from the forthcoming 1
National Research Institute of Science and Technology for Environment and Agriculture, Montpellier, France, 2School of
SWOT mission without any in situ
observation Earth Sciences, The Ohio State University, Columbus, Ohio, USA, 3Byrd Polar and Climate Research Center, The Ohio State
 2D SWOT radar measurements are University, Columbus, Ohio, USA, 4School of Civil Engineering, Department DICAM, University of Bologna, Bologna, Italy
simulated and mapped onto river
centerlines, producing realistic errors,
temporal, and spatial sampling Abstract Space-borne instruments can measure river water surface elevation, slope, and width. Remote
 Discharge accuracy was 12.1% and
11.2% on the Po and Sacramento
sensing of river discharge in ungauged basins is far more challenging, however. This work investigates the
Rivers, respectively, illustrating estimation of river discharge from simulated observations of the forthcoming Surface Water and Ocean
potential for ungauged basins Topography (SWOT) satellite mission using a variant of the classical variational data assimilation method
‘‘4D-Var.’’ The variational assimilation scheme simultaneously estimates discharge, river bathymetry, and
Supporting Information: bed roughness in the context of a 1.5 D full Saint-Venant hydraulic model. Algorithms and procedures are
 Supporting Information S1
 Data Set S1
developed to apply the method to fully ungauged basins. The method was tested on the Po and Sacra-
mento Rivers. The SWOT hydrology simulator was used to produce synthetic SWOT observations at each
Correspondence to: overpass time by simulating the interaction of SWOT radar measurements with the river water surface and
H. Oubanas, nearby land surface topography at a scale of approximately 1 m, thus accounting for layover, thermal noise,
[Link]@[Link] and other effects. SWOT data products were synthesized by vectorizing the simulated radar returns, leading
to height and width estimates at 200 m increments along the river centerlines. The ingestion of simulated
Citation: SWOT data generally led to local improvements on prior bathymetry and roughness estimates which
Oubanas, H., Gejadze, I.,
Malaterre, P.-O., Durand, M., Wei, R., allowed the prediction of river discharge at the overpass times with relative root mean squared errors of
Frasson, R. P. M., et al. (2018). 12.1% and 11.2% for the Po and Sacramento Rivers, respectively. Nevertheless, equifinality issues that arise
Discharge estimation in ungauged from the simultaneous estimation of bed elevation and roughness may prevent their use for different appli-
basins through variational data
assimilation: The potential of the cations, other than discharge estimation through the presented framework.
SWOT mission. Water Resources
Research, 54, 2405–2423. [Link]
org/10.1002/2017WR021735

Received 26 AUG 2017

1. Introduction
Accepted 5 MAR 2018
Accepted article online 8 MAR 2018 As a primary source of fresh water, rivers are among the most important natural resources and are the nerve
Published online 30 MAR 2018 system of ecology and human society. Since the earliest civilizations, continental waters have been central
to society, featuring in applications such as water supply, irrigation, drainage, navigation, fisheries, flood
control, hydropower generation, wastewater treatment, pollution abatement, and wildlife protection (Benke
& Cushing, 2011; Chow, 1964; Tockner et al., 2009).
River discharge represents the flow of water from continental to oceanic environments, and is one of the pri-
mary quantities of interest in characterizing fluvial environments. Despite the undeniable importance of global
rivers, the availability of the world’s in situ gauge stations have been steadily declining since the late 1970s due
to political, economical, and geographical reasons (Fekete & Vo €ro
€smarty, 2002; Sneeuw et al., 2014; Tourian
et al., 2013). Moreover, the reluctance from countries to share data in a timely fashion makes the monitoring of
discharge and the flood forecasting in international rivers a daunting task (Biancamaria et al., 2011a; Hossain
et al., 2014). Therefore, the need for alternative and/or complementary inland water measuring techniques,
such as space-borne sensors, has become a primary concern of the scientific community and space agencies.
In this respect, several nadir altimeters that observe the water height are becoming widely available for the
land surface applications, such as JASON-1, JASON-2, ENVISAT, CryoSat-2, SARAL/ALTIKA and the recent
near-real time operational missions, JASON-3 and Sentinel-3 (Alsdorf & Lettenmaier, 2003; Cretaux & Birkett,
2006; Calmant & Seyler, 2006; Cretaux et al., 2011; Papa et al., 2006). Most of them are focused on observing
C 2018. American Geophysical Union.
V water surface elevation (WSE) in large rivers and lakes (Bartsch et al., 2009; Birkett, 1995, 1998; de Oliveira
All Rights Reserved. Campos et al., 2001; Frappart et al., 2006; Kouraev et al., 2004; Medina et al., 2008; Papa et al., 2003).

OUBANAS ET AL. 2405

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

The forthcoming SWOT mission, scheduled for launch in 2021, will provide a global mapping of the conti-
nental water bodies using the wide swath radar interferometer KaRIN (Ka-band Radar INterferometer).
SWOT will observe two ground swaths of 50 km separated by a 20 km nadir gap and will provide observa-
tions of water surface elevation, slope, and river width of rivers wider than 100 m, with a possibility of pro-
viding accurate observations of rivers as narrow as 50 m (Rodriguez, 2016). The SWOT repeat cycle is 21
days, during which sites will be revisited typically 2–4 times at irregular intervals. The temporal frequency is
dependent on the latitude of the region of interest, with low observation frequency near the equator and
higher observation frequency in higher latitude (Biancamaria et al., 2010, 2016).
Several algorithms designed to use future SWOT observations for river discharge estimation have been pro-
posed. Among the proposed methodologies, data assimilation (DA) has become increasingly popular within
the hydrological and hydraulic communities. DA methods allow the combination of available observations,
prior knowledge and expertise, and natural system dynamics, with a treatment of the associated errors, to
provide the best estimates of the unknown system variables and/or parameters. Early DA studies assumed
the river bathymetry and the bed roughness to be known a priori (Andreadis et al., 2007; Biancamaria et al.,
2011b). Durand et al. (2008) and Yoon et al. (2012) estimated the river bathymetry, but assumed that the
bed roughness was well known. The complexity of implementing DA methods with dynamical models and
the difficulty of estimating discharge together with the river bathymetry and the bed roughness have led to
the use of simplified models derived from the Saint-Venant equations. Four such algorithms have been pre-
sented in the literature, namely Garambois and Monnier (GaMo) (Garambois & Monnier, 2015), Metropolis
Manning (MetroMan) (Durand et al., 2014), the Mean-Annual Flow and Geomorphology (MFG) (Durand
et al., 2016), and Mean Flow with Constant Roughness (MFCR) (Durand et al., 2016). These algorithms were
assessed in Durand et al. (2016) for 19 rivers with various hydraulic conditions. At least one algorithm
(although not always the same one) had a discharge estimate relative root mean squared error less than
35%, on the 14 nonbraided rivers included in the study.
Among DA methods, variational DA has been the preferable approach in operational geophysical applications,
such as numerical weather prediction (NWP) (Courtier et al., 1998; Fischer et al., 2005; Gauthier et al., 1999,
2007; Rabier et al., 2000). Within the meteorology and oceanography communities, these methods are com-
monly named in their classical form ‘‘3D-Var’’ and its temporal extension ‘‘4D-Var.’’ The problem is formulated as
an optimal control problem and can be considered as a special case of the maximum a-posteriori probability
estimator (MAP) in the Bayesian framework. The optimal estimates of the unknown model variables and/or
parameters such as initial and/or boundary conditions, source terms (forcing), distributed and/or lumped coeffi-
cients (e.g., roughness coefficient, bathymetry, etc.), are obtained via the minimization of a well-defined cost
function. Gradient-based optimization methods resolve the minimization problem (e.g., quasi-Newton, conju-
gate gradient) and require the derivation of the tangent linear and adjoint models associated with the consid-
ered direct model. References on variational DA applied to the 1D full Saint-Venant model are rather scarce
(Ding & Wang, 2012). The adjoint model in the latter reference is derived analytically, then implemented
numerically, i.e., ‘‘optimize-then-discretize’’ approach, which yields an ‘‘inconsistent’’ adjoint. The disadvantage
of the latter approach is lower accuracy of the gradient and difficulty in applying in the framework of con-
strained optimization. The consistent adjoint for the full Saint-Venant hydraulic model Simulation and Integra-
tion of Control for Canals (SIC2) has been reported only recently in Gejadze and Malaterre (2017, 2016) and
Oubanas et al. (2018). This model has been developed at the National Research Institute of Science and Tech-
nology for the Environment and Agriculture (IRSTEA). The difficulty behind the computation of the stable
adjoint model is one reason why variational DA has not been very popular within the hydraulic research com-
munity. In this case, alternative filtering methods, i.e., the Kalman filter and its ensemble extensions, are widely
used instead. Note that for two-dimensional (2D) Shallow Water equations-based models, variational DA has
been reported in Lai and Monnier (2009) and Hostache et al. (2010) for different problem setups, e.g., single
variable estimation, data types, temporal, and spatial scales.
In the present paper, a variant of the classical ‘‘4D-Var’’ method, presented in Gejadze and Malaterre (2017,
2016), has been adapted to the general and realistic case of ungauged rivers observed from space. This
method, based on the iterative regularization technique, is more suitable for the nonlinear systems where het-
erogeneous variables, e.g., the state (discharge), the parameters (roughness), and the domain (bathymetry), are
estimated simultaneously. Such configuration, involving nonlinear operators, may compromise the robustness
of the minimization problem (due to large norm differences among gradient components), which requires

OUBANAS ET AL. 2406

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

appropriate weight assignment to different contributions. Using synthetic SWOT observations produced by the
SWOT simulator, we demonstrate the simultaneous estimation of discharge, river bathymetry, and roughness
for 1 year over a 133 km section of the Po River and for a 6 month period over a 153 km section of the Sacra-
mento River. Note that the prior knowledge of the variables of interest, needed for the optimization method, is
derived from the SWOT observations and globally available ancillary information only which enables the appli-
cation of the proposed method for ungauged as well as gauged basins.
The outline of paper is as follows. In section 2, we present the SWOT simulator and the study areas and
period. The methodology is described in section 3 introducing the 1.5 D full Saint-Venant hydraulic model,
the variational data assimilation method with its sequential version, and handling of a priori information.
Experimental design is presented in section 4. Section 5 presents and discusses the results of numerical
experiments. Finally, main findings of this study are summarized in section 6.

2. Study Area, Period, and Simulated Data

2.1. The Po River
The Po River is the longest river that flows entirely in the Italian peninsula, running from the northeastern
Alps to the Adriatic Sea. The mean channel is about 650 km long with over 141 tributaries, draining an area
of approximately 71,000 km2, corresponding to the largest catchment in Italy. The Po River is also the larg-
est Italian river in terms of the streamflow with a maximum historical discharge of 13,000 m3 s21 observed
at Pontelagoscuro in 1951. The Po Valley has experienced intensive agricultural and industrial development
during the twentieth century, which have led to an increasing vulnerability to hydrological hazards (Monta-
nari et al., 2017).
The modeled stretch of the river is 133 km in length, and is located between the gage of Borgoforte at the
upstream boundary and the city of Corbola at the beginning of the Po River delta (Figure 1). It is observed
by three passes during each of the 21 day SWOT cycle, for which the corresponding SWOT simulations are
available, i.e., the left and right swaths of the pass 0560, the left swath of the pass 0211, and the right swath
of the pass 0489 (see Figure 1). Note that ascending passes are labeled with pair identification number
while odd numbers refer to descending passes. The study period is 1 year long, ranging from May 2008 to
April 2009 during which a hydrodynamic model that solves the Saint-Venant equations were applied to the
Po River. It uses the flow hydrograph as an upstream boundary condition to simulate the flow at a daily
scale and is conditioned downstream by the observed water surface elevation. The flow hydrograph
reported in Figure 2 shows a discharge variability from 565 to 6,850 m3 s21. In this figure, the SWOT

Figure 1. The Po River study area (V

C Google Earth) with the ground track of SWOT overpasses (50 km swaths at each side of the 20 km nadir gap) 0560 (grey),

0211 (blue), and 0489 (yellow).

OUBANAS ET AL. 2407

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

Figure 2. Daily discharge data at the Borgoforte streamgage (black) at the Po River during the study period with vertical
dashed lines indicating the timing of the SWOT overpasses.

temporal sampling is illustrated by the vertical lines, which represent the time of the SWOT overpasses,
with the line color identifying the overpass identity as defined in Figure 1. Note that these lines serve as an
indication of the time of the SWOT observations that are assimilated to estimate the upstream discharge at
Borgoforte. They are spatially distributed dependent on the corresponding swath coverage.

2.2. The Sacramento River

The Sacramento River rises in the Klamath Mountains and drains an area of 69,000 km2 in northern Califor-
nia including the Coast Range, the Sierra Nevada, the Modoc Plateau, and the Trinity Mountains. This river is
known to be the largest in the state flowing for 640 km long and the major freshwater source for the San
Francisco Estuary. Since early 1940s, the flow of the Sacramento River has been regulated by dams to con-
trol the magnitude of the flow discharge and the frequency of the flooding events (Buer et al., 1989).
The region of interest is 153 km long and is located between Hamilton and Tyndall Landing cities, upstream
of the city of Sacramento. The study area is observed entirely by the right swath of the pass 0014, while the
far range portion of the left swath of the pass 0292 intersects a small downstream section of the study area
(see Figure 3). The study period is a 6 month simulation from January 2009 to June 2009 conducted with a
one-dimensional hydraulic model of the basin released in 2013 by the California Department of Water
Resources as part of the Central Valley floodplain evaluation and delineation program (Rogers, 2014). The
model has been simplified to remove diversions, tributaries, and storage cells. The hydrograph during this
period is presented in Figure 4, covering a range of discharge from 115.31 to 1541.97 m3 s21. Similarly, the
SWOT temporal sampling is illustrated by the vertical lines, with the line color identifying the overpass iden-
tity as defined in Figure 3.

2.3. Deriving 2D Hydraulic Properties From 1D Hydraulic Simulations

Accurate SWOT simulation requires high-resolution elevation maps, whereas the hydraulic simulations for
both the Po and Sacramento were performed in one dimension using different models. The Po River was
simulated by a quasi-two-dimensional hydraulic model built from a combination of a 1 or 2 m resolution
LiDAR Digital Elevation Model (DEM) of the area and boat surveys of the river (Castellarin et al., 2011; Dome-
neghetti et al., 2014) leading to 111 cross-sections irregularly spaced, with an average distance between the
cross sections of 1.2 km. The quasi-2D scheme depends on the presence of floodplain areas, which are pro-
tected from frequent inundations by a system of minor embankments and are connected to the main chan-
nel by means of lateral structures. The model considers these floodplains as storage areas in which the

OUBANAS ET AL. 2408

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

Figure 3. The Sacramento River study area (V

C Google Earth) with the ground track of SWOT overpasses (50 km swaths at each side of the 20 km nadir gap) 0014

(blue), 0292 (yellow).

water level is controlled by a level-volume curve estimated from the DEM and depends on the flow
exchange between the main channel and the storage areas (see Castellarin et al. (2011) for more details
about the model). Similarly, the Sacramento River hydraulic model was built using a combination of air-
borne LiDAR with land and boat surveys of the area leading to 601 cross sections, irregularly spaced with an
average distance of 258 m. In both cases, the estimation of the 2D water surface elevation required by
SWOT simulator has been performed interpolating the water elevation associated with all cross sections to
a regular grid. In particular, the interpolation refers to water elevations mapped laterally perpendicular to
the river centerline and has been performed with tools suitable for such spatial interpolation (such as, e.g.,
HEC-GeoRAS). These operations resulted in temporally dynamic high-resolution maps of water inundation
and elevation, combined with topographic features.

Figure 4. Daily discharge data at Hamilton streamgage (black) during the study period with vertical dashed lines indicat-
ing the timing of the SWOT overpasses.

OUBANAS ET AL. 2409

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

2.4. SWOT Simulator

The SWOT simulator designed by the Jet Propulsion Laboratory (JPL) is used for the first time in this work
for DA purposes to simulate the expected performance of the KaRIn instrument on-board the SWOT satel-
lite. It was used for the first time by our group in this paper for DA purposes. The SWOT simulator takes as
inputs the SWOT orbit and radar parameters, the expected water-land radar contrast, the DEM of the study
area including water and terrain, and a water mask, i.e., a two-dimensional map that distinguishes between
inundated areas and dry land. The water DEM, i.e., a digital elevation model that contains only the elevation
of the water surface, and the water mask are built based on the results of the dynamic hydraulic simulations
of the rivers at the time of the satellite overpass. As a first step, the SWOT simulator builds synthetic radar
interferograms of the scene, containing no noise and only subjected to terrain layover errors, i.e., errors that
happen when the radar return from two or more distinct targets, generally water and surrounding terrain at
higher elevations, reach the satellite at the same time, leading to overestimation of the water surface eleva-
tion (for more information of terrain layover, refer to Fjortoft et al., 2014). Subsequently, the simulator adds
noise to the synthetic interferograms, which are modeled as correlated circular Gaussian noise. The noise
added to the target depends on the surface type, i.e., land or water, the position in the swath, and radar
parameters. As the next step, the simulator processes the noisy interferograms to create pixel cloud, which
is composed of a pixel cloud with target class, elevation, and area associated with the pixel.
The processing of the noisy interferograms entails the following steps: a smoothing procedure called multi-
looking, target classification, and geolocation. The multilooking step is a procedure that averages the
returned power from consecutive pixels in the along-track and cross-track dimensions, effectively decreas-
ing noise at the cost of reduced spatial resolution (Cuchi, 1986; Ulaby et al., 2014). The number of looks
used in the present study was 4, which is currently envisioned as the lowest level of smoothing for SWOT
data products. Next, the simulator classifies the target based on the returned power. The classification relies
on the expectancy that for SWOT band and antenna characteristics, water will appear brighter than land tar-
gets (Alsdorf & Lettenmaier, 2003; Fjortoft et al., 2014). Ulaby and Dobson (1989) provide a list of expected
backscatter coefficient for different targets when viewed by a range of viewing angles. Given SWOT’s view-
ing angle, the contrast between water and land is conservatively on the order of 10 dB (Biancamaria et al.,
2016), which would allow differentiation based on thresholding. However, sharp corners in buildings may
appear as bright as water, as seen the images published by Fjortoft et al. (2014). Filtering of false detections
can be done by searching for nearby bodies of water or by implementing more advanced classification
techniques. Finally, the simulator geolocates the targets, translating the pixels from radar coordinates into
elevation, latitude, and longitude (for more information of terrain layover, refer to Fjortoft et al., 2014).

2.5. Deriving 1D Simulated Observations From 2D SWOT Simulations

The pixel clouds were processed with the RiverObs package developed at JPL. RiverObs aggregates the 2D
pixel clouds produced by the SWOT simulator into regularly spaced points, called nodes, located at the river
centerline and estimates node-averaged height, width, and associated observational uncertainties. RiverObs
produces node statistics by assigning pixels located within a user-defined search window to the nearest
river node. The search window is a polygon with outer boundaries running parallel to the river centerline. In
the present work, we utilized a search window of 1,200 m for the Sacramento River which corresponds 6–
12 times the averaged width of the river. For the Po River, we used a 1,600 m wide search window for all
but the three highest flow overpasses, for which we increased the width to 5,000 m to account for overbank
flow. This represents 4 times the averaged width during the low flow and 9 times during high flow (Frasson
et al., 2017). The node height assumes the value of the average of all water pixel heights associated with
that node whereas the node width is estimated by dividing the inundated area associated with the node by
the node spacing. We use the synthetic SWOT node-averaged products within our data assimilation frame-
work to estimate river discharge as described in sections 3.2 and 3.3. The Po River is observed 52 times dur-
ing the 1 year study period while the Sacramento study area is observed by 18 overpasses during the 6
month study period. The simulated node-averaged water surface elevation and the top river width, along
the Po and Sacramento Rivers, are presented in Figures 5 and 6. In the case of the Sacramento River, the
width variations are smaller after 90 km due to the presence of levees, which confine the river during our
simulations. In the first 90 km, the river is wider which allows to decrease the random error in the WSE
measurements. The errors are more important in the downstream part because the river is narrower. More-
over, the presence of the levees increases the layover errors.

OUBANAS ET AL. 2410

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

Figure 5. Plots A, B, and C show the river width measured at the nodes for the Po River using the simulated SWOT passes 0560 (18 superimposed revisits), 0211 (17
revisits), and 0489 (17 revisits), respectively. Plots D, E, and F show the water surface elevation at nodes produced from the passes 0560, 0211, and 0489, respectively.

The simulated observations are subjected to biases, although not all the sources of systematic errors are
taken into account. The bias depends on the orientation of the river with respect to the satellite swath; it is
more important when the river is parallel to the swath (e.g., case of the Sacramento River) and is negligible
when the river is perpendicular to the swath (e.g., case of the Po River). For more information about the
characteristics of the SWOT data, one can refer to Fernandez et al. (2013) and Frasson et al. (2017). It should
be mentioned that the SWOT measurements are also subject to systematic and spatially correlated errors
introduced by different factors (e.g., satellite roll, wet troposphere, etc.). These types of errors are not
addressed in the present paper and will be subject to future investigation.
RiverObs provides the standard deviation of pixel heights associated to each node, at each SWOT overpass.
No information about correlation between noise in different pixels is currently available, thus we cannot
properly define the nodes-associated observation error covariance matrix, which is actually required for a
classical DA algorithm. This is why a simplified representation of the observation error covariance matrix
must be used. Let us note that in our DA scheme, presented in section 3.2, it is sufficient to know the stan-
dard deviation at nodes up to a scaling factor. That is why we can rely on the standard deviation of the

OUBANAS ET AL. 2411

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

Figure 6. Plots A and B show the river width measured at the nodes for the Sacramento River using the simulated SWOT passes 0014 (9 revisits) and 0292 (9 revis-
its), respectively. Plots C and D show the water surface elevation at nodes produced from the passes 0014 and 0292, respectively.

pixel. Moreover, if it is assumed a constant for all nodes, its value does not affect the result of DA. When
more detailed information (correlation between pixels/nodes) are available, we will introduce the observa-
tion covariance matrix that handles the temporal and spatial variability of the error.
Figure 7 presents an average in space for each overpass during the study period, for the Po and the Sacra-
mento Rivers. The averaged standard deviation in space and time is about rz 5 2.5 m in the case of the Po
River, while the averaged uncertainty is higher in the case of the Sacramento River, r 5 4.5 m. In fact, the
observations of the pass 0292 are located in the far range of the swath and have higher errors due to the
low signal return. The corresponding errors are not representative of the observation uncertainty, therefore,
we only consider the average standard deviation associated with the observations of the pass 0014, which
is about rz 5 2.5 m.

Figure 7. The standard deviation associated to the simulated WSE averaged over the RiverObs nodes for each pass, and
the mean of all the passes, for the (left) Po and the (right) Sacramento Rivers.

OUBANAS ET AL. 2412

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

3. Methodology
3.1. Hydraulic Model SIC2
The 1.5 D hydraulic model SIC2, based on the full Saint-Venant equations, has been under development at
IRSTEA for about 30 years, succeeding the former CEMAGREF hydraulic models (Talweg-Fluvia-Sirene)
([Link] The longitudinal and transversal hydraulic effects are described at a finite number of
cross sections along the reach using the compound channel bathymetry representation. A simplified repre-
sentation of the out-of-bank flow involves connected storage areas. For each cross section, the hydraulic
variables such as the wetted area A(z, pg), the wetted perimeter P(z, pg), the hydraulic radius R(z, pg), and the
top width L(z, pg) are computed, for a given water surface elevation z. The pg(x) refers to the parameters
which define the geometry of the corresponding computational cross section, i.e., bottom width l, bank
slope b, and the bed elevation Zb with respect to a chosen reference level, in the case of trapezoidal approx-
imation. Note that SIC2 can handle any type of cross sections including irregular ones. Later on, when talk-
ing about the bathymetry estimation we mean the estimation of the distributed vector Zb(x).
The flow dynamics in the longitudinal direction x at time t are described by the Saint-Venant equations:
@A @q
1 5QL ; (1)
@t @x
@q @q2 =A @z
1 1gA 52gASf 1CL QL v; (2)
@t @x @x
where q(x, t) is the local discharge, z(x, t) is the WSE, QL(x, t) is the lateral discharge, CL ðxÞ 2 ½0; 1 is the lateral
discharge coefficient, v(x, t) 5 q/A is the mean longitudinal velocity, and Sf is the friction slope defined as

q2
Sf 5 ; (3)
Ks2 A2 R4=3

which is related to the Strickler coefficient Ks(x) (inverse of the Manning coefficient) by equation (3). The ini-
tial state (z0, q0) and the boundary conditions (i.e., the inflow discharge Q(t) at the upstream node and a rat-
ing curve, defined by the rating curve parameters prc, at the downstream node) are needed to solve the
equations (1)–(3). Note that the SIC2 model supports different types of boundary conditions; here we
present only those relevant to the chosen test cases. The four-point implicit finite-difference method,
called the Preissmann scheme (Cunge et al., 1980; Novak et al., 2010), is utilized for discretizing the problem
(1)–(3). The fixed-point iterations are used to resolve nonlinearity. For more details on the solver, refer to
Malaterre et al. (2014).
In order to use the hydraulic model SIC2 in a chosen river system, some inputs need to be provided. These
include the cross-sectional geometrical parameters, roughness coefficients, upstream, and downstream
boundary conditions. In the framework of ungauged basins, we propose a methodology, presented in sec-
tion 3.3, to generate initial approximations of these inputs from SWOT observations (WSE, width, and slope)
and globally available ancillary information. From these inputs, the model SIC2 simulates the WSE and the
local discharge fields. Next, we solve the corresponding inverse problem using variational DA to improve
approximations of the inputs of interest.

3.2. Variational Data Assimilation

Dynamical systems can be described by a numerical model, here denoted M : U ! X , which maps the
model inputs U 2 U, also called the ‘‘control vector,’’ into the model state X 2 X:
MðUÞ5X: (4)

The model state X is assumed to be related to observations via an observation operator H : X ! Y, where
Y is the observation space:
Y5HðXÞ 2 Y: (5)

Therefore, the control-to-observation nonlinear mapping G : U ! Y can be defined as:

HðXÞ5HðMðUÞÞ : 5GðUÞ5Y: (6)

OUBANAS ET AL. 2413

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

Both observations Y and model inputs U embody uncertainties no and nb that arise from various sources
(e.g., instrumental noise, parameterization, discretization, etc.) and should be taken into account:

Y  5Y t 1no ; (7)

where Y t 5GðUt Þ is the ‘‘true’’ observation vector. And:

Ub 5Ut 1nb ; (8)

where Ut is the true model input vector and Ub is its best available first guess, known as a ‘‘background’’ in
variational DA or as a ‘‘prior’’ in Bayesian statistics.
The full input vector of the hydraulic model presented above consists of the following variables:
 T
U5 z0 ; q0 ; Q; prc ; QL ; CL ; Ks ; Zb ; pg ; pnm ; (9)

where pnm are the numerical scheme parameters. For a given U, we obtain the flow fields by solving the
model equations (1)–(3) such that:
ðz; qÞ5fðzðxi ; tÞ; qðxi ; tÞÞ; i51; . . . ; Ng; t 2 ½0; T: (10)

In certain components of the input vector V  U, the uncertainty is significant and strongly affects the
model predictions. The aim of data assimilation is to estimate V using observations Y (i.e., to improve the
first guess Vb), whereas the remaining components U0 5U n V are fixed at their background value Ub0 .
In the present study, the inputs of interest are:
V5ðQðtÞ; Zb ðxÞ; KS ðxÞÞ; (11)

and the assimilated observations are:

Y5fzðxi ; tj Þ; i51; . . . ; NN ; j51; . . . ; NP g; (12)

where NN and NP are, respectively, the numbers of RiverObs nodes and SWOT overpasses. In this case, the
dynamical model M refers to the code of the hydraulic model SIC2, while H is an additional module repre-
senting the observation operator.
The variational DA method gives the best estimate of the control vector by minimizing a cost function.
Assuming that errors associated with the observation and the background information are Gaussian, i.e., nb
 Nð0; BÞ and no  Nð0; RÞ, where B and R are the corresponding covariance matrices, the conventional for-
mulation of the DA problem is as follows:

V^ 5arg min JðVÞ; (13)

V

where
1 1
JðVÞ5 k R21=2 ðGðV; Ub0 Þ2Y  Þk2 1 k B21=2 ðV2Vb Þk2 : (14)
2 2
A variant of the conventional variational DA method useful for the estimation of the model variables and
parameters affected by uncertainties in hydraulic applications is described in detail in Gejadze and Mala-
terre (2017). First, instead of (14), we consider a modified cost function:
1 a
JðV; aÞ5 k R21=2 ðGðV; Ub0 Þ2Y  Þk2 1 k B21=2 ðV2Vb Þk2 ; (15)
2 2
where a > 0 is a regularization parameter. This is done to reduce the impact of possible errors in assigning
B. The cost function (15) is used in the Tikhonov regularization method (Tikhonov et al., 1977). Second, we
consider the change of variables V5Vb 1B1=2 W, in which case the above cost function takes the form:
1 a
JðW; aÞ5 k R21=2 ðGðVb 1B1=2 W; Ub0 Þ2Y  Þk2 1 k Wk2 : (16)
2 2
If the control vector V is composed of heterogeneous components (i.e., flow variables, physical parameters, and
parameters describing the domain geometry), then the corresponding parts of the gradient J 0V could have a

OUBANAS ET AL. 2414

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

very different norm. Since G is a nonlinear operator, this may compromise the robustness of the minimization
process. The change of variables helps to assign the appropriate weights to the gradient, see Gejadze and
Malaterre (2017). This is the advantage of the latter formulation as compared to the classical one where at first
iterations the contribution of the background term is negligible (zero at first iteration). Let us mention that this
change of variables serves as a preconditioning in the framework of the incremental approach of variational
DA. To implement (16) one needs B1/2 instead of B21 (or B21/2) in the original formulation (14). This is obtained
using the approach presented in Gejadze and Malaterre (2017), based on the assumption that the control vari-
able is a distributed function of space or time which belongs to the Sobolev space of the second order. The
square root of the matrix B is obtained using Cholesky decomposition. This approach can be performed using
modest computational resources in the case of one-dimensional hydraulic problems.
The optimal choice of the regularization parameter a is a key issue in the Tikhonov regularization method.
For example, a can be chosen from the residual principle as follows:
^ Þ  v2 ðM; dÞ;
JðW (17)

where M and d are the observation space dimension and the confidence level, respectively. This implies
solving the minimization problem for different values of a to satisfy the condition (17). This could be com-
putationally expensive. An alternative is to use the iterative regularization method (Kaltenbacher et al.,
2008), in which the cost function to be minimized is reduced to its residual term:
1
JðWÞ5 k R21=2 ðGðVb 1B1=2 W; Ub0 Þ2Y  Þk2 : (18)
2
This should be combined with a stopping criterion in order to obtain a regularized solution by early termi-
nation of the iterations. We employ the residual principle (17) to stop the iterations in the iterative form:
^ i Þ  v2 ðM; dÞ;
JðW (19)

where i is the iteration number. The equivalence of these two approaches has been proved in Kaltenbacher
et al. (2008) for a class of the so-called ‘‘regular’’ iterative methods which include the steepest descent, con-
jugate gradient, BFGS or L-BFGS. Let us underline that the equivalence is only valid for the cost function in
the form (16). Note that the residual principle can only be used in case the statistical properties of the obser-
vation error are well known. Otherwise, one should use different approaches such as L-curve (Hansen &
O’Leary, 1993), cross-validation (Golub et al., 1979), etc.
The update step of the L-BFGS algorithm reads as follows:

~ 21 B1=2 J 0 ðVi Þ; W0 50;

Wi11 5Wi 1bi H (20)
i

Vi11 5Vb 1B1=2 Wi11 ; V0 5Vb ; (21)

where ~ 21
H is the approximated inverse of the Hessian built by the algorithm. The gradient of the cost func-
i
tion is given by:

J 0 ðVi Þ5ðG0 ðVi ÞÞ R21 ðGðVi ; Ub0 Þ2YÞ; (22)


where ðG0 ðVÞ and ðG0 ðVÞÞ are, respectively, the tangent linear and adjoint counterparts of the nonlinear
operator ðGðV; Ub0 Þ given by the following Gateaux derivative:

GðV1tw; Ub0 Þ2GðV; Ub0 Þ

G0 ðVÞw5 lim ; (23)
t!0 t
hw; ðG0 ðVÞÞ w  iU 5hG0 ðVÞw; w  iY ; 8w; w  : (24)

The tangent linear and adjoint operators associated with the SIC2 model have been produced using the
automatic differentiation tool TAPENADE developed at INRIA (Hasco€et & Pascual, 2004).
The gradient-based minimization involves the use of the direct and adjoint models. While the direct model
represents the forward time integration, the adjoint model is regarded as the corresponding backward time
integration. For the latter, the system trajectory needs to be saved at each forward integration time step.
Thus, for a long assimilation window (e.g., a few months) and a time step consistent with the Preissmann

OUBANAS ET AL. 2415

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

numerical scheme (e.g., 20 min), the process memory requirements become prohibitively expensive. To
overcome this limitation, we consider a sequential implementation of the variational DA approach which
operates with assimilation subwindows.
Let us introduce a time lag dT consistent with the characteristic time of the dynamical system. At the start
of the DA process, i.e., for the first assimilation subwindow, the initial condition is given by the steady state
flow solution consistent with the initial guess, i.e., QðtÞ5Qb ðt0 Þ5const. For any subsequent subwindow, the
initial condition is given by the estimated state at time T2dT from the previous assimilation subwindow,
^ ðx; T2dTÞ, whereas the background is given by the estimated discharge
i.e., z0 ðxÞ5^z ðx; T2dTÞ and q0 ðxÞ5q
^
at time T2dT, i.e., Qb ðtÞ5QðT2dTÞ5const. This sequential implementation is an additional option to the
method presented in Gejadze and Malaterre (2017). The conducted study is only possible using the sequen-
tial version of DA, when using modest computational resources, i.e., a laptop with 16–32 GB of RAM. Note
that all operational variational DA systems are based on sequential (cyclic) technologies.

3.3. Computing First Guess Model Inputs

The optimization method requires a first guess on the variables of interest Q(t), Zb(x), and KS(x). Therefore,
we suggest a way of generating this initial background using the SWOT observations only and globally
available background information applicable to ungauged basins. Note that the background values of Q, Zb,
and KS, generated using this method, are only introduced during the first DA subwindow then will be later
improved during the data assimilation process.
This method is based on a prior knowledge of discharge Qb, which for the current test cases was taken
as an output of the global Water Balance Model (WBM), that uses monthly atmospheric forcing data
(mean air temperature and precipitation), obtained from the Climatic Research Unit (CRU, [Link]
[Link]/) (Wisser et al., 2010), and SWOT observations of the water surface elevation z, the top
@z
width L, and the slope S5 @x . The first guesses on A and KS are calculated using the Manning’s equation,
assuming trapezoidal cross sections at each river node. The WBM-based prior discharge Qb is 841.8 m3
s21 throughout the study area for the Po River and 377 m3 s21 for the Sacramento River. Note that any
available first guess on discharge, although not reliable, can be used as long as the resulting dynamics
are supported by SIC2 model. As we aim to apply the presented methodology to ungauged basins, no in
situ gauge data were considered.
An approximate bathymetry is then built by estimating the bottom width l, the depth h, and the bank slope
b, assuming trapezoidal cross sections. To do so, let us consider the discharge formula derived from the
2 1
Manning equation QM 5KS R3 S2 A, where the index M refers to the ‘‘Manning,’’ the wetted perimeter PM can
be then expressed as:
 232
1 1
PM 5A Qb S22 A21 : (25)
KS
Therefore, the top width L, the area A and the wetted perimeter P are given by:
L5l12bh; (26)
h
A5 ðL1lÞ; (27)
2
pffiffiffiffiffiffiffiffiffiffiffi
P5l12h 11b2 : (28)

From the above formulas, P can be written as a function of l as:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
  
2A 2 L2l 2
PðlÞ5l12 1 : (29)
l1L 2
The solution to this problem was computed using an iterative numerical method to find the estimates of l, h,
and b. Finally, the bed elevation is retrieved by subtracting the depth h from the h from the WSE z; Zb 5 z 2 h.

4. Experimental Design
The challenge of this study stems from different factors such as low data accuracy and temporal frequency,
unknown river bathymetry and bed roughness, the need for estimating a heterogeneous control vector,

OUBANAS ET AL. 2416

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

etc. Therefore, our aim is to demonstrate the performance and robustness of the proposed variational DA
method for discharge estimation under uncertainties from the simulated SWOT observations, with no in
situ information. Since the long-term hydraulic behavior is governed by the boundary conditions and the
lateral inflow, we focus on estimating the upstream inflow discharge hydrograph Q(t). Note that lateral
inflow discharge QL(t) is not considered in this study. Future work will investigate more complex hydraulic
configurations involving tributaries and their interactions with the main channel; a configuration that is sup-
ported by the hydraulic model SIC2.
We perform simultaneous estimation of the inflow discharge Q(t), the bed level Zb(x), and the Strickler coef-
ficient KS(x) at the nodes scale, i.e., at every 200 m. The background information is subjected to bias, which
is the systematic part of the original uncertainty in the model inputs. Hence, the corresponding discharge
estimation error can be arbitrarily large. It is the purpose of DA to remove/reduce this uncertainty. This is
exactly the reason why Zb and KS are included into the control vector V. Removing the bias in observations
is more tricky and requires a special analysis of the residuals and their derivatives in a series of independent
experiments (reanalysis) (Dee & Da Silva, 1999; Desroziers et al., 2005). This is not implemented in the cur-
rent version of the algorithm.
In order to deal with long study periods, the sequential version of the variational DA method presented in
section 3 is considered. The size of the subwindow can be arbitrary as long as the required computational
resources are available (long subwindow needs more memory). The priors of discharge, river bathymetry,
and roughness coefficient, during the first DA subwindow, are taken as described in section 3.3. For subse-
quent subwindows, the prior of discharge is taken as the final time minus time-lag discharge estimate from
the previous subwindow according to the sequential approach described in section 3.2. In order to verify
the usefulness of the combination of the estimates of Zb(x) and KS(x), these variables are simultaneously esti-
mated with Q(t) during the first subwindow only then are fixed at their estimated value for subsequent
subwindows.
The estimation quality is measured based on different error metrics:

 Root Mean Square Error (RMSE) defined by:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1 X ^ 2
RMSE½Q5 QðtÞ2Qt ðtÞ ; (30)
T t

 relative Root Mean Square Error (rRMSE):

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
!2ffi
u
u1 X QðtÞ2Q ^ t
ðtÞ
rRMSE½Q5 t ; (31)
T t
Qt ðtÞ

 Normalized Root Mean Square Error (NRMSE):

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1 1 X ^ 2
NRMSE½Q5 QðtÞ2Qt ðtÞ ; (32)
Qt T t

 Nash-Sutcliffe Efficiency (NSE):

X 2
^
QðtÞ2Q t
ðtÞ t
NSE½Q512 X  2 ; (33)
t
Qt 2Qt ðtÞ

 Volumetric Efficiency (VE):

X
^
jQðtÞ2Qt
ðtÞj
tX
VE½Q512 ; (34)
t
Qt ðtÞ

OUBANAS ET AL. 2417

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

Figure 8. The Po River discharge estimation at Borgoforte station (upper) together with the (lower left) Strickler coefficient and (lower right) bed level along the study
area.

^
where Qt(t) and Qt are the ‘‘true’’ (reference) inflow discharge and its time-average value, respectively, QðtÞ
is the inflow discharge estimate and T is the full time-length of the experiment.

5. Results and Discussion

Despite the highly uncertain first guess on the Po River discharge, with a relative Root Mean Square Error of
rRMSE½Qb 581:1%, simultaneous estimation of the upstream discharge Q(t), the bed level Zb(x) and the
Strickler coefficient KS(x) was successfully performed using eight subwindows 42 days wide. The comparison
between the daily discharge estimated with our method and the assumed true discharge from the Po River
hydraulic model are presented in Figure 8. Summary statistics for the estimated daily discharge and the esti-
mated discharge strictly during the SWOT overpasses are presented in Table 1, using the error metrics intro-
duced in section 4. Discharge errors were lower at the time of the SWOT observations, amounting to a

Table 1
Error Metrics of Discharge Estimation for the Po River During the Full Assimilation Windows for Daily Discretization and
Irregular Observation Sampling

Computational DT RMSE (m3 s21) rRMSE (%) NRMSE (%) NSE VE

Daily 657.9 29.8 36.5 0.54 0.77
Irregular DTobs 221.4 12.6 12.1 0.95 0.91

OUBANAS ET AL. 2418

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

Table 2
Error Metrics of Discharge Estimation for the Sacramento River During the Full Assimilation Windows for Daily Discretization
and Irregular Observation Sampling

Computational DT RMSE (m3 s21) rRMSE (%) NRMSE (%) NSE VE

Daily 167.2 21.2 75.0 0.18 0.75
Irregular DTobs 22.6 11.2 12.0 0.87 0.91

relative Root Mean Square Error (rRMSE) of 12.1%. The estimate of discharge at the daily scale increases
rRMSE to 29.8% due to unobserved dynamics between the SWOT overpasses.
Similarly, the first guess on the Sacramento River discharge had high uncertainty, amounting to a relative
error of rRMSE[Qb] 5 84.4% (Table 2). Here, two assimilation subwindows of 63 and 105 days are considered
for DA. Figure 9 shows the estimates of Q(t), Zb(x), and KS(x) over the Sacramento study area, during the 6
month study period. The flow dynamics at the Sacramento River change are significantly faster than at the
Po River, which exemplifies how the limited SWOT temporal resolution might affect the understanding of
narrower rivers. As the Sacramento River study area is only revisited twice per SWOT cycle and due to the
faster discharge dynamics, the flood peaks were not sampled, with one of the most significant flood events
having no observations. For example, the maximum discharge 1541.97 m3 s21 was recorded on 17 February,
after the 13 February overpass. By the time SWOT returned, on 22 February, the flood wave had already left the

Figure 9. The Sacramento River discharge estimation at Hamilton station (upper) together with the (lower left) Strickler coefficient and (lower right) bed level
along the study area.

OUBANAS ET AL. 2419

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

area. Nevertheless, discharge was successfully estimated, with rRMSE computed solely during SWOT overpasses
being as low as 11.2%, whereas daily discharge rRMSE reached a higher value of 21.2% (Table 2). The increased
rRMSE were driven by the unobserved higher discharge events, however, the model showed considerable skill
estimating the preflood flows (before day of year 40), despite the errors in the discharge background.
The Sacramento River example illustrates the fact that the characteristic time of the hydraulic system, i.e.,
the dynamical transition time between two possible equilibrium states, should be consistent with the tem-
poral frequency of observations in order to ensure that the river dynamics between satellite overpasses can
be accurately recovered. The characteristic time of the Po and Sacramento study Rivers under consideration
are 6 and 5 days, respectively. Therefore, longer study areas may be more suitable for this type of sparse
temporal sampling. Moreover, the river orientation with respect to the satellite ground tracks has a key role
in the temporal and spatial observation sampling, i.e., across-track rivers have a higher temporal frequency
(the Po River) while the along-track rivers are observed with better spatial distribution (the Sacramento
River). The temporal and spatial gaps between the SWOT observations can be completed by supplementary
sources of data such as the virtual stations of the other existing satellite missions.
The Strickler coefficient KS and the bed level Zb may be locally improved during the minimization process
(Figures 8 and 9). However, this combination of optimal estimates (Zb, Ks) only allow accurate estimation of
discharge Q in both study areas despite the use of highly uncertain background information. This behavior
is typical of ill-posed estimation problems (the so-called ‘‘equifinality issue’’) (Gejadze & Malaterre, 2016).
Thus, the estimates of Zb and KS may not provide a reliable global information of river bathymetry and bed
roughness.
The performance of the variational DA method is illustrated in Figure 10. The cost function J and the norm
of its gradient jjrJjj are presented with respect to the number of iterations at each subwindow, for the Po
and the Sacramento Rivers. The minimization algorithm has converged within less than 10 iterations, i.e., 5
iterations for the Po River, while few more iterations, 6–8, were required for the Sacramento River case.

1800
1600
1400
Gradient norm

1200
1000
800
600
400
200
0
0 2 4 6 8 10 12 14 16
Iteration number

Figure 10. The minimization process for the Po River (upper plots) and the Sacramento River (lower plots) study cases. (left) The cost function J and (right) the
norm of the gradient jjrJjj.

OUBANAS ET AL. 2420

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

Note that in the first subwindow, the three variables Q, Zb, and KS are estimated simultaneously while during
subsequent subwindows, only Q is estimated which explains why fewer iterations can be required for the
convergence.
In case the statistical properties of the observation error are not very well known, the iterations are stopped
using the criterion Ji11 2Ji < , where  is a threshold equal to 1024, which is a simplified implementation of
the L-curve criterion. Moreover, Figure 10 indicates the actual level of the observation noise standard devia-
tion, which is significantly lower than the one predicted by RiverObs. This level can be roughly assessed
from the residual principle:

r 2z =r2z ;
J  m~ (35)

where m is the total number of observation points and rz and r ~ z are the estimated standard deviations
from RiverObs and simplified L-curve, respectively. For example, the first subwindow of data assimilation
has a cost function equal to J 5 122.3 at the end of the iterations, in the case of the Po River. The total num-
ber of observations during the subwindow is m 5 2,938, while the initial estimation of the standard devia-
tions is rz 5 2.5 m. Therefore, following the residual principle (35), the estimated actual error standard
deviation is about r ~ z 50:5 m. Future work will include simultaneous estimation of the model inputs
together with parameters of the observation noise.

6. Conclusions
Our work presented a methodology suitable for the estimation of discharge on ungauged basins using
solely remote sensing information from platforms such as the upcoming SWOT satellite. The present article
describes necessary modifications to the classical variational data assimilation method applied to the full
Saint-Venant-based hydraulic model SIC2. Moreover, the method offers the flexibility of also assimilating in
situ data in addition to satellite observations when those are available. The implemented sequential version
enables better use of computer resources, overcoming memory limitations, which hindered the applicability
of the classical data assimilation methods to long-time periods.
We demonstrated the applicability of the proposed methodology using synthetic SWOT overpasses gener-
ated with the SWOT simulator developed at the JPL over the Sacramento and the Po Rivers. In our two
study cases, discharge was successfully estimated at the time of the overpasses, with rRMSE of 12.1% and
11.2% for the Po and the Sacramento Rivers, respectively. The estimated river discharge at a daily scale
increases the rRMSE to 29.8% and 21.2%, respectively. The higher degradation of the model skill for the Sac-
ramento River when estimating daily discharge happened mostly due to poor temporal sampling during
the observed flood waves, showing the importance of having temporal sampling that is compatible with
the characteristic time of the hydraulic system.
The presented method demonstrates good accuracy and robustness, while requiring modest computational
resources. In both cases, the method converged after less than 10 iterations when estimating Q, Zb, and KS
simultaneously during the first subwindow, while fewer iterations (2–5) were required to estimate discharge
during the subsequent subwindows. Future work will report multimission variational data assimilation
including the existing satellites missions such as JASON, ENVISAT, and Sentinel in addition to the future
Acknowledgments SWOT platform.
The presented work was undertaken
with the financial support of IRSTEA
and Collecte Localisation Satellite (CLS) References
and takes part of the PhD of Hind
Oubanas. The authors would like to Alsdorf, D. E., & Lettenmaier, D. P. (2003). Tracking fresh water from space. Science, 301(5639), 1491–1494.
thank Dr. Franck Mercier, CLS, for his Andreadis, K. M., Clark, E. A., Lettenmaier, D. P., & Alsdorf, D. E. (2007). Prospects for river discharge and depth estimation through assimila-
implication in the project. The tion of swath-altimetry into a raster-based hydrodynamics model. Geophysical Research Letters, 34, L10403. [Link]
experimental data used in this study 2007GL029721
are attached as supporting Bartsch, A., Wagner, W., Scipal, K., Pathe, C., Sabel, D., & Wolski, P. (2009). Global monitoring of wetlands—The value of ENVISAT ASAR
information files ‘‘[Link].’’ global mode. Journal of Environmental Management, 90(7), 2226–2233.
Complementary data can be Benke, A. C., & Cushing, C. E. (2011). Rivers of North America. Cambridge, MA: Academic Press.
requested from Dr. Alessio Biancamaria, S., Andreadis, K. M., Durand, M., Clark, E. A., Rodriguez, E., Mognard, N. M., et al. (2010). Preliminary characterization of SWOT
Domeneghetti ([Link]@ hydrology error budget and global capabilities. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 3(1), 6–19.
[Link]) and Rui Wei (wei.263@osu. Biancamaria, S., Durand, M., Andreadis, K., Bates, P., Boone, A., Mognard, N., et al. (2011b). Assimilation of virtual wide swath altimetry to
edu). improve arctic river modeling. Remote Sensing of Environment, 115(2), 373–381.

OUBANAS ET AL. 2421

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

Biancamaria, S., Hossain, F., & Lettenmaier, D. (2011a). Forecasting transboundary river water elevations from space. Geophysical Research
Letters, 38, L11401. [Link]
Biancamaria, S., Lettenmaier, D. P., & Pavelsky, T. M. (2016). The SWOT mission and its capabilities for land hydrology. Surveys in Geophysics,
37(2), 307–337.
Birkett, C. (1995). The contribution of TOPEX/Poseidon to the global monitoring of climatically sensitive lakes. Journal of Geophysical
Research, 100(C12), 25179–25204.
Birkett, C. M. (1998). Contribution of the TOPEX NASA radar altimeter to the global monitoring of large rivers and wetlands. Water Resources
Research, 34(5), 1223–1239.
Buer, K., Forwalter, D., Kissel, M., & Stohlert, B. (1989). The middle Sacramento River: Human impacts on physical and ecological processes
along a meandering river. In Abell, D. L. (Ed.), Proceedings of the California Riparian Systems Conference: Protection, management, and res-
toration for the 1990s (Gen. Tech. Rep. PSW-GTR-110, pp. 22–32). Berkeley, CA: Pacific Southwest Forest and Range Experiment Station,
Forest Service, U.S. Department of Agriculture.
Calmant, S., & Seyler, F. (2006). Continental surface waters from satellite altimetry. Comptes Rendus Geoscience, 338(14), 1113–1122.
Castellarin, A., Domeneghetti, A., & Brath, A. (2011). Identifying robust large-scale flood risk mitigation strategies: A quasi-2D hydraulic
model as a tool for the Po River. Physics and Chemistry of the Earth, Parts A/B/C, 36(7), 299–308.
Chow, V. T. (1964). Handbook of applied hydrology. New York, NY: McGraw-Hill Book Company.
Courtier, P., Andersson, E., Heckley, W., Vasiljevic, D., Hamrud, M., Hollingsworth, A., et al. (1998). The ECMWF implementation of three-
dimensional variational assimilation (3D-var). I: Formulation. Quarterly Journal of the Royal Meteorological Society, 124(550), 1783–1807.
Cretaux, J.-F., Berge-Nguyen, M., Leblanc, M., Del Rio, R. A., Delclaux, F., Mognard, N., et al. (2011). Flood mapping inferred from remote
sensing data. International Water Technology Journal, 1, 48–62.
Cr
etaux, J.-F., & Birkett, C. (2006). Lake studies from satellite radar altimetry. Comptes Rendus Geoscience, 338(14), 1098–1112.
Cuchi, K. (1986). Multi-look processing of synthetic aperture radar data from dynamic ocean surfaces. Pattern Recognition Letters, 4(4),
305–314.
Cunge, J. A., Holly, F. M., & Verwey, A. (1980). Practical aspects of computational river hydraulics. Ann Arbor, MI: University of Michigan
Press.
de Oliveira Campos, I., Mercier, F., Maheu, C., Cochonneau, G., Kosuth, P., Blitzkow, D., et al. (2001). Temporal variations of river basin waters
from TOPEX/Poseidon satellite altimetry: Application to the amazon basin. Comptes Rendus de l’Academie des Sciences, 333(10), 633–643.
Dee, D. P., & Da Silva, A. M. (1999). Maximum-likelihood estimation of forecast and observation error covariance parameters. Part I: Method-
ology. Monthly Weather Review, 127(8), 1822–1834.
Desroziers, G., Berre, L., Chapnik, B., & Poli, P. (2005). Diagnosis of observation, background and analysis-error statistics in observation
space. Quarterly Journal of the Royal Meteorological Society, 131(613), 3385–3396.
Ding, Y., & Wang, S. S. (2012). Optimal control of flood diversion in watershed using nonlinear optimization. Advances in Water Resources, 44,
30–48.
Domeneghetti, A., Tarpanelli, A., Brocca, L., Barbetta, S., Moramarco, T., Castellarin, A., et al. (2014). The use of remote sensing-derived water
surface data for hydraulic model calibration. Remote Sensing of Environment, 149, 130–141.
Durand, M., Andreadis, K. M., Alsdorf, D. E., Lettenmaier, D. P., Moller, D., & Wilson, M. (2008). Estimation of bathymetric depth and slope
from data assimilation of swath altimetry into a hydrodynamic model. Geophysical Research Letters, 35, L20401. [Link]
2008GL034150
Durand, M., Neal, J., Rodrıguez, E., Andreadis, K. M., Smith, L. C., & Yoon, Y. (2014). Estimating reach-averaged discharge for the River Severn
from measurements of river water surface elevation and slope. Journal of Hydrology, 511, 92–104.
Durand, M., Gleason, C., Garambois, P.-A., Bjerklie, D., Smith, L., Roux, H., et al. (2016). An intercomparison of remote sensing river discharge
estimation algorithms from measurements of river height, width, and slope. Water Resources Research, 52, 4527–4549. [Link]
10.1002/2015WR018434
Fekete, B. M., & V€or€osmarty, C. J. (2002). The current status of global river discharge monitoring and potential new technologies comple-
menting traditional discharge measurements. In Predictions in ungauged basins: PUB kick-off (proceedings of the PUB kick-off meeting
held in Brasilia, 20–22 November 2002) (IAHS Publ. 349). London, UK: International Association of Hydrological Sciences.
Fernandez, D. E., Pollard, B., & Vaze, P. (2013). SWOT mission performance error budget: A revision (Tech. Rep. JPL Doc. D-79804). Jet Pro-
pulsion Laboratory.
Fischer, C., Montmerle, T., Berre, L., Auger, L., & Ştefanescu, S. E. (2005). An overview of the variational assimilation in the ALADIN/France
numerical weather-prediction system. Quarterly Journal of the Royal Meteorological Society, 131(613), 3477–3492.
Fjortoft, R., Gaudin, J.-M., Pourthie, N., Lalaurie, J.-C., Mallet, A., Nouvel, J.-F., et al. (2014). Karin on SWOT: Characteristics of near-nadir Ka-
band interferometric SAR imagery. IEEE Transactions on Geoscience and Remote Sensing, 52(4), 2172–2185.
Frappart, F., Calmant, S., Cauhop e, M., Seyler, F., & Cazenave, A. (2006). Preliminary results of ENVISAT RA-2-derived water levels validation
over the amazon basin. Remote Sensing of Environment, 100(2), 252–264.
Frasson, R. P. D. M., Wei, R., Durand, M., Minear, J. T., Domeneghetti, A., Schumann, G., et al. (2017). Automated river reach definition strate-
gies: Applications for the surface water and ocean topography mission. Water Resources Research, 53, 8164–8186. [Link]
1002/2017WR020887
Garambois, P.-A., & Monnier, J. (2015). Inference of effective river properties from remotely sensed observations of water surface. Advances
in Water Resources, 79, 103–120.
Gauthier, P., Charette, C., Fillion, L., Koclas, P., & Laroche, S. (1999). Implementation of a 3D variational data assimilation system at the Cana-
dian Meteorological Centre. Part I: The global analysis. Atmosphere-Ocean, 37(2), 103–156.
Gauthier, P., Tanguay, M., Laroche, S., Pellerin, S., & Morneau, J. (2007). Extension of 3Dvar to 4Dvar: Implementation of 4Dvar at the meteo-
rological service of Canada. Monthly Weather Review, 135(6), 2339–2354.
Gejadze, I., & Malaterre, P.-O. (2017). Discharge estimation under uncertainty using variational methods with application to the full Saint-
Venant hydraulic network model. International Journal for Numerical Methods in Fluids, 83(5), 405–430.
Gejadze, I. Y., & Malaterre, P.-O. (2016). Design of the control set in the framework of variational data assimilation. Journal of Computational
Physics, 325, 358–379.
Golub, G. H., Heath, M., & Wahba, G. (1979). Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics,
21(2), 215–223.
Hasco€et, L., & Pascual, V. (2004). Tapenade 2.1 user’s guide (technical report RT-0300, 78 pp.). France: INRIA.
Hansen, P. C., & O’leary, D. P. (1993). The use of the l-curve in the regularization of discrete ill-posed problems. SIAM Journal on Scientific
Computing, 14(6), 1487–1503.

OUBANAS ET AL. 2422

 19447973, 2018, 3, Downloaded from [Link] by National Institutes Of Health Malaysia, Wiley Online Library on [22/09/2024]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Water Resources Research 10.1002/2017WR021735

Hossain, F., Siddique-E-Akbor, A., Mazumder, L. C., ShahNewaz, S. M., Biancamaria, S., Lee, H., et al. (2014). Proof of concept of an altimeter-
based river forecasting system for transboundary flow inside Bangladesh. IEEE Journal of Selected Topics in Applied Earth Observations
and Remote Sensing, 7(2), 587–601.
Hostache, R., Lai, X., Monnier, J., & Puech, C. (2010). Assimilation of spatially distributed water levels into a shallow-water flood model. Part
II: Use of a remote sensing image of Mosel River. Journal of Hydrology, 390(3), 257–268.
Kaltenbacher, B., Neubauer, A., & Scherzer, O. (2008). Iterative regularization methods for nonlinear ill-posed problems (Vol. 6). Berlin, Ger-
many: Walter de Gruyter.
Kouraev, A. V., Zakharova, E. A., Samain, O., Mognard, N. M., & Cazenave, A. (2004). Ob’river discharge from TOPEX/Poseidon satellite altime-
try (1992–2002). Remote Sensing of Environment, 93(1), 238–245.
Lai, X., & Monnier, J. (2009). Assimilation of spatially distributed water levels into a shallow-water flood model. Part I: Mathematical method
and test case. Journal of Hydrology, 377(1), 1–11.
Malaterre, P., Baume, J., & Dorchies, D. (2014). Simulation and integration of control for canals software (SIC2), for the design and verifica-
tion of manual or automatic controllers for irrigation canals. In USCID conference on planning, operation and automation of irrigation
delivery systems (pp. 377–382). Arizona, AZ: Phoenis.
Medina, C. E., Gomez-Enri, J., Alonso, J. J., & Villares, P. (2008). Water level fluctuations derived from ENVISAT Radar altimeter (RA-2) and in-
situ measurements in a subtropical waterbody: Lake Izabal (Guatemala). Remote Sensing of Environment, 112(9), 3604–3617.
Montanari, A., Ceola, S., Baratti, E., Domeneghetti, A., & Brath, A. (2017). Chapter 116 - Po River Basin. In Handbook of applied hydrology
(2nd Edn., pp. 1–4). New York City, NY: McGraw-Hill Education.
Novak, P., Guinot, V., Jeffrey, A., & Reeve, D. E. (2010). Hydraulic modelling—An introduction: Principles, methods and applications. Boca
Raton, FL: CRC Press.
Oubanas, H., Gejadze, I., Malaterre, P. O., & Mercier, F. (2018). River discharge estimation from synthetic SWOT-type observations using vari-
ational data assimilation and the full Saint-Venant hydraulic model. Journal of Hydrology. [Link]
Papa, F., Legr esy, B., & R
emy, F. (2003). Use of the TOPEX-Poseidon dual-frequency radar altimeter over land surfaces. Remote Sensing of
Environment, 87(2), 136–147.
Papa, F., Prigent, C., Durand, F., & Rossow, W. B. (2006). Correction to ‘‘Wetland dynamics using a suite of satellite observations: A case study
of application and evaluation for the Indian subcontinent’’. Geophysical Research Letters, 33, L15401. [Link]
2006GL026943
Rabier, F., J€arvinen, H., Klinker, E., Mahfouf, J.-F., & Simmons, A. (2000). The ECMWF operational implementation of four-dimensional varia-
tional assimilation. I: Experimental results with simplified physics. Quarterly Journal of the Royal Meteorological Society, 126(564), 1143–
1170.
Rodriguez, E. (2016). SWOT science requirements document (JPL Doc. d-61923). Jet Propulsion Laboratory.
Rogers, W. (2014). Central valley floodplain evaluation and delineation, subtask 5, combined sacramento river system model Rep. Sacra-
mento, CA: California Department of Water Resources.
Sneeuw, N., Lorenz, C., Devaraju, B., Tourian, M. J., Riegger, J., Kunstmann, H., et al. (2014). Estimating runoff using hydro-geodetic
approaches. Surveys in Geophysics, 35(6), 1333–1359.
Tikhonov, A. N., Arsenin, V. I., & John, F. (1977). Solutions of ill-posed problems (Vol. 14). Washington, DC: Winston.
Tockner, K., Uehlinger, U., & Robinson, C. T. (2009). Rivers of Europe. Cambridge, MA: Academic Press.
Tourian, M., Sneeuw, N., & Bardossy, A. (2013). A quantile function approach to discharge estimation from satellite altimetry (ENVISAT).
Water Resources Research, 49, 4174–4186. [Link]
Ulaby, F. T., & Dobson, M. C. (1989). Handbook of radar scattering statistics for terrain (Artech House Remote Sensing Library). Norwood, MA:
Artech House.
Ulaby, F. T., Long, D. G., Blackwell, W., Elachi, J. C., Fung, A., Ruf, K., et al. (2014). Microwave radar and radiometric remote sensing (Vol. 4).
Ann Arbor, MI: University of Michigan Press.
Wisser, D., Fekete, B., V€ or€
osmarty, C., & Schumann, A. (2010). Reconstructing 20th century global hydrography: A contribution to the global
terrestrial network-hydrology (GTN-H). Hydrology and Earth System Sciences, 14(1), 1–24.
Yoon, Y., Durand, M., Merry, C. J., Clark, E. A., Andreadis, K. M., & Alsdorf, D. E. (2012). Estimating river bathymetry from data assimilation of
synthetic SWOT measurements. Journal of Hydrology, 464, 363–375.

OUBANAS ET AL. 2423