Large scale hydrological model river storage and discharge correction using satellite altimetry-based discharge product

. Land surface models (LSMs) are widely used to study the continental part of the water cycle. However, even though their accuracy is increasing, inherent model uncertainties can not be avoided. In the meantime, remotely sensed observations of the continental water cycle variables such as soil moisture, lakes and river elevations are more frequent and accurate. Therefore, those two different types of information can be combined, using data assimilation techniques to reduce a model’s uncertainties in its state variables or/and in its input parameters. The objective of this study is to present a data assimilation platform that assimilates into the large-scale ISBA-CTRIP LSM a punctual river discharge product, derived from ENVISAT nadir altimeter water elevation measurements and rating curves, over the whole Amazon basin. To deal with the scale difference between the model and the observation, the study also presents an initial development for a localization treatment that allows one to limit the impact of observations to areas close to the observation and in the same hydrological network. This assimilation platform is based on the ensemble Kalman ﬁlter and can correct either the CTRIP river water storage or the discharge. Root mean square error (RMSE) compared to gauge discharges is globally reduced until 21 % and at Óbidos, near the outlet, RMSE is reduced by up to 52 % compared to ENVISAT-based discharge. Finally, it is shown that localization improves results along the main tributaries.


Introduction
The continental part of the water cycle is commonly studied, at large scale, with hydrological modelling.These models are generally issued from the coupling of a land surface model (LSM) with a river routing model (RRM).The LSM determines the water and energy budget at the surface by spreading precipitations between the soil and the canopy.Meanwhile, the RRM transfers water mass through the basin to the outlet and gives an estimate of river discharge.
RRMs are mainly based on kinematic (e.g.Oki and Sud, 1998) or sometimes also diffusive wave models (e.g.Yamazaki et al., 2011).Several RRMs have been developed since the 1990s and they differ mainly in their modelling of the flow velocity and the inclusion or not of groundwater and floodplain dynamics (Oki and Sud, 1998;Coe, 1998;Hagemann and Dümenil, 1998;Ducharne et al., 2003;Ngo-Duc et al., 2007;Paiva et al., 2011;Yamazaki et al., 2011;Decharme et al., 2012).The modelling of the river flow velocity is addressed in several ways in the literature.Coe (1998) and Oki and Sud (1998)

considered a uniform and
Published by Copernicus Publications on behalf of the European Geosciences Union.2136 C. M. Emery et al.: Large-scale hydrological model river storage and discharge correction constant flow velocity over the entire basin in SWAM and TRIP RRMs, respectively.Other studies rather use a spatially distributed, but still constant in time, flow velocity based on topography and river channel characteristics (Vörösmarty et al., 1989;Hagemann and Dümenil, 1998;Ducharne et al., 2003).However, most recent models rely on a time-varying and spatially distributed flow velocity estimation based on the Manning formula (Manning, 1891), e.g.Arora et al. (1999), Ngo-Duc et al. (2007), Lucas-Picher et al. (2010) and Decharme et al. (2012).The first RRMs only modelled water flowing in the river channel (Vörösmarty et al., 1989;Coe, 1998;Oki and Sud, 1998).Subsequent RRM developments included the modelling of the groundwater inflow to the river (Hagemann and Dümenil, 1998;Arora et al., 1999;Ducharne et al., 2003;Decharme et al., 2008;Lucas-Picher et al., 2010) as well as the floodplains-river dynamics (Decharme et al., 2008;Paiva et al., 2011;Yamazaki et al., 2011).The routing network is derived from a digital elevation model (DEM): some of them remain defined on a regular mesh grid (Oki and Sud, 1998;Decharme et al., 2012), while others use an irregular discretization by sub-catchments, such as MGB-IPH (Paiva et al., 2011) and CaMa Flood (Yamazaki et al., 2011), or by river reaches, such as RAPID (David et al., 2011).More information on some global LSMs and their associated RRMs (with the corresponding references) could be found, for example, in Schellekens et al. (2017).
However, even if hydrological models become more and more accurate, inherent model uncertainties are unavoidable.They originate from several sources: simplification and lack of knowledge in the real physics, numerization and discretization-induced errors and uncertainties in the input parameters and forcing.All these uncertainties impact a model's outputs.In the worst case, all those uncertainties could accumulate and result in the collapse of the model.The model gives therefore an approximate view of the system's real state.
Observations of the system can be used to calibrate and/or validate the model and reduce its errors.These observations can be obtained from in situ or remote techniques.In situ techniques mainly focus on measuring river water elevations at a gauge station.Another important variable of interest in river hydrology is the river discharge, which is sparsely measured compared to water elevation.Based on river discharges and elevations measured at the same time and at the same location, it is possible to build a rating curve that represents the elevation-discharge relationship.This rating curve is then applied to water elevation to set continuous discharge time series.Institutions delivering in situ data provide mainly discharge.Even though in situ measures are generally quite accurate with a high time sampling (i.e.sub-daily), their main limitation is their local and spatially sparse sampling over the river network.Furthermore, nowadays, remotely sensed data from satellite missions are more and more available and provide useful observations of rivers.The most straightforward and used instrument to measure river water elevations is the nadir altimeter.
DA aims to improve model skills to forecast/simulate the physical system evolution.To do so, DA techniques focus on either correcting the model's input parameters (parameter estimation) or the model's outputs (state estimation).State estimation (SE) consists in using observations to directly correct the model output state.It is based on the assumption that the model (and the observations) are known to be imperfect.So, SE aims at correcting model outputs, whose errors result from all sources of uncertainties previously described.
SE has been widely used in oceanography and meteorology (Evensen and Leeuwen, 1996;Houtekamer and Mitchell, 2001;Gustafsson et al., 2012;Tanajura et al., 2015).However, DA of remotely sensed observations to correct hydrological model states is more recent.Moreover, it is more developed for LSMs than RRMs, as shown for example by the global-scale Land Data Assimilation System of the NASA Goddard Earth Observing System (Reichle et al., 2014).This platform assimilates simultaneously the SMOS soil moisture product, MODIS snow cover extent fraction and integrated GRACE terrestrial water storage variations into an ensemble Kalman filter (EnKF) to correct the states of several LSMs.Other studies assimilate similar kinds of observations, along with in situ data, into smaller-scale hydrological models (Trudel et al., 2014;López López et al., 2016).As for RRMs, to the authors' knowledge, there are few studies where remotely sensed and/or in situ data are assimilated into global-scale RRMs.However, in the liter- ature, there are several studies that used assimilation techniques at smaller and local scales with finer spatial resolution than global RRMs, using mostly in situ data (Schumann and Domeneghetti, 2016).For example, Clark et al. (2008) applied the EnKF to correct soil, aquifer and surface water storage in a small river in New Zealand (the Wairu).More particularly, they used gauge discharge data from four gauges to correct water storages in the 380 sub-catchments dividing the study zone.Paiva et al. (2013a) also used an EnKF over the Amazon basin for three different experiments, which assimilate, first, ENVISAT water surface anomalies from 212 virtual stations, and then discharge data from 109 gauges and finally remotely sensed discharges from 287 stations obtained from Getirana and Peters-Lidard (2013).This study aimed at correcting discharge estimated by the MGB-IPH hydrological model over more than 5000 sub-catchments comprising the Amazon River basin.Moreover, in two different studies, Michailovsky et al. (2013) and Michailovsky and Bauer-Gottwein (2014) assimilated for the Brahmaputra River in Asia and the Zambezi River in Africa, respectively, using an extended Kalman filter, ENVISAT water surface elevation measurements from six and nine virtual stations in Michailovsky et al. (2013) and Michailovsky and Bauer-Gottwein (2014), respectively, to correct simulated water volumes in 18 and 37 sub-catchments, respectively.
The objective of the present study is to investigate the contribution of remotely sensed data, and in particular measurements derived from nadir altimeters that provide local information, to improve a large-scale RRM via DA.The scale difference between the observations and the model leads us to also study the need to use localization methods within our DA framework.We used an ensemble Kalman filter, to which we added a simple localization module, to assimilate discharges derived from ENVISAT water surface ele-vation measurements.These observations are used to correct the state of the large-scale Total Runoff Integrated Pathways (TRIP, Oki and Sud, 1998) RRM version included in the Surfaces Externalisées land surface modelling platform (SurfEx, Masson et al., 2013) and developed at the Centre National de Recherches en Météorologie (CNRM, France).This particular version is denoted by the CTRIP acronym hereinafter.CTRIP is coupled with the Interactions-Soil-Biosphere-Atmosphere (ISBA, Boone et al., 1999) LSM at a resolution of 0.5 • × 0.5 • .In Sect.2, we present the study domain along with the ISBA-CTRIP model version and remotely sensed product used in this study.Section 3 provides first a general presentation of the ensemble Kalman filter (EnKF) DA method.Then we introduce the special features associated with the study and the description of the assimilation strategy.Then, in Sect.4, we present results for a series of DA experiments testing the ensemble generation strategy and the correction of different state variables.Finally, Sect. 5 discusses these results and some perspectives.The last section gives the conclusions and some perspectives of the study.
2 Study domain, model and data used

Study domain: the Amazon basin
The study is focused on the Amazon River basin (see Fig. 1a).It is the world's largest river in terms of averaged discharge (2 × 10 5 m 3 s −1 ) and drainage area (6.15 × 10 6 km 2 ).The discharge at its mouth represents 30 % of total freshwater inflow to the Atlantic Ocean (Wisser et al., 2010) and its catchment area covers about 40 % of South America.The river source is located in the Peruvian Andes and flows through the Brazilian rainforest while receiving water from several important tributaries: first, the Ucayali, the Japurá River, the Purus River and, at Manaus, the Negro River (14 % of the total discharge).At this point, the river has reached 56 % of its total discharge.From Manaus to its mouth, it receives water from the Madeira River (17 % of the total discharge), the Tapajós River and the Xingu River (11 % of the total discharge all together) (Molinier et al., 1993).
The Amazon basin's geology can be divided into three major morpho-structural units: the western Andean Cordillera, the central Amazon trough and the shields at the eastern part of the basin (Guiana shield to the north and the Brazilian shield to the south).The northern and southern regions of the basin are under a tropical climate with a dry and a wet season, but the maximum rainfall season for the two parts occurs at different periods during the year (Meade et al., 1991).This implies that annual peak discharge in southern tributaries occurs a few months earlier than in northern tributaries.Meanwhile, the central basin is under an equatorial climate zone, implying high surface temperatures, air humidity and, especially, precipitation.Thus, a vast floodplain along the mainstream is filled every year, leading to the damping of discharge extremes.

Model presentation
The ISBA model (Noilhan and Planton, 1989) is a relatively standard land surface model (LSM) defined over a regular mesh grid at global scale.The model's equations are solved for each grid cell separately from the others.All grid cells are only correlated through the spatial patterns of atmospheric (especially precipitation) and radiative inputs, vegetation cover and soil composition.By taking into account the heterogeneity in precipitation, topography and vegetation within each grid cell (Decharme and Douville, 2006) and based on the force-restore method (Blackadar, 1976), ISBA gives a diagnosis of the water and energy budgets in each grid cell.Especially the ISBA-3L configuration (Boone et al., 1999), used in Alkama et al. (2010) and Decharme et al. (2012), has been chosen for the present study.In this version, the soil is divided into three layers: the superficial layer, the root zone and the sub-root zone.Precipitation can either fall directly on the soil surface or be intercepted by the canopy.The soil water content varies with canopy dripping, surface infiltration, soil evaporation, plant evapotranspiration, surface runoff and deep drainage (for more details, see Alkama et al., 2010;Decharme et al., 2012).Then, ISBA gives a diagnostic of each water budget component, in particular the surface runoff (Q ISBA,sur ) and gravitational drainage (Q ISBA,sub ) which are the main inputs for CTRIP.
The CTRIP RRM is also defined over a regular mesh grid.In this study, it is run at the same resolution as ISBA (0.5 • × 0.5 • ).CTRIP is dedicated to the lateral transfer of water from one cell to the other, up to the continent-ocean interface following a river network (Oki and Sud, 1998).The CTRIP version (Decharme et al., 2010) used in this study is coupled with the ISBA LSM and was subsequently developed by Decharme et al. (2010Decharme et al. ( , 2012)).It consists of a system of three reservoirs (see Fig. 1b): the surface reservoir S (kg) modelling the river, the groundwater reservoir G (kg) and the floodplain reservoir F (kg).
Only the surface reservoir S sends water from cell to cell based on the TRIP routing network.A cell can receive water from several upstream cells, but sends water into a unique downstream cell based on a space and time-varying flow velocity v(t) estimated with the Manning formula (Manning, 1891).For any given cell, TRIP inputs are the TRIP outflow from upstream cells Q S in,TRIP (t) and the ISBA surface runoff for that cell Q ISBA,sur .Moreover, S receives water from the groundwater reservoir G, Q G out (t), and can exchange water mass with the floodplain . G receives inflows from ISBA gravitational drainage Q ISBA,sub and outflows to the river reservoir S.This outflow represents more a delayed contribution of the gravitational drainage to the river than a real groundwater dynamic.
The floodplain scheme activates when the water height in the river, h S , exceeds a given critical bankful height H c .Then, part of the precipitation is intercepted by the floodplain (P F (t)) and the water in the floodplain can either evaporate (E F (t)) or infiltrate into the soil (I F (t)).A detailed description of the floodplain scheme is given in Decharme et al. (2008Decharme et al. ( , 2010Decharme et al. ( , 2012)).

Model implementation over the Amazon
ISBA-CTRIP is run in offline mode.This implies that external atmospheric data are needed to force the model.Here, the atmospheric data from the Global Soil Wetness Projet 3 (GSWP3, http://hydro.iis.u-tokyo.ac.jp/GSWP3) are used.The project consists of three global-scale experiments with the objective of investigating long-term changes in the energy-water-carbon cycle components and their interactions.The 3-hourly resolution atmospheric boundary conditions used in the present study were generated by dynamically downscaling the global 2 • -resolution 20th Century Reanalysis (Compo et al., 2011).This reanalysis assimilates several atmospheric observations into the Climate Forecast System (CFS) operational model from NCEP (National Centers for Environmental Prediction).
For ISBA-CTRIP, the Amazon basin is composed of a total number of 2028 cells.A sensitivity analysis (SA) of the ISBA-CTRIP has been conducted by Emery et al. (2016).In this analysis, the basin was divided into nine hydrogeomorphological zones which are shown in Fig. 2.These zones were designed to take into account different components: (1) a hydrological component (the main course is separated from the tributaries which have their own zones); and (2) a geological component (the three major morpho- structural units are distinguishable).The nine zones are the following: (1) the upstream Andean part of the basin until the city of Iquitos, Peru; (2) the mainstream from Iquitos to Óbidos; (3) the mainstream from Óbidos to the river mouth; (4) left-bank tributaries from the Napo River to the Japurá River; (5) left-bank tributaries from the Japurá River to Óbidos, including the Negro River and its drainage area; (6) right-bank tributaries from Iquitos to the Purus River confluence at Anamã; (7) right-bank tributaries from Anamã to Óbidos, including the Madeira River; (8) right-bank tributaries exiting in zone 3, including the Tapajós River and the Xingu River; and (9) left-bank tributaries exiting in zone 3.This subdivision will be used within the DA platform.

Altimetry-based discharge product
The altimetry-based discharge product used in this study is derived from water surface elevations measured by the EN-VISAT Radar Altimeter-2 altimeter instrument at Virtual Station (VS).VS is computed where the altimeter track crosses the river.The ENVISAT mission operated from September 2002 to October 2010 on its nominal orbit, which has a 35-day repeat period and an 80 km inter-track distance at the Equator.The water surface elevations measured over the Amazon basin were initially generated by Silva et al. (2010).The final product was referenced to the EGM2008 geoid (Palvis et al., 2012) and the vertical precision ranged from 12 to 30-40 cm for most of the stations (and can reach several metres for the worst stations).
Turning water surface elevation measures into an equivalent discharge requires the use of elevation-discharge rating curves.The rating curves used in this study have been built and validated by Paris et al. (2016), using water surface elevations from ENVISAT (Silva et al., 2010(Silva et al., , 2012(Silva et al., , 2014) ) and discharges simulated by hydrological-hydrodynamic model MGB-IPH (Model de Grandes Bacias-Instituto de Pesquisas Hidráulicas, Collischon et al., 2007).The model's original version, developed over the Amazon River basin, consists of a large-scale distributed hydrological model coupled with a hydrodynamic module that uses a simple storage scheme for floodplains (Paiva et al., 2011).The entire basin is divided into 5765 elementary catchments with an area varying between 100 and 5000 km 2 .A surface scheme is applied for each mini-basin to estimate the main flows and a routing network is used to direct the flows from one elementary catchment to another, down to the outlet.Two approaches are used to estimate the river discharge: 1 -the Muskingum-Cunge method (MC) for basin heads and small tributaries; 2 -the Saint-Venant equations (HD for hydrodynamic) for the main stem and main tributaries.The DEM used is SRTM (Farr et al., 2007), and parameters such as the river width and depth are determined using geomorphological relationships calibrated over the sub-basins (Paiva et al., 2013a).Moreover, the model version used to determine the rating curves is the version developed by Paiva et al. (2013b), where gauge discharges are assimilated into the model via an ensemble Kalman filter (EnKF, Evensen, 2003), over the period between 1998 and 2010.The assimilated discharges allow one to correct the simulated discharges over both gauged and ungauged elementary catchments.With a better estimation of discharges, Paiva et al. (2013b) also provide an estimation of discharge uncertainty (modelled as a white noise) for each elementary catchment.
The MGB-IPH discharges were used by Paris et al. (2016) as a baseline to estimate the altimetric rating curves such that where -H alti,i is the altimetric water surface elevation at the ith virtual station which is available at time t alti,i , ) is the discharge estimated by the MGB-IPH model, at time t alti,i , in the j th mini-basin which contains the ith virtual station, and a i , z 0,i and b i are the rating curve parameters to be determined.Those parameters are constant in time but vary from one virtual station to the other.
To calculate those parameters at each virtual station, a global optimization algorithm, the Shuffled Complex Evolution Metropolis developed by Vrugt et al. (2003), was used.It allowed determination of rating curves for 767 ENVISAT virtual stations.More details about the rating curve computation can be found in Paris et al. (2016).Once rating curve parameters are determined, altimetric water surface elevations are easily converted into equivalent "altimetric discharges".
Moreover, the altimetric discharges are provided with an estimation of their uncertainty including the normalized deviation from the MGB discharge.
The quality assurance of the discharge product has been made by constraining the rating curve coefficients within a physical range of values (Paris et al., 2016).Paris et al. (2016) also conducted a sensitivity analysis that showed a small sensitivity of the coefficient estimation to their first guess value.The quality check was done by comparing the satellite-derived discharge to the modelled discharge over a validation time period distinct from the calibration period used to derive rating curves.Discharge was also compared to some in situ gauges.Satellite-derived discharge is of course heavily correlated with the model accuracy.Overall, a comparison to 51 gauge measurements led to a mean Nash-Sutcliffe coefficient of around 0.8 and a normalized root mean square error of around 10 % over the validation period (Table 8 in Paris et al., 2016).
Altimetric discharges have then to be compared to ISBA-CTRIP discharges.However, while the virtual stations are irregularly distributed over the entire basin, the model is defined over a coarse regular mesh grid of 0.5 • × 0.5 • .A preliminary treatment of the virtual stations is applied to associate each ENVISAT virtual station with an ISBA-CTRIP cell with respect to their localization and the drainage network.The following algorithm has been used.
-The CTRIP river network is compared to a realistic river system (produced with GoogleEarth) to properly associate ISBA-CTRIP cells with a given tributary in the basin.
-Then, each virtual station is coupled with the closest ISBA-CTRIP cell along the same tributary.It may be the cell containing the virtual station or an adjacent cell according to the river network.
This algorithm allowed association of most of the virtual stations with a unique CTRIP cell.However, some particular cases have been treated.First, some virtual stations were located on tributaries too small to be represented on the CTRIP river network.In this case, the virtual station was not included in the study.Then, there were several very close virtual stations associated with the same ISBA-CTRIP cell.In this second case, the virtual stations with the lowest deviation from the MGB discharge were conserved.Finally, over the 767 ENVISAT virtual stations initially available, 368 ENVISAT virtual stations were kept and associated with an ISBA-CTRIP cell.Among them, 19 % (69 virtual stations out of 368) have been associated with an adjacent cell.

In situ discharge product
At a national or basin scale, water agencies can share discharge time series, such as the Agencia Nacional de Agua (ANA, hidroweb.ana.gov.br) in Brazil for the Amazon River basin.For the present study, we retrieved discharge time series from 145 ANA in situ stations over the entire basin.These gauge discharges have then been used to evaluate the performances of the DA (but they have not been assimilated into ISBA-CTRIP).

Method
The purpose of the SE DA is to correct model outputs using observations while taking into consideration uncertainties in both the model and the observations.In this work, as observed data correspond to discharge estimates, we chose to correct model output variables such as discharge or river storage.Indeed, following the results from the ISBA-CTRIP sensitivity analysis (SA; Emery et al., 2016), discharge is mainly sensitive to water inflow.Figure 3 presents the general DA method in the present study.The figure reads from top to bottom and from left to right.Three types of state variables will be considered: the river initial storage, the river final storage (which are both the main ISBA-CTRIP state variables) and the river discharge that will all be compared to the observed discharge.All three can be corrected through assimilation with specific treatment that will be detailed in the following sections.The DA will use several operators (in Fig. 3, Z k and H k ) that link state variables with each other.

Data assimilation variables
The DA technique implemented in the present study is a sequential EnKF (Evensen, 2003).Here we shortly give the mathematical formalism used in the rest of the paper and a brief description of the EnKF method.First of all, the the term "assimilation window" used hereafter corresponds to the period during which a complete assimilation cycle is conducted.It is delineated by two consecutive observation times and will be denoted by [k − 1, k].From now on, the kth assimilation cycle will be the cycle starting at time k − 1 and assimilating the available observation(s) at time k.

Control variables
The vector x k ∈ R n x is called the control vector.It includes the n x uncertain variables to be estimated during the kth DA cycle (within the time interval [k − 1, k]).As stated before, control variables are prognostic or diagnostic variables of the ISBA-CTRIP model.Prognostic variables are the physical unknown in the differential equations' system that describes the model's behaviour.Diagnostic variables are also physical variables, but they are estimated from the prognostic variables.The choice of the control variables determines the observation operator H k that maps the control variables into the observation space: where y k are the control variables equivalent in the observation space, also called model observations.They are then compared to the measured observations y o k (described in Sect.3.1.2)during the DA step.
Unlike hydrodynamic models, which directly solve Saint-Venant equations and for which discharge is a model state variable (or prognostic variable), the hydrological model ISBA-CTRIP solve differential equations describing the time evolution of water stock in the river (S), the groundwater (G) and the floodplain (F ).Then, water elevation and river discharge are diagnostic variables derived from these prognostic variables.In CTRIP, river discharge Q S out is computed as follows: with L (m) river section length, v the flow velocity (estimated from the Manning formula) and S the surface water mass.Therefore, three types of variables can be considered as control variables in the data assimilation scheme: the discharge Q S out (denoted Q in the remaining of the study to simplify notation, which is a diagnostic variable), the river final water stock S end (a prognostic variable) or the river initial water stock 5also a prognostic variable).Definition and complexity of the observation operator H k , that maps the control space into the observation space, depends on the nature of the control variable.These three options are presented below.
-Option 1: correcting ISBA-CTRIP discharges: for this option, the control variables, gathered into the vector x k , are the ISBA-CTRIP discharges Q i,k , i = 1 . . .n x = 2028 (number of TRIP cells in the Amazon basin) estimated for all 2028 cells in the TRIP Amazon basin, at the assimilation cycle k.
The observation operator H k resumes to a selection operator S k which selects the observed TRIP cells at the current assimilation cycle: This operator is linear.The difficulty with this operator is that, once the assimilation analysis is produced, it is necessary to convert the analysis discharge Q a i,k , i = 1 . . .n x (i.e. the corrected discharge obtained after assimilation) into the equivalent final water stock S a end,i,k .Indeed, as already stated, in ISBA-CTRIP, discharge is not a prognostic variable.Correction from the assimilation step needs to be propagated to the model prognostic variables, here the river final stock.Moreover, the analysis of the final water stock S a end,i,k will be used as an initial condition for the model run until the next assimilation cycle: ∀i = 1 . . .n x , ∀k, S a end,i,k = S b ini,i,k+1 .However, the ex-act relationship linking discharge to the final river stock is unknown.
A possible solution consists in inverting Eq. ( 3).Assuming that the discharge estimated by ISBA-CTRIP Q a i,k is the instantaneous flow at the final time of the integration window, We obtain that (for more details on this approximation, see Appendix A) with ρ (m 3 kg −1 ) the water density, L (m) the river section length, W (m) the river width, s (-) the riverbed slope and n (-) the Manning coefficient in the riverbed.
Then, for experiments with discharges as control variables, the formula in Eq. ( 5) will be used to convert corrected discharges into river stock and then propagate the model to the next observation time.
-Option 2: correcting ISBA-CTRIP final water stock: for this option, the control variables, gathered into the vector x k , are the ISBA-CTRIP final water stock S end,i,k , i = 1 . . .n x estimated for all 2028 cells in the TRIP Amazon basin, at the assimilation cycle k.
The computational cost for this option is the same as for the first option but, now, the observation operator is defined as where Z k is the operator (implicitly defined within TRIP) that turns the river final stock S end,i,k into equivalent discharge Q i,k .Even though H k is not linear any more, this option presents the advantage of correcting the river final stock S a k,end that can be directly used for the next assimilation cycle and no additional uncertainties are introduced.However, the corresponding analysis discharge Q a i,k is now unknown as the explicit expression of Z k is also unknown.A potential formula to determine Q a i,k can be deduced from Eq. ( 5).Such a formula will be necessary to make comparative statistics to ENVISAT and in situ discharges and be able to evaluate the assimilation performances.
-Option 3: correcting ISBA-CTRIP initial water stock: for this final option, the control variables, gathered into the vector x k , are the ISBA-CTRIP initial water stock S ini,i,k , i = 1 . . .n x estimated for all 2028 cells in the TRIP Amazon basin, at the assimilation cycle k.
The discharge observations are used to correct the surface water stock at the time prior to the observation time or, in other words, at the initial time of the integrating window.Therefore, the observation operator is written as the composition of the model operator M [k−1,k] with the observation operator defined in Eq. ( 6): This operator is highly non-linear as it contains the full model operator.However, it is the only option where no uncertainties are added from the use of an external formula to compute corrected discharge at the observation time.Uncertainties in the stock-discharge relationship are only due to the model uncertainties.Indeed, once the analysis initial water stock S a ini,k,i is determined, the control variable update must be propagated through the next assimilation time by re-integrating the ISBA-CTRIP model over the assimilation window.The initial water stock S a end,k,i and the analysis discharge Q a k,i are automatically determined during this run.

Observation variables
In the framework of the state estimation, the observation variables, at a given day within the Amazon basin, are the discharge estimates derived from ENVISAT water surface elevations at the virtual stations associated with an ISBA-CTRIP cell.The ENVISAT repeatability is 35 days, and therefore a given virtual station will provide an observation every 35 days at best.During the data assimilation experiments, all virtual stations will be used simultaneously.Because of the ENVISAT orbit, the number of available observations at a given day will vary between 0 and 15, and these observations will be assimilated daily via the EnKF.Then, the observation vector y o k at the assimilation cycle k (equivalently, at the simulation day k) is composed of the n y,k discharge measures available at day k: where y o k,j , j = 1 . . .n y,k is the j th observation among the n y,k at cycle k.
Measurement errors m k come from errors in ENVISAT water surface elevations, errors in MGB discharges and uncertainties in the rating curve parameters used to turn water surface elevation into discharge.Sorooshian and Dracup (1980), Clark et al. (2008) and Paris et al. (2016) noticed that the concavity of the elevation-discharge relationship implies that the higher a water elevation, the more uncertain the corresponding discharge.Therefore, the observation error standard deviation σ o k,j , associated with the j th observation at cycle k, is defined with respect to the instantaneous discharge measure y o k,j , i.e.: where η o j ∈ [0, 1] is a constant depending on the virtual station, such that η o j models the relative error.[%] and determined by Paris et al. (2016).As MGB discharges were used to determine ENVISAT discharges from ENVISAT water elevations, σ mgb j represents the deviation between the two discharges data.Second, to take into account that MGB discharge is not perfect (in other words, to take into account some deviation from the real discharge), 0.05 is added to σ mgb j and the sum is rounding up to the nearest whole number, giving where the function E is the ceiling function.Finally, η o j is equal to the first multiple of 5, above η mgb j .At the end, η o j varies from 0.10 to 0.35 over the entire basin and is constant in time.Besides, the representativeness error r k induced when a virtual station is associated with cells of the coarse TRIP mesh grid is neglected here.
Moreover, for a given assimilation cycle and also between different cycles, the observations are considered uncorrelated in space and time.The observation error covariance matrix at cycle k R k is then a diagonal definite positive square matrix.

The ensemble Kalman filter (EnKF) for
ISBA-CTRIP state estimation

The EnKF sequential estimation
In the EnKF framework, the model and observation operator are not linear.Therefore, the main idea is to use stochastic ensembles to represent the control variables PDFs along with the error models (Evensen, 1994(Evensen, , 2003)).First, the background control variables x b k are stochastically represented by an ensemble of n e members: Each member is estimated separately through the EnKF prediction step.For each control variable case (see Sect. 3.1.1),each member of the control ensemble X b e,k is estimated by integrating the model operator from the corresponding analysis member at the previous assimilation cycle, while adding external uncertainties (see Sect. 3.2.3): Then, the background control ensemble must be compared to the observations.Depending on the control variables nature, the model operator is already included (option 3) or not (option 1 and 2) within the observation operator.Besides, following Burgers et al. (1998), an additional noise o k is added to the observation vector y o k to avoid ensemble collapse.The observation ensemble thus obtained is noted: Finally, the EnKF analysis step is applied to each member of the ensemble such that where the direct non-linear observation operator is applied to convert the control variables into equivalent model observations.
The particularity of the EnKF is that the Kalman gain (K e,k ) is stochastically estimated from the different control and model observation ensemble, as follows:

Localization of the error covariance matrices
In the framework of state estimation, the sampling error can introduce artificial correlations into the background/analysis error covariance matrices, and generate spurious correlations between two distant grid cells in the mesh (Anderson, 2007).
The ensemble size n e is limited by computational constraints.Therefore, before the EnKF analysis step, a numerical processing of the matrices [PH T ] e,k and [HPH T ] e,k is necessary to suppress spurious correlations that can potentially degrade the analysis.Localization methods are designed to reduce these problems.
There exist two types of localization techniques (Greybush et al., 2011;Sakov and Bertino, 2011).The first one is called B-localization.It is based on explicitly modifying the background error covariance matrix P b e,k .It consists in multiplying the matrix P b e,k by a correlation matrix generated from a radial function, namely a function of the two/three spatial dimensions which monotonously decrease with the distance between control variables (Hamill et al., 2001;Houtekamer andMitchell, 2001, 2005).This modified matrix replaces P b e,k in the calculation of the Kalman gain matrix K e,k .The other common localization technique is called R-localization or local analysis.This one consists in performing the analysis step in characteristic sub-spaces of the overall problem space.
However, all these localization techniques described above have been developed for atmospheric modelling where problems are in two or three dimensions.The use of localization in hydrology is more limited.Several studies exist to improve subsurface flow modelling (Devegowda et al., 2010;Delijani et al., 2014), but these approaches have a dimensionality close to atmospheric approaches, as they take place in a continuous medium in two or three dimensions.Other studies using localization allow estimation of better model parameters, still continuously defined in two or three dimensions (Sun et al., 2009;Rasmussen et al., 2015).
The localization method used with the CTRIP river routing model is of the B-localization type.However, it can not be simply defined on a two-dimensional radial function.Indeed, the river flow is along several one-dimensional flow directions, modelled by the routing network.The localization technique must consider the routing network to decorrelate adjacent cells on the mesh grid but located in two different sub-catchments.Nevertheless, along a same flow direction, the correlation between two distinct cells depends on the distance between the two cells.Then, for each assimilation cycle, the localization consists in a localization mask delimiting an influence area for each observation.These influence areas gather a limited number of neighbouring downstream and upstream cells around the observed cell with respect to the river routing network.We chose a fixed localization scale for simplicity and as a first step in the feasibility study of the development of a localization method for a hydrology application.
To determine the number of cells defining the influence area, the basin subdivision into nine hydrogeomorphological zones is used with a mean flow velocity for each zone.The influence area, for a given observed cell, is given by the criteria below.For an influence area of size p cells, the area is composed of the following.
-The observed cell.
-The p downstream cells according to the river routing network.
all the cells upstream the observed one covering p upstream levels.However, the going up stops when the hydro-geomorphological zone is different from the one of the observed cell.
The number of cells within the influence area depends on the mean flow velocities (averaged over a year of simulation) in the zone in which the considered cell is situated.Those mean velocities are calculated from the free run simulation, namely the ISBA-CTRIP simulation realized without any assimilation step.The ISBA-CTRIP resolution is 0.5 • × 0.5 • , or approximately 50 km × 50 km.Given the river meanderings, the minimal covered distance through a cell is of 50 km.Furthermore, by comparing the free run discharge to in situ and ENVISAT discharges, it seems that the free run underestimate discharge (and so the flow velocity).Consequently, to fix the number p of cells defining the influence area in each hydro-geomorphological zone, the following steps have been performed: 1. the mean velocity for the cells into a given zone is converted into an equivalent distance in km, 2. the maximal distance within the zone is kept and rounded to the closest higher multiple of 50, 3. the number p determining the size of the influence area is the number of cells covered by the maximal rounded distance, knowing that 50 km = 1 cell.
The number of cells into the influence area is presented, for each zone, in the Table 2.
The final localization mask is presented into a matrix of size n x × n y,k (with n x , the number of control variables, and n y,k , the number of observation variables, at the assimilation cycle k) containing only 0 and 1.The localization mask S has the same dimension as the matrix [PH T ] e,k .So, the localization matrix is term-to-term multiplied (sign "×" in Eq. 15) to the error covariance matrix [PH T ] e,k such as in Moore (1973) and Biancamaria et al. (2011): To then extend the localization to the error covariance matrix [HPH T ] e,k , the lines in [PH T ] e,k corresponding to the observed cells are extracted to form the second localization matrix.This second matrix is also term-to-term multiplied to [HPH T ] e,k .This localization step is applied in each assimilation cycle with respect to the activated ENVISAT virtual stations at the current assimilation cycle.

Generating the ensembles
The background error covariance matrices ] e,k are estimated from the control variable ensemble using the definition suggested by Evensen (2004), Moradkhani et al. (2005) and Durand et al. (2008).To get a large ensemble, while maintaining a reasonable computational time, the ensemble size n e has been set to a hundred members.Details on how they are exactly calculated are given in Appendix B. These matrices have a n x × n y,k and n y,k × n y,k size, respectively.The elements in the error covariance matrices, depend only on the definition of the background ensemble stored in the control matrix X b e,k and the parameter uncertainties taken into consideration to generate H(X b e,k ).In the framework of state estimation, we choose to consider uncertainty into the initial condition and uncertainty into the precipitation forcing (Nijssen and Lettenmaier, 2004;Andreadis and Lettenmaier, 2006;Clark et al., 2008;Paiva et al., 2013b).
-Perturbation of the initial condition: the vector containing initial surface reservoir storage for each n x = 2028 CTRIP cell at the assimilation cycle k is called c k .To ease the notations, we will omit, as much as possible, the assimilation cycle k subscript, knowing that, for all randomly perturbed variables and constants, a new ensemble is generated at each cycle.
We used the Amazon basin division into n s = 9 hydrogeomorphological zones (see Fig. 2).Initial conditions are perturbed by applying a multiplying factor over each zone η where c s,k is the reduction of the initial condition c k to the only zone s.For the perturbation to vary from one member to another, the value η c, [l]   s is the realization of a Gaussian distribution, different for each member [l] = 1 . . .n e and for each hydrogeomorphological zone s.The Gaussian distributions used have a mean value of 1 and a standard deviation of σ η c s,k whose values are detailed in Table 1.
The η c, [l]   i values depend on the assimilation cycle k and on the hydrogeomorphological zone in which the ith cell is.
-Firstly, a more important perturbation is applied to cells situated on the river mainstream (zones 2 and 3), as we assume that the uncertainties are more important in those zones.Indeed, discharges in these zones are the highest of the entire basin.
Besides, several cells are confluence cells and are subject to backwater effects.As ISBA-CTRIP does not model the backwater effects, the water stock uncertainties in these cells are increased.
-Secondly, at the first assimilation cycle, the initial condition before perturbation c 1 is identical for every member.At this particular cycle, the ensemble variance after perturbation depends only on the perturbation method presented in Eq. ( 16) while, for the other assimilation cycles, the successive previous assimilation cycles introduced an additional variability between members, before the perturbation step in Eq. ( 16).Therefore, the initial condition variance is more important at the second assimilation cycle and after.Then, to generate a larger variability at the first assimilation cycle, the standard deviation σ η c s,k of the variable η c, [l]   i is more important at the first cycle than for the others.
-Perturbation of precipitations: another source of uncertainties for the generation of the ensemble H k (X b e ) lies in the precipitation fields.Atmospheric forcing comes from the GSWP3 product (see Sect. 2.2.2).Precipitation forcing F has been perturbed using the procedure presented by Clark et al. (2008).The ensemble of perturbed precipitation fields F e is defined such that where -F is the two-dimensional field of precipitation forcing before perturbation (with a time step of 3 h), -F [l] , for l = 1 . . .n e , is the lth perturbed precipitation field, ϕ [l]  p , for l = 1 . . .n e , is the lth multiplying uniformly distributed field of F to generate F [l] .More details on how the fields ϕ [l]  p have been generated are given in Appendix C.

Assimilation diagnostics
During the assimilation experiment, it is necessary to quantify the assimilation performances.The quality of the assimilation will be evaluated in a given cell i by estimating the root mean square error (RMSE) between the simulated discharge Q * i and the observed discharge K † represents the total number of available observed discharges Q † i at the studied cell for the study period.The " * " superscript can be either the superscript "f" for the free run Table 3. Presentation of the different state estimation experiments.The "SE" acronym stands for "state estimation", indexes "1", "2" or "3" are to differentiate the control variables ("1": initial river storage, "2": final river storage and "3": discharge) and the suffixes "direct", "diag" and "local" indicate the localization scheme ("direct": without localization, "diag": diagonal error covariance matrices and "local": with localization).(without assimilation) or the superscript "a" for the analysis run (after assimilation).The " †" superscript can be either "o" for the observed ENVISAT discharge or "situ" for the gauge discharge.

Exp. name Control variable
Based on this definition, the assimilation performance will be estimated at each cell with the normalized RMSE (RM-SEn) defined by where Q † i,• corresponds to the mean of Q † i,k averaged over the available time steps k.
Also, to evaluate the global performance of the assimilation over the entire basin, a global RMSEn (RMSEn global ) will be determined by RMSEn * , † global = 100 × 1 where I † is the total number of stations and w i a weighting constant depending on the maximal discharge at the ith station (max k (Q † i,.)) and the maximal discharge over the basin (max i,k (Q † .,. )) such that With this weighting, cells along the mainstream and the largest tributaries (with the highest discharges) have more importance in the global statistics than cells located in basin heads.
Besides, the analysis run is available as an ensemble.The statistics will then be estimated for each member of the ensemble and the mean (see Eq. 22) of the ensemble will be presented, such as RMSEn a, † i = where RMSEn a,[l], † i is the normalized root mean square deviation of the lth member of the analysis discharge ensemble.

Assimilation strategy
The state estimation experiments have the objective of testing the different control variables described in Sect.3.1.1.Also, another objective is to test, validate and criticize the localization mask introduced in Sect.3.2.2.In the following, experiment names using the localization will have the "-local" suffix, ones without any localization will have the "-direct" suffix and ones with no correlation between cells will have the "-diag" suffix.The objective of this study is to determine the best strategy for assimilating ENVISAT discharges into the ISBA-CTRIP model using the EnKF to correct the model state variables.The different experiments are presented in Table 3.After analysing these five elementary simulations over a single year, a last experiment will be run over the entire ENVISAT observing period (from September 2002 to June 2010), based on the best configurations.
For all the DA experiments, the observation errors are those described in Sect.3.1.2,and the model errors are those presented in Sect.3.2.3.Moreover, each experiment in Table 3 lasts 365 assimilation cycles of 1 day (so 1 year of assimilation) from 1 January to 31 December 2009.We chose this period as it overlaps with another altimetry mission (namely JASON-2), and future works may include comparison of the two datasets' contribution.The numerous ensemble ISBA-CTRIP simulations were realized with the CALMIP high-performance computation platform (https: //www.calmip.univ-toulouse.fr/spip/)with the EOS supercomputer.
In the SE1-direct experiment, ENVISAT discharges are assimilated to correct the initial surface reservoir storage in TRIP (and inherently TRIP simulated discharges).For this first experiment, a classical EnKF, without any localization treatment of the error covariance matrices

Free run performances
The current section briefly presents the model performance without assimilation called the free run.As all in situ and EN-VISAT VS have been associated with a unique ISBA-CTRIP cell, it is possible to compare observed discharge at these stations to corresponding ISBA-CTRIP simulated discharge.To begin with, a sample of 12 in situ stations, spread over the entire basin (over the mainstream and the main tributaries), is selected.The location and the name of these stations are represented in Fig. 4. Figure 5 compares ISBA-CTRIP free run discharges to in situ and ENVISAT discharges at the 12 stations over 1 year of simulation (year 2009).From this comparison the following observations can be drawn.
-Over the majority of cells where there are both an in situ station and a virtual station, the two discharge time series are similar (but not identical; see Fig. 5, panels 1-3, 5-6, 8-9, and 12).These results are due to the fact that gauge discharges were directly assimilated into the MGB-IPH hydrological model to correct the MGB-IPH estimated discharges (Paiva et al., 2013a).Then, those same estimated discharges were used to calculate parameters of the rating curves between ENVISAT water elevations and MGB discharges (Paris et al., 2016).Even though these rating curves have been derived from a model that assimilated in situ data, ENVISAT-derived discharges depend essentially on the remotely sensed water surface elevation variations (Paris et al., 2016).Therefore, ENVISAT discharges remain independent enough of in situ data.
-Finally, in most cases, the free run discharge is quite different from the observed discharge.At downstream mainstream stations (at Manacapuru and Óbidos in Fig. 5, panels 2 and 3), the ISBA-CTRIP model is not able to reproduce flooding occurring between June and August.Therefore, in the free run, the discharge peak occurs earlier in the year and the discharge variations in this period are faster than the observed discharge variations.Similarly, at most of the right-bank tributary stations, the free run discharge peak is higher than the observed discharge peak (see Fig. 5, panels 7-12).However, the seasonal cycle is well reproduced for all these stations.These results illustrate the necessity of conducting the DA experiments.
Then, Fig. 6 displays the global performances of the free run.For each ENVISAT virtual station (see Fig. 6a) and each in situ station (see Fig. 6b), the RMSEn (defined in Eq. 19) between the simulated and observed discharges is calculated and its value is indicated by a colour at the location of the station over the basin.The results are similar between EN-VISAT and gauge discharges, confirming good concordance between the two discharge datasets.RMSEn shows important deviations in basin heads on most of the tributaries as well as at the confluence between right-bank tributaries and the mainstream.Apart from the confluence and basin heads, the largest tributaries, such as the Negro and the Madeira, are well represented.Concerning global statistics (see Eq. 20), RMSE f,o  global is equal to 71.12 % compared to ENVISAT dis- charges and RMSE f,simu global is equal to 68.96 % compared to gauge discharges.These deviations are likely due to atmospheric forcing, parametrization and modelling errors, especially floodplain parametrization.The DA experiments will focus on correcting the model outputs which result from those uncertainties.

Evaluation of the localization method
The first series of experiments assimilates ENVISAT discharges to correct the ISBA-CTRIP initial river stock (see the three first rows in Table 3).They differ on the definition of the background error covariance matrices [PH T ] e,k and [HPH T ] e,k .The SE1-direct experiment uses the com-plete stochastic matrix defined in Eqs. ( B1) and (B2).In the SE1-diag experiment, these matrices are processed such that covariance between two different CTRIP cells is set to 0 if the two variables are situated in two different CTRIP cells.Lastly, SE1-local is based on the localized version of the matrices presented in Sect.3.2.2.So, Table 4 displays the global RMSEn (see the definition in Eq. 20) for the three experiments compared to the free run global statistics.From Table 4, we can see that the RMSE between the free run discharge and both the ENVISAT and gauge discharges is reduced for all experiments, showing that the data assimilation platform is working correctly.The SE1-diag experiment gives the worst results when compared to both the EN-VISAT discharge and the gauge discharge.Then, compared to ENVISAT discharges, SE1-local gives the best results by reducing the global RMSEn by more than 56 % (49 % for SE1-direct), while SE1-direct presents slightly better global statistics than SE1-local when compared to gauge discharges (RMSEn is reduced by 16.5 % for SE1-direct and by 15.25 % for SE1-local).Overall, the global statistics are more reduced when compared to ENVISAT discharges than to gauge discharges.This is due to the fact that gauge discharges are not directly assimilated, unlike ENVISAT discharges.The next subsections present and analyse in more detail results from each experiment.

SE1-direct results
Figure 7 displays the mean analysis RMSEn (RMSEn a, † i defined in Eq. 22) for each ENVISAT virtual station (Fig. 7a) and for each in situ station (Fig. 7b).First of all, results between the ENVISAT RMSEn and the in situ RMSEn are similar, due to the similarity between ENVISAT and gauge discharge time series at most stations.According to Fig. 7a, the assimilation worked quite well along the main-stream and the main left-bank tributaries, namely the Negro River, the Japurá and the Icá, with several stations where RMSEn is below 20 %.The assimilation performances are more moderate over right-bank tributaries, where RMSEn is mostly between 20 and 60 %.Over the entire basin, RM-SEn remains high in all basin heads, along small tributaries and also at most confluences; see for example the Jutaí-Solimões confluence (RMSEn above 60 %), the Purus-Solimões/Madeira-Solimões/Tapajós-Amazon confluences (RMSEn above 40 %) or the Xingu-Amazon confluence (RMSEn above 80 %).
Figure 8 compares the mean analysis discharge in red line at the 12 stations previously introduced in Sect.4.1.For most stations, we can see that the mean analysis discharges recovers a seasonal cycle closer to the observations than the free run.It is especially true for stations along the mainstream, namely Sao Paulo de Olivenca, Manacapuru and Óbidos (Fig. 8, panels 1-3).Also, for stations along right-bank tributaries (Fig. 8, panels 7-12), the analysis seasonal discharge peak is lowered compared to the free run seasonal discharge peak and fits better the observations.This shows the good functioning of the assimilation platform.Nevertheless, mean analysis discharge for all displayed stations presents a chaotic behaviour with numerous local minima and maxima.We can assume that this behaviour is present for all CTRIP cells in the basin.Moreover, for a given cell, most of these sudden variations are asynchronous with ENVISAT observation dates for this cell.For example, at Serrinha in the left panel in Fig. 9, an ENVISAT observation is available on the 4th day of the 35-day repeat period when big off-peaks appear on the 25th day of the same 35-day repeat period.The right panel in Fig. 9 displays the Serrinha station (red circle) with all ENVISAT observations available during the 25th day (yellow circles).On inspection of the contribution of all these observations to the analysis control variable at Serrinha (not shown here), we find that it is observation number 4 that has the highest impact on the analysis (and not observation number 5, as could be expected).This observation 4, located on a very small Negro tributary, has a low discharge value and is responsible for the low corrected discharge at Serrinha after the assimilation step.
These abrupt variations are completely artificial and directly result from the assimilation processing.Indeed, for days with unrealistic peaks/off-peaks, there are multiple EN-VISAT observations available on the basin, which impact many cells all over the basin, even if they are located on other sub-catchments or tributaries.This is due to the construction of the error covariance matrices [PH T ] e,k and [HPH T ] e,k .As these matrices are generated from the ensemble with a limited number of members, some spurious elements may appear in the matrices and link two cells that are very distant in the basin or even situated on different sub-basin.This first experiment highlights the necessity to treat the error covariance matrices to limit such spurious elements.

SE1-diag results
In the SE1-diag experiment, the error covariance matrices are forced to be diagonal.The objective of such processing on the error covariance matrices is to limit the impact of a given observation only to the observed cell.According to Table 4, the assimilation experiment allowed one to reduce the global ENVISAT, an in situ RMSE, when compared to the free run.However, among all three experiments, it is the one which gives the worst global performances.In this experiment, the chaotic behaviour of the mean analysis discharge is not present any more (not shown here).Nevertheless, the mean analysis discharge remains close to the free run discharge, except for regular peaks/off-peaks at an observation time when it is closer to the observed discharge.Therefore, the information brought by only one local observation is not enough.With the localization (see next section), whose results are presented in the following section, the information of several neighbouring VS is used and should constrain the analysis discharge more.

SE1-local results
SE1-local uses the localization treatment presented in Sect.3.2.2. Figure 10 displays the RMSEn evolution from the SE1-direct to SE1-local experiments for both ENVISAT and gauge discharge.Green colours indicate that the SE1local experiment reduced the RMSEn compared to the SE1direct experiment, while yellow to red colours indicate that the SE1-local experiment increased them.The RMSEn is mostly improved over the entire basin and more particularly along major right-bank tributaries.However, the RMSEn is generally degraded along the mainstream.These maps show the good performances of the localization method over tributaries.Now, it appears that, compared to gauge discharges (see Fig. 10b), the SE1-direct experiment gives better results, especially along the mainstream.Indeed, Table 5 details the local RMSEn a, † i at ENVISAT/in situ stations located along the mainstream and confirms that SE1-direct gives better results.As the global RMSEn is defined as the mean of the RMSEn weighted by the maximum discharge at the station (see Eq. 20), this explains why the SE1-direct experiment gives a better global RMSEn compared to gauge discharge.
Then, Fig. 11 displays the mean analysis discharge for SE1-local compared to the free run discharge, with corresponding ENVISAT discharge and gauge discharge, at the 12 in situ stations already used in Figs. 5 and 8. Except for stations along the mainstream (Fig. 11, panels 1-3) and also at Boa Sorte (Fig. 11,panel 12), the analysis discharge shows less sharp variations.From these results, we can say that the Hydrol.Earth Syst.Sci., 22, 2135Sci., 22, -2162Sci., 22, , 2018 www.hydrol-earth-syst-sci.net/22/2135/2018/ localization scheme is necessary and improves the assimilation.

Discussions on the localization scheme
The localization mask has been built to avoid the effect of spurious correlations between distant cells or ones situated on different sub-basins.The current localization scheme meets this constraint.Indeed, results for the SE1-local experiment are globally improved compared to the previous experiments.
Nevertheless, along the mainstream, the initial experiment without localization gives better results.We can interpret that by the fact that discharge along the mainstream integrates hy-drological processes from all the upstream basins.So when, in the SE1-local experiment, we limit the impact of the observation to only close cells, we suppress part of the information brought by distant cells to mainstream cells.
Therefore, the current localization mask should be improved.The main difficulty here is to determine the size of the influence area for each observation.Currently, this size is predetermined and is constant in time according to averaged flow velocity.A potential development is to consider an influence area size that can vary in time, according to the hydrological season (high-flow/low-flow season).For example, during high-flow season, the flow velocity is higher so is the size of the influence area.Thus, the error covariance matrices would depend on the river time and space dynamic  Table 5. Local RMSEn a,situ i for the SE1-direct and SE1-local experiments at in situ stations along the mainstream (from the most upstream to the most downstream).

Station
RMSEn (as if there were defined from a well-sampled and significant ensemble).

Impact of the chosen control variables
In the second series of experiments, all of them uses the localization scheme (see Sect. 3.2.2) to correct different types of state variables.After assimilating ENVISAT discharge to correct river initial storage (SE1-local experiment), we are testing in a second experiment the assimilation of ENVISAT discharge to correct river final storage (SE2-local) and, in a last experiment, the assimilation of ENVISAT discharge to directly correct river discharge (SE3-direct).These two other experiments need to use an empirical relationship (see Eq. 5) linking simulated river final storage to simulated discharge.
For the SE2-local experiment, the formula is used to convert analysis final rive storage to discharge.Indeed, experiment statistics are based on discharge and, when correcting the final river storage, we do not have an equivalent discharge.
For the SE3-local experiment, the formula is used during the assimilation steps to convert the analysis discharge into an equivalent river storage to propagate in time the corrected discharge.
Table 6 displays the global RMSEn for the three experiments compared to the free run global statistics.For all experiments, the assimilation enables to improve the RMSEn compared to the free run.Also, compared to both ENVISAT and gauge discharge, SE3-local experiment (discharge is the control variable) gives the best results, followed by SE1-local experiment (initial river storage is the control variable).Finally, it is SE2-local experiment (final river storage is the control variable) that gives the worst results, even if it is still improving the RMSEn compared to the free run.
Figures 12-13 display, for each ENVISAT (Figs. 12a  and 13a) and for each in situ stations (Figs.12b and 13b), mean RMSEn difference (in percent) between SE1-local experiment and SE2-local (Fig. 12) and between SE1-local experiment and SE3-local experiment (Fig. 13).Figure 12 shows a slight increase of the RMSEn in SE2-local experiment globally over the Amazon basin except for some basin heads.Also, the upstream part of the mainstream is more degraded (RMSEn increased of more than 60 %).These degraded results imply that the assimilation of discharges may not be adapted to correct the final river stock.However, we need to keep in mind that the analysis discharges are determined from the analysis final river stock using Eq. ( 5).The bad SE2-local experiment results can either be due to bad assimilation results or to an unadapted formula to convert the final river storage into discharge.On the other hand, SE3local experiment gives better general results.As in Table 6 and Fig. 13, the SE3-local experiment shows a global improvement of the RMSEn compared to the SE1-local experiment (apart from a few cells upstream the Amazon mainstream).Indeed, even if Eq. ( 5) is still used to convert the analysis discharges back into river stock, it is used within the assimilation experiment (and not afterwards as for the SE2local experiment).Therefore, the formula uncertainties are accounted for within the EnKF.Also, as the observed discharges are directly used to correct the simulated discharge, it appears logical that the assimilation gives better results as the link between the observed and the simulated variables is immediate.

Discussions
From the different approaches tested in this paper, it appears that there is no one specific configuration that gives the best results for all rivers, when compared to both ENVISAT and gauge discharges.In contrast, the most effective configuration depends on the size and location of the rivers.Along the river mainstream (the Solimões and the Amazon in Fig. 1a), the SE1-direct experiment clearly gives the best results (see the three first rows in Table 7).When the contribution of observations on tributaries is suppressed with the localization, the assimilation is less effective along the mainstream cells (see panels 1 to 3 in Fig. 8 for SE1-direct and compare to the same panels in Fig. 11 for SE1-local).This could be due to the fact that discharge along the mainstream is the result of hydrological processes from the entire drainage area.So, using all available observations helps the EnKF to correct the most efficiently discharge on the mainstream.However, it is different for cells along tributaries.As presented in Table 7, the localization method improves assimilation results for most cells along tributaries compared to the SE1-direct experiment.Along these cells, the localization allows the impact of observation from different sub-basins to be suppressed, especially the ones that are not connected to these cells.Finally, the comparison between the two experiments with localization (i.e.SE1-diag and SE1-local) shows that the local experiment (SE1-local) performs better than SE1-diag.This result was expected, as the localization mask is more realistic in SE1-local, because there is more than one cell impacted by the correction from one observation, in contrast to SE1-diag.
Nevertheless, among all experiments (see Table 3), the one producing the best results globally is SE3-local, where the localization method is used to directly correct the discharge.Therefore, the SE3-local configuration is used for an 8-year experiment, from 25 September 2002 (first date with an ENVISAT observation of the study domain) to 24 September 2010 (last date with an ENVISAT observation).At the basin scale, RMSEn between model outputs and gauge discharges is reduced by 27.11 % (it decreases from 96.71 to 70.49 %) and RMSEn between model outputs and ENVISAT discharges is reduced by 63.28 % (it decreases from 75.10 to 27.58 %).RMSEn f,situ global is high, because a large fraction of in situ stations (25 out of 108) are situated along very small tributaries or at basin heads, where the local RMSEn f,situ is largely over 100 %.These very high RMSEn f,situ have a huge impact on the global RMSEn f,situ global (despite the weighting used to calculate RMSEn f,situ global ).If the statistics are computed using only cells with a RMSEn f,situ below 100 %, we find that the global RMSEn f,situ global is reduced by 14.66 % (it goes from 49.80 to 42.50 %) and the RMSEn f,o global is reduced by 50.21 % (it goes from 51.74 to 25.76 %).This shows the limitation of this assimilation scheme, as ISBA-CTRIP resolution (roughly 50 km by 50 km) does not simulate basin heads well (rivers are too small to be correctly represented in a coarse grid).
Figure 14 displays, for the 12 in situ stations (see Fig. 4 for their locations) already used in Figs. 5 and 8, the mean analysis discharge over the whole experiment time period (red line), which is compared to the free run discharge at the station (blue line), the ENVISAT discharge (green markers) and the gauge discharge (black markers).Overall, analysis discharge is quite close to observed discharges (ENVISAT and in situ).
However, despite the use of the localization, the analysis discharge keeps presenting a quite chaotic behaviour: more     particularly at Sao Paulo de Olivenca (Fig. 14, panel 1), Manacapuru (Fig. 14, panel 2) and during high flow seasons along right-bank tributaries (Fig. 14, panels 7 to 12).This shows the limit of assimilating 35-day repeat period ENVISAT observations.If no data are missing at a given VS, it means that there will be, at most, 11 available observations during 1 year.Moreover, in a state estimation context, only the model output state is corrected and not the model parametrization or, in our set-up, forcings.Therefore, if the model is not constrained by direct or neighbouring observations, it naturally goes back to free run discharge.The performance of the assimilation, with respect to the daily in situ data, is therefore often limited by the low ENVISAT observation frequency.In future works, it will be interesting to study the assimilation of similar data with a higher frequency, such as the JASON-2 altimeter data (which has a 10-day repeat period but a coarser spatial sampling).

Conclusions and perspectives
This study presents, over the Amazon basin, the assimilation of a satellite-derived discharge product into a large-scale hydrological model to correct its state variables.The remotely sensed discharge data are derived from the ENVISAT nadir altimeter and are assimilated into the ISBA-CTRIP model using an ensemble Kalman filter.Five experiments were carried out over the year 2009.For all experiments, the assimilations were able to reduce the modelling errors compared to both observed and gauge discharges.
The first experiments tested different definition of the background error covariance matrices, where the influence of a given observation is either reduced to the only observed cell (SE1-diag), or limited to a few close cells on the hydrological network (SE1-local), or not limited and can potentially impact the entire basin (SE1-direct).Results showed that the complete stochastic matrices gave the best results along the mainstream and the localization treatment appeared neces-sary along the tributaries.The need for the localization is explained by the spurious elements in the error covariance matrix due to the limited ensemble size and the methodology used to generate it.
The last tests compared the corrections of different state variables: the river initial storage (SE1-local), or the river final storage (SE2-local), or the river discharge (SE3-local).The main difficulty with these different types of variables is, on the one hand, the relationship between the control and the observed variables (gathered in the observation operator) and, on the other hand, the reciprocal relationship to generate inputs for the next DA cycle.Results showed that correcting river discharge gives the best global results over the entire basin, as the link between the observed and corrected variables is the most straightforward.Therefore, the final experiment (SE3-local-long) uses the SE3-local configuration over the whole ENVISAT observation period (from September 2002 to 2010) and confirms the possibility of using such low-resolution remotely sensed data in a large-scale model.
These experiments offer several perspectives.First, the localization treatment could be improved by combining the three tested approaches according to the cell's position on the river: discharge correction for cells along the mainstream should be impacted by all upriver observations, while correction for cells on tributaries should be impacted only by close observations along the same sub-catchment.Moreover, the size of the area of influence for a given observations could also vary in time according to the season (high flow/low flow).With ulterior developments of the localization method, new challenges may appear such as the risk of imbalance, already studied in the field atmospheric DA (e.g.Greybush et al., 2011).The analysis of imbalance may need to be considered in future works.
A main limitation of assimilating ENVISAT data is their low repeat period (one observation every 35 days, at best).Indeed, corrected discharges often present strong sudden variations between unobserved and observed dates, as the model goes back to its free run when it is not constrained by an observation.However, there are other satellite altimetry missions with different repeat periods, for example JASON-2 (10-day repeat period from June 2008 to October 2016), JASON-3 (10-day repeat period, launched in January 2016), Sentinel-3A (27-day repeat period, launched in February 2016) or Sentinel-3B (27-day repeat period, which should be launched in 2018).Also, the incoming SWOT (Surface Water and Ocean Topography, launch scheduled for 2021) wide-swath altimetry mission will provide a remotely sensed discharge product.SWOT will have a 21-day repeat period, with an almost global spatial coverage thanks to its two 50 km swaths.All these data could be combined with ENVISAT data (during the overlapping period) within the assimilation scheme to have a denser network of observation over the study domain, to get a better estimate of discharge (similar to a reanalysis) over a multi-decadal time frame (Tourian et al., 2017).
To improve these DA results, several aspects could be investigated.For example, one could study whether a more realistic ensemble method generation could be helpful.In the present study, only the model initial condition and the precipitation forcing are perturbed to generate the background forecast ensemble.More uncertainties in this ensemble could be added by also perturbing CTRIP parameters and/or ISBA outputs.Another DA aspect to look into is the potential use of a smoothing data assimilation algorithm, such as the ensemble Kalman smoother (Evensen and Leeuwen, 2000).A smoother could help to have less "variability" in the corrected discharge.Finally, the assimilation scheme presented in this study could be applied to other river basins in the world, as ISBA-CTRIP is a global LSM.However, more work is needed to apply the DA platform at a global scale.

Figure 1 .
Figure 1.(a) The Amazon basin main tributaries and rivers with the underlying topography from SRTM.(b) Schematic representation of the ISBA-CTRIP system for a given grid cell.ISBA surface runoff (Q ISBA,sur ) flows into the river/surface reservoir S; ISBA gravitational drainage (Q ISBA,sub ) feeds the groundwater reservoir G.The surface water is transferred from one cell to another following the TRIP river routing network.

Figure 2 .
Figure 2. Map of hydro-geomorphological zones defined over the Amazon basin.

Figure 3 .
Figure 3. General framework of the DA method at a kth assimilation cycle.The figure reads from top to bottom and from left to right.The three main variables involved are the river initial storage, the river final storage and the river discharge.M [k−1,k] is the model operator that maps the initial river storage into final river storage, Z k is the diagnostic operator and H k is the observation operator that maps simulated discharge into observed discharge for assimilation.

Figure 5 .
Figure5.Comparison between the ISBA-CTRIP free run (blue line), ENVISAT-derived observed discharges (green markers) and ANA gauge discharges (black dots) over the year 2009.For each panel, the x-axis represents time (in days) and the y-axis represents discharge (in m 3 s −1 ).

Figure 6 .Figure 7 .
Figure 6.RMSEn for the free run simulation compared to the ENVISAT discharges (a) and the gauge discharges (b).

Figure 8 .
Figure8.SE1-direct ensemble mean analysis discharge (red line) compared to the free run discharge (blue line), the ENVISAT observed discharges (green markers) and the measured gauge discharges (black dots) over the year 2009.For each panel, the x-axis represents time (in days) and the y-axis represents discharge (in m 3 s −1 ).

Figure 9 .
Figure 9. (a) Same as Fig. 8, panel 5 but only over the 35 first days of simulation.(b) Location of all active ENVISAT VS on the 25th day of the assimilation (yellow circles) compared to the location of the Serrinha stations (red circle).

Figure 10 .
Figure10.Analysis RMSEn difference between SE1-direct and the SE1-local experiment with respect to (a) the ENVISAT discharge and (b) gauge discharge.Negative RMSEn differences (green colours) mean that the results of the SE1-local experiment are better than the SE1direct results at the given CTRIP cell.Positive RMSEn differences (yellow, orange and red colours) mean that the results of the SE1-direct experiment are better that the SE1-local results at the given CTRIP.

Figure 11 .
Figure11.SE1-local ensemble mean analysis discharge (red line) compared to the free run discharge (blue line), the ENVISAT observed discharges (green markers) and the measured gauge discharges (black dots) over the year 2009.For each panel, the x-axis represents time (in days) and the y-axis represents discharge (in m 3 s −1 ).

Figure 12 .
Figure 12. Analysis RMSEn differences between SE1-local and SE2-local experiments with respect to (a) the ENVISAT discharge and (b) gauge discharge.

Figure 13 .
Figure 13.Analysis RMSEn differences between SE1-local and SE3-local experiments with respect to (a) the ENVISAT discharge and (b) gauge discharge.

Figure 14 .
Figure14.SE3-local ensemble mean analysis discharge (red line) compared to the free run discharge (blue line), the ENVISAT observed discharges (green markers) and the measured gauge discharges (black dots) from 25 September 2002 and for 8 years.In each panel, the x-axis represents time (in days) and the y-axis represents discharge (in m 3 s −1 ).
error standard deviation σ o k,j is then a fraction of the instantaneous discharge.For each virtual station, the value of η o j depends on, first, the deviation from the MGB discharge, noted σ The observation Hydrol.Earth Syst.Sci., 22, 2135-2162, 2018 www.hydrol-earth-syst-sci.net/22/2135/2018/

Table 1 .
Constant values used to generate the background control ensemble X b e,k , the observation ensemble Y o e,k and the model observation ensemble H(X b e,k ).k is the assimilation index and s is the basin zone index.

Table 2 .
Size of the influence area for the localization process.
pared to the two next experiments, SE1-diag and SE1-local, which will highlight the contributions and/or limitations of the chosen localization approach.Finally, the last two experiments, SE2-local and SE3-local, will test the other possible control variables and the reliability of the operator Z k .

Table 4 .
Global statistics for experiments with different localization schemes.

Table 6 .
Emery et al.:Large-scale hydrological model river storage and discharge correction Global statistics for experiments with different types of control variables.

Table 7 .
Statistics between analysis and in situ stations for the different assimilation experiments.Italic values indicate the best result among the three experiments for all tested gages.