Figure 2 . Study area: Purus River Basin. (a) drainage network (in blue), location of the discharge gauge (Canutama, triangle in red), tracks of the spatial altimetry mission Jason 2 (dashed black lines), location of the altimetry virtual station (circle, in black), and the area used for extraction of flood extent (Lower Purus, pink polygons); (b) Hydrological Response Units (Fan et al., 2015); (c) Bare Earth Digital Elevation Model (O’Loughlin et al., 2016); (d) Köppen-Geiger Climatic Zones (Kottek et al., 2006).

Hydrological-hydrodynamic model: MGB

The MGB (“Modelo de Grandes Bacias”, a Portuguese acronym for “Large Basin Model”) is a semi-distributed, hydrological-hydrodynamic model (Collischonn et al., 2007; Pontes et al., 2017). It was chosen for this study because (1) it has been widely and successfully applied in several South American basins (e.g., Paiva et al., 2013; Siqueira et al., 2018); (2) it is representative and similar to other conceptual hydrological models like VIC (Liang et al., 1994) and SWAT; and (3) the hydrological component is tightly coupled to a hydrodynamic routing scheme, allowing the simulation of complex flat, tropical basins. Moreover, the source code of MGB is freely available at www.ufrgs.br/lsh.
Within the model structure, basins are discretized into unit-catchments, which are further divided into Hydrological Response Units (HRU’s) based on soil type and land use. A vertical water balance is performed for each HRU, considering canopy interception, soil infiltration, evapotranspiration, and generation of surface, subsurface and groundwater flows. Soil is represented as a bucket model with a single layer. Flow generated in each HRU is routed to the outlet of the unit-catchment with linear reservoirs. Outflow from each unit-catchment is then propagated through the stream network by using a 1D hydrodynamic model based on the inertial approximation proposed by Bates et al. (2010). The stream network is derived from Digital Elevation Model (DEM) processing. The model has 19 parameters, which are further detailed in the next section. Other model inputs are precipitation, climate data, soil type and land use maps, which are further described in section 2.6 Model Setup.

A priori uncertainty of model parameters

Within MGB model, there are parameters related to vegetation cover (albedo, leaf area index, vegetation height and Penman-Monteith surface resistance), river hydraulics (Manning’s roughness, and width and depth parameters related to geomorphological relationships), and conceptual parameters related to soil water budget (Wm, b, Kbas, Kint, XL, CAP, Wc, CI, CS, CB), which are further detailed in Supporting Information (Table S2). Out of the 19 model parameters, six are fixed and 13 are calibrated.
The a priori uncertainty of MGB model parameters is estimated based on their variability as reported in literature (references in Table S2). Supporting Information (Table S2) presents the calibration parameters, their initial values, range, and the references that support these assumptions.

Sensitivity analysis

In order to understand how different parameter sets (river hydraulic, soil, vegetation) affect model output variables (river discharge, flood extent, river water level, soil moisture, evapotranspiration and terrestrial water storage), multiple model runs were conducted considering four uncalibrated model setups: (1) varying only soil parameters; (2) varying only vegetation parameters; (3) varying only hydraulic parameters; (4) varying all parameters together. One hundred runs were conducted in triplicate to ensure that convergence is not dependent on the initial parameter sets, thus resulting in 300 runs for each setup. In this step, no RS-based calibration is performed yet.
Parameters were varied considering a uniform distribution, and results were analyzed in terms of mean RMSD (root mean square deviation) of each variable, by comparing each run with a reference one (i.e., the initial run with the initial parameter set as defined in Supporting Information Table S2). This was performed in order to understand the sources of model uncertainties related to different sets of parameters (e.g., are flood extent estimates sensitive to vegetation parameters, or are ET estimates sensitive to hydraulic parameters?). The dispersion of model outputs was also compared to uncertainty in the observations, as derived from literature.
To understand which variables are inter-related in the model, we further analyzed the results of setup “(4) varying all parameters together”. This was done by firstly computing the Kling-Gupta Efficiency metric (KGE; Gupta et al., (2009)) between the perturbed runs and a reference one (i.e., run with the initial parameter set) for each variable, and then calculating the Pearson correlation (r) between KGE values for each pair of variables (e.g., discharge and water level, discharge and flood extent, and so forth). This experiment is relevant to evaluate whether two variables get improved or get worsened together, or whether a variable improvement impacts on the deterioration of another. In other words, this approach allows to evaluate the correlation between the variables.

Model setup

The Bare Earth Digital Elevation Model (DEM; O’Loughlin et al., 2016) (Figure 2c) was used for stream network computation and basin discretization with the IPH-HydroTools GIS package (Siqueira et al., 2016). The original DEM resolution is 90 m, and it was resampled to 500 m to facilitate GIS processing. An upstream area threshold of 100 km2 was adopted to delineate the drainage network, and unit-catchments were discretized by dividing the stream network into fixed reach length of 10 km, resulting in 2957 unit-catchments for the whole basin. Soil type and land cover maps were extracted from the HRU discretization developed by Fan et al. (2015) (Figure 2b): (1) deep and (2) shallow forested areas, (3) deep and (4) shallow agricultural areas, (5) deep and (6) shallow pasture, (7) wetlands, (8) semi-impervious areas, and (9) open water, where “deep soils” refer to soils with high water storage capacity, and “shallow soils” are those with low water storage capacity. In the Purus River Basin, 57.4% of the region is covered by forest with deep soils, 26.9% by forest with shallow soils, and 13.7% by wetlands (i.e., river floodplains). Daily precipitation data were derived from TMPA 3B42 (version 7), with spatial resolution of 0.25º x 0.25º (Huffman et al., 2007; available at: <https://gpm.nasa.gov/data-access/downloads/trmm>), and were interpolated with the nearest neighbor method for the centroid of each unit-catchment. Long term climate averages for mean surface air temperature, relative humidity, insolation, wind speed and atmospheric pressure were obtained from the Climatic Research Unit database (New et al., 2000; available at: <http://www.cru.uea.ac.uk/data>), at a spatial resolution of 10’, and also interpolated with the nearest neighbor method.

Model calibration

The MOCOM-UA calibration algorithm (Yapo et al., 1998; Multi-objective global optimization for hydrologic models) was adopted due to its satisfactory performance when coupled with hydrological models (e.g., Collischonn et al., 2008; Maurer et al., 2009; Naz et al., 2014). MOCOM-UA is an evolutionary algorithm, based on SCE-UA (Duan et al., 1992), that simultaneously optimizes a model population with respect to different objective functions. The model population consists of randomly distributed points within the parameter search space, and it reflects the a priori uncertainty of model parameters. Here, the population size was set to 100 individuals. The altered model parameters and their respective ranges are described in Supporting Information Table S2. All calibration experiments are repeated three times with differing initial parameter sets to ensure that convergence is not dependent on the initial parameter sets.
Objective functions to be optimized depend on the calibration setup. In the single-variable calibration, for each variable, three objective functions (OF) that summarize the agreement between simulated and observed (RS) time-series are simultaneously optimized: Pearson correlation (r), ratio of averages (\(\text{\ μ}_{\text{sim}}\ /\ \mu_{\text{obs}}\ \)), and ratio of standard deviations (\(\ \sigma_{\text{sim}}\ /\ \sigma_{\text{obs}}\ \)), which are associated to the individual terms of KGE metric. These 3 objective functions are depicted in Equations 1 to 3, where X denotes the assessed variables (Q, h, A, TWS, ET or W).
\begin{equation} \text{OF}_{1}=\ \left(\frac{\text{\ μ}_{\text{sim}}}{\mu_{\text{obs}}}\right)_{\text{X\ \ }}\left(1\right);\ \ \ \ \ \ \text{OF}_{2}=\ \left(\frac{\sigma_{\text{sim}}}{\sigma_{\text{obs}}}\right)_{\text{X\ \ }}\left(2\right);\ \ \ \text{\ \ OF}_{3}={\ r}_{\text{X\ \ }}(3)\nonumber \\ \end{equation}
For the multi-variable calibration, the objective functions are the KGE of each variable considered: firstly, five objective functions were considered (KGE of all variables except discharge), as depicted in Equations 4 to 8.
\begin{equation} \text{OF}_{1}=\text{KGE}_{h}\ \left(4\right);\ \ \ \ \ \text{OF}_{2}=\text{KGE}_{A}\ \left(5\right);\text{\ \ \ \ OF}_{3}=\text{KGE}_{\text{TWS}}\ \left(6\right);\text{\ \ \ \ \ OF}_{4}=\text{KGE}_{\text{ET}}\ \left(7\right);\text{\ \ \ \ \ \ \ OF}_{5}=\text{KGE}_{W}(8)\nonumber \\ \end{equation}
Secondly, two objective functions were adopted (KGE of selected variables 1 (x) and 2 (y)), as depicted in Equations 9 and 10.
\begin{equation} \text{OF}_{1}=\text{KGE}_{x}\ \left(9\right)\ \ \ \ \ \ \ \ \ ;\text{\ \ \ \ \ \ \ \ \ \ OF}_{2}=\ \text{KGE}_{y}(10)\nonumber \\ \end{equation}
Results are expressed in terms of a Skill Score (S) (Equation 11; Zajac et al., 2017), in order to evaluate the improvement (or deterioration) in the representation of a variable when the model is calibrated with a given variable, compared to the uncalibrated setup.
\begin{equation} S=\frac{\text{KGE}_{\text{calibrated}}-\text{KGE}_{\text{initial}}}{1-\text{KGE}_{\text{initial}}}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (11)\nonumber \\ \end{equation}
KGEcalibrated is the mean KGE resulting from running the model with the calibrated parameters. KGEinitial is the mean KGE resulting from running the model with the a priori parameter sets.

Calibration/Evaluation Data

In the next paragraphs we introduce the data used for model calibration and evaluation, as well as how MGB outputs were evaluated in comparison to them.
-In-situ discharge measurements were obtained from the Brazilian Water Agency Hidroweb database (available at: <http://www.snirh.gov.br/hidroweb/publico/apresentacao.jsf>), at the gauge “Canutama” (code: 13880000; location: S ° 32’ 20.04”; W 64° 23’ 8.88”; drainage area: 236,000 km², period of data availability: 1973 to 2016). Uncertainty in discharge observations can be estimated as ranging from 6.2% to 42.8% at the 95% confidence level, with an average of 25.6% (Di Baldassarre & Montanari, 2009). Discharge was evaluated on a daily basis.
- Remotely sensed water level data were obtained from Jason-2 mission, which presents an orbit cycle of approximately 10 days, and tracks separated by approximately 300 km at the equator (Lambin et al., 2010). It presents an accuracy of approximately 0.28 m (Jarihani et al., 2013), and data are available since 2008. The virtual station presented in Figure 1 corresponds to Track number 165. Processed data for this study were downloaded from the Hydroweb/Theia database (available at: <http://hydroweb.theia-land.fr>). Water level was computed in MGB at the unit-catchment associated to the altimetry virtual station, being an advantage of using the hydrodynamic scheme for flood routing instead of the Muskingum simplification. Simulated and RS water level data were compared every 10 days in terms of anomaly (values subtracted from long term average).
- Satellite flood extent data were derived from ALOS-PALSAR imagery, which presents a ground resolution of 100 m (Rosenqvist et al., 2007). Images were downloaded from Alaska Satellite Facility (available at: <https://www.asf.alaska.edu/>) in processing level 1.5, which already presents geometric and radiometric corrections. A 3 x 3 median filter was used to remove speckle noise (Lee et al., 2014). Images were classified into water (backscattering coefficient less than -14 dB), non-flooded forest (between -14 dB and -6.5 dB), and flooded forest (higher than -6.5 dB) classes, according to Hess et al. (2003) and Lee et al. (2014). The uncertainty of flood extent estimates was estimated based on the RMSE between the resulting classification of this study, and the dual-season mapping developed by Hess et al. (2003). Simulated and RS flood extent data were compared for the pink area depicted in Figure 1, in order to avoid spurious flood extent data in regions that are known to be not subject to flooding. ALOS-PALSAR presents a recurrence cycle of 46 days (from 2006 to 2011), and flood extent data were available and compared to MGB for 21 dates.
- Satellite-based terrestrial water storage (TWS) anomalies were extracted from GRACE mission, launched in March 2002. GRACE provides monthly TWS estimates based on anomalies in gravitational potential, at a resolution of 300-400 km, with a uniform accuracy of 2 cm over the land and ocean regions (Tapley et al., 2004). TWS anomalies were retrieved from three processing centers - GFZ (Geoforschungs Zentrum Potsdam, Germany), CSR (Center for Space Research at University of Texas, USE), and JPL (Jet Propulsion Laboratory, USA), available at <https://grace.jpl.nasa.gov/>, and then the mean value based on the three products was averaged for the whole basin. In MGB, TWS values were computed as the sum of water storage of all hydrological compartments: river, floodplains, soil, groundwater and vegetation canopy. Simulated and RS-based TWS were compared in terms of anomaly (values subtracted from long term average) at a monthly time-scale.
- Satellite-based evapotranspiration estimates were retrieved from the MOD16 product, derived by an algorithm presented by Mu et al. (2011) based on the Penman-Monteith equation. The dataset covers the period 2000-2010 with a spatial resolution of 1 km for global vegetated land areas. Because of that, even though MGB evapotranspiration is calculated for flooded areas (open water evaporation in main channel and floodplains) and vegetation for water balance purposes, only the vegetation-ET output was compared to MOD16. MOD16 products are provided in 8-days, monthly and annual intervals. Monthly intervals were used here and averaged for the whole basin. Accuracy of MOD16 along the Amazon basin is estimated as 0.76 mm/day (Gomis-Cebolla et al., 2019). MOD16 data is available at: < https://www.ntsg.umt.edu/project/modis/mod16.php>. In MGB, evapotranspiration is computed via Penman-Monteith equation, based on the climate input variables.
- Satellite-based soil moisture is derived from the SMOS mission (Kerr et al., 2001), processed by the Centre Aval de Traitement des Données SMOS (CATDS), and downloaded in processing level 4, which combines lower level products with data from other sensors and modeling/data assimilation techniques. The daily L4 root zone soil moisture product at 0-1 m soil depth (Al Bitar et al., 2013) were used (available at: <https://www.catds.fr/Products/Available-products-from-CEC-SM/L4-Land-research-products>), and data from ascending and descending orbits were averaged for the whole basin. In MGB, soil moisture as a saturation degree was computed as the water in the soil compartment divided by the maximum water capacity of the soil (Wm parameter). Since MGB estimates saturation degree values for a soil bucket reservoir, SMOS values were rescaled for the range 0 - 100% according to the Min/Max Correction method described by Tarpanelli et al. (2013) and applied by some studies (e.g., Rajib et al., 2016; Silvestro et al., 2015), and them compared to MGB on a daily time-scale as an average for the whole basin.