The continental part of the water cycle is commonly represented with hydrological models. Yet, there are limits in their capacity to accurately estimate water storage and dynamics because of their coarse spatial resolution, simplified physics, and an incomplete knowledge of atmospheric forcing and input parameters. These errors can be diminished using data assimilation techniques. The model’s most sensitive parameters should be identified beforehand. The objective of the present study is to highlight key parameters impacting the river-routing scheme Total Runoff Integrating Pathways (TRIP) while simulating river water height and discharge as a function of time focusing on the annual cycle. Thus, a sensitivity analysis based on the decomposition of model output variance (using a method called ANOVA) is utilized and applied over the Amazon basin. Tested parameters are perturbed with correcting factors. First, parameter-correcting coefficients are considered uniform over the entire basin. The results are specific to the TRIP model and show that geomorphological parameters explain around 95% of the water height variance with purely additive contributions, all year long, with a dominating impact of the river Manning coefficient (40%), the riverbed slope (35%), and the river width (20%). The results also show that discharge is essentially sensitive to the groundwater time constant that makes up more than 90% of the variance. To a lesser extent, in rising/falling flow period, the discharge is also sensitive to geomorphological parameters. Next, the Amazon basin is divided into nine subregions and the sensitivity analysis is carried out for regionalized parameter-correcting coefficients. The results show that local-region parameters impact water height, while upstream-region parameters affect discharge.
The earth’s climate is undergoing changes in response to natural variability and also increasing concentrations of greenhouse gases (Stocker et al. 2013). The continental part of the water cycle is an active component of this system (Alkama et al. 2008). Global climate change will affect the water cycle, likely creating perennial drought in some areas and frequent flooding in others (Trenberth 2011). It is therefore of prime importance to study and understand these changes to continental waters.
Continental water can be studied at the global scale using land surface models (LSMs) coupled with a river-routing model (RRM). LSMs provide lower boundary conditions to atmospheric general circulation models, while RRMs compute river discharge, which can be potentially used as boundary conditions to ocean general circulation models.
Several RRMs have been developed at a global scale (Vörösmarty et al. 1989; Coe 1998; Hagemann and Dümenil 1997; Oki and Sud 1998; Arora et al. 1999; Olivera et al. 2000; Ducharne et al. 2003; Ngo-Duc et al. 2007; Lucas-Picher et al. 2003; Yamazaki et al. 2011). In those models, the hydrological network is generally derived from a digital elevation model (DEM) at coarse resolution (from 0.25° × 0.25° to 4° × 4°). This network is then used by the RRM to transfer water mass from grid cell to grid cell until the river outlet at the continent–ocean interface. Water mass transfer is based on simple budget equations representing the temporal variations of the water mass stored in each grid cell. The currently available RRMs differ from each other in their method of dealing with river flow velocity, surface parameterization, groundwater, and floodplains.
The description of river flow velocity is addressed in several ways in the literature. For instance, Coe (1998) and Oki and Sud (1998) considered a single surface reservoir with a uniform and constant flow velocity. It is also possible to use a time-constant but spatially distributed flow velocity based on topography and river channel characteristics (Ducharne et al. 2003; Hagemann and Dümenil 1997; Vörösmarty et al. 1989). In contrast, several works rely on a time-varying and spatially distributed flow velocity estimation based on the Manning formula (Manning 1891), for example, Arora et al. (1999), Decharme et al. (2012), Lucas-Picher et al. (2003), and Ngo-Duc et al. (2007).
When a single surface reservoir is considered (Coe 1998; Ngo-Duc et al. 2007; Oki and Sud 1998; Vörösmarty et al. 1989), input water mass to each gridcell reservoir corresponds to surface runoff within the grid cell and river outflow from upstream grid cells. The computed water mass is then transferred to the downstream grid cell along the river network. Several studies (Arora et al. 1999; Decharme et al. 2008; Ducharne et al. 2003; Hagemann and Dümenil 1997; Lucas-Picher et al. 2003) highlighted the necessity of modeling groundwater inflow to the surface reservoir in order to simulate the delayed contribution of groundwater in the river discharge estimation. Inland water bodies such as floodplains are rarely taken in consideration. When it exists, a flooding scheme is triggered when the water stage (Decharme et al. 2008, 2012) or the river discharge (Vörösmarty et al. 1989) exceeds a given threshold and the floodplain within a grid cell exchanges water mass with the associated river reservoir.
Several input parameters are used by RRMs. The parameters can be topographical (such as riverbed slope) and geomorphological (such as river width or roughness) and can also be based on empirical relationships, in which case other parameters are usually required. These relationships correspond to multiplying and powering coefficients linking parameters between each other or with the model state (Leopold and Maddock 1953; Moody and Troutman 2002; Gleason and Smith 2014). Parameters can be either spatially or uniformly distributed over river catchments according to assumptions and simplifications adopted in each study.
It seems obvious that to better represent observations in hydrological modeling, model complexity can be increased (Grayson and Blöschl 2001). In this context, Sieber and Uhlenbrook (2005) point out that sensitivity analysis (SA) can be a powerful tool to both identify the most sensitive (and therefore important) parameters and to understand the hydrological model structure and its response.
SA can be used for two purposes: either to explore multidimensional parameter spaces and understand the source of model uncertainties (Hornberger and Spear 1981) or to identify key parameters and their role in the model response (Saltelli et al. 2008). The first option has been widely investigated, most notably by using the Generalized Likelihood Uncertainty Estimation (GLUE) method (Beven and Binley 1992) both in catchment modeling (Freer et al. 1996) and rainfall–runoff modeling (Blazkova and Beven 2004). The GLUE method assumes that, for a given catchment numerically represented by a given model, several sets of this model’s input parameters can be equally accepted as a simulator of the catchment’s real response. By comparing an ensemble of simulated responses with observations of the same nature, a value of likelihood is assigned to each set of parameters of this ensemble. All simulations giving an acceptable—or behavioral—likelihood are resampled to build a distribution function of the parameter sets. The likelihood values are finally used to estimate the uncertainty associated with the model predictions. Other common methods are regional sensitivity analysis (Wagener et al. 2003) for dynamic identifiability analysis and Bayesian total error analysis method (Kavetski et al. 2006) for calibration and uncertainty estimation. Alternatively, to identify the key parameters, SA can be seen as the study of how uncertainty in the model output can be apportioned to different sources of uncertainty in the model inputs (Saltelli et al. 2008). Methods based on analysis of variance (ANOVA) are particularly suited for this purpose (Francos et al. 2003; Hall et al. 2005; van Griensven et al. 2006). Contrary to GLUE, ANOVA does not require system observation. The model output variance is synthetically generated by considering uncertain parameters as random variables. Using an ensemble of parameters sets, ANOVA determines the contribution of each parameter to the unconditional variance.
Several SAs have been carried out on various RRMs (Apel et al. 2004; Aronica et al. 1998; Bates et al. 2004; Pappenberger et al. 2004; Werner et al. 2005). These studies, mainly based on the GLUE method, were usually performed on short-time events (flood) at catchment scale. To the authors’ knowledge, only Pappenberger et al. (2010) have performed a global-scale SA of the river-routing Total Runoff Integrating Pathways-2 (TRIP2) model (Ngo-Duc et al. 2007), coupled with the Hydrology Tiled ECMWF Scheme of Surface Exchanges over Land (HTESSEL; Balsamo et al. 2009). The study is based on a 10-yr reanalysis run (1986–95) of ECMWF, which was then compared with an ensemble of in situ discharges from more than 400 worldwide river gauging stations provided by the Global Runoff Data Centre (GRDC; http://grdc.bafg.de) by means of a performance measure (Pappenberger and Beven 2004). The ensemble of time-averaged performances was then used to apply both GLUE and ANOVA methods. The sensitivity of crucial parameters for the routing scheme such as slope, river length, groundwater time delay constant, and Manning formula parameters was studied. Results from both GLUE and ANOVA revealed that the groundwater time constant—which models the delayed contribution of drained water into the river—is the most sensitive parameter.
The present study displays an SA using a global method of the TRIP RRM version included in the land surface modeling platform Surface Externalisée (SURFEX; Masson et al. 2013) developed at the Centre National de Recherches Météorologiques (CNRM). For hydrological purposes, TRIP describes river discharge at catchment scale and is coupled to the well-established Interactions between Soil, Biosphere, and Atmosphere (ISBA; Boone et al. 1999) LSM used to simulate thermal and hydrological exchanges between the soil and the atmosphere.
In section 2, we present the ISBA–TRIP version used in this paper along with the Amazon River basin description. Section 3 provides a basic background of variance-based SA methods and details the present ANOVA SA configuration. In sections 4 and 5 we present results for two SA studies differing in terms of number and spatial distribution of perturbed parameters, which are referred to as AMA8 and AMA45. The last section presents our conclusions and perspectives.
2. The hydrological model and study zone
a. The ISBA–TRIP model classical configuration
The ISBA model (Boone et al. 1999) is a relatively standard LSM, including the force–restore method (Blackadar 1976), which is adopted in the current study. The model is globally defined over a regular mesh grid at a 0.5° × 0.5° resolution. The model’s equations are solved in each grid cell independently from the others. All grid cells are only correlated through the spatial patterns of vegetation cover, soil composition, and external forcing, that is, atmospheric (especially precipitation) and radiative inputs. By taking into account heterogeneity in precipitation, topography, and vegetation within each grid cell (Decharme and Douville 2006), ISBA gives a diagnosis of the water and energy budgets in each grid cell (Noilhan and Planton 1989).
The present study relies on the ISBA configuration used in Alkama et al. (2010) and Decharme et al. (2012). The soil is divided into three layers: the superficial layer, the root zone, and the subroot zone (Boone et al. 1999). Precipitation can either fall directly on the soil surface or be intercepted by the canopy. When the canopy is saturated, water drips from the canopy to the soil and then either flows on the surface or infiltrates. Hydraulic conductivity through soil layers is determined following an exponential profile of the saturated hydraulic conductivity (Decharme et al. 2006). The soil water content varies with surface infiltration, soil evaporation, plant evapotranspiration, surface runoff, and deep drainage [for more details, see Alkama et al. (2010) and Decharme et al. (2012)]. ISBA then gives a diagnostic of each water budget component, in particular the surface runoff Qs and gravitational drainage Qsub.
The TRIP RRM is also defined over a regular mesh grid. The present study is at the same resolution of 0.5° × 0.5° as ISBA. Because ISBA works for each grid cell separately, TRIP is dedicated to the lateral transfer of water from one cell to another up to the continent–ocean interface following a river network. The original version of TRIP was developed by Oki and Sud (1998) and consisted of a single linear surface reservoir. Surface runoff data from any LSM were needed to force TRIP and were then converted into river discharge under the assumption of a 0.5 m s−1 constant flow velocity. Subsequent developments by Decharme et al. (2010, 2012) first coupled TRIP with ISBA by taking ISBA simulated surface runoff as input. Then, a groundwater and floodplain reservoir was added and the constant flow velocity was changed to a time- and space-dependent flow velocity based on the following Manning formula:
where s is the channel slope, n is the friction coefficient at the bottom of the channel, and R is the hydraulic radius depending on the channel geometry. In the present work, s, n, and the channel geometry parameters are constant in time but spatially distributed. The actual TRIP model (Fig. 1), consists then in a system of three reservoirs: the surface reservoir S modeling the river, the groundwater reservoir G (Decharme et al. 2010), and the floodplain reservoir F, solved using Eq. (2) through the estimation of the water mass stored in each reservoir (kg):
Only the surface reservoir 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. For any given cell, TRIP inputs are the TRIP outflow from upstream cells (kg s−1) and the ISBA surface runoff for that cell Qs (kg s−1). Moreover, S receives water from groundwater (kg s−1) and can exchange water mass with the floodplain (kg s−1; Fig. 1).
With the implementation of G in TRIP, ISBA Qsub (kg s−1) now flows into the groundwater reservoir, whose outflow goes to river reservoir. This simple modeling assumption does not represent the real groundwater dynamics, but only the delay of drainage contribution to the river. The groundwater outflow (kg s−1) is then estimated by
where τG (s) is the time delay factor. When introduced in Decharme et al. (2010), τG was taken as a constant over river catchment with values between 30 and 60 days.
The floodplain scheme activates when the water height in the river hS (m) exceeds a given critical bankfull height Hc (m). The terms PF, IF, and EF are the precipitation intercepted by the floodplain, the reinfiltration, and the direct free water surface evaporation over the floodplain, respectively. A detailed description of the flood scheme is given in Decharme et al. (2008, 2010, 2012).
b. TRIP specific parameters
Within a cell, S corresponds to an equivalent river that may represent in reality several river branches. The river section is rectangular and is characterized by a slope s (unitless), width W (m), depth Hc (m), length L (m), and a Manning coefficient n (unitless) that quantifies friction as well as channel resistance at the bottom of the river. A more extensive description of these TRIP specific parameters is given below.
The section length is the arc length between the cell center and the direct downstream cell center multiplied by a meandering factor μ, equal to 1.4 in order to take the river meanders into account.
The riverbed slope is deduced from the elevation difference between cells (E − Enext) and based on the Simulated Topological Network (STN-30p) DEM provided at 0.5° × 0.5° resolution by the International Satellite Land Surface Climatology Project, Initiative II (ISLSCP II) database (http://daac.ornl.gov/ISLSCP_II/islscpii.shtml) and the section length as follows:
The river width is globally computed via an empirical geomorphologic relationship between W and the mean annual discharge Qyr (m3 s−1; Moody and Troutman 2002):
with β = 32 (m−1/2 s1/2) for an equatorial basin such as the Amazon River (Decharme et al. 2012).
The river depth is computed from the river width using the following nonlinear function:
The Manning coefficient can be defined as channel resistance to the flow and is generally difficult to estimate. It takes low values (0.025–0.03) for natural streams with deep pools and higher values for small, mountainous streams and floodplains (0.75–0.1) (Maidment 1993). Previous global studies used a global constant value taken equal to 0.035 (Arora et al. 1999). Still, it is commonly assumed that the Manning coefficient should be spatially distributed over an entire catchment. In TRIP, n is therefore spatially distributed according to a simple relationship considering that upstream grid cells contain narrow rivers with high Manning coefficients and that the Manning value decreases as the cells become closer to the river mouth. Arbitrarily, n was chosen to vary between nmin = 0.04 and nmax = 0.06 according to the following linear relation:
with SO being a measure of the relative size of stream in the current cell, SOmax being the same measure at the river mouth, and SOmin (=1) being the measure at source cells.
All these parameters are critical to determine the flow velocity with the Manning formula:
where the river water height (m) is estimated from its actual water mass S and river geometry as follows
where ρW (kg m−3) is the water density.
The system in Eq. (2) is solved using the fourth-order Runge–Kutta method. Once the water storage (kg) of each reservoir is computed, TRIP determines in each cell a diagnosis of the quantities of the river water height (m), following Eq. (9), and the river discharge (m3 s−1), defined as follows:
c. Description of the Amazon basin
In terms of average discharge (2 × 105 m3 s−1) and drainage area (6.15 × 106 km2), the Amazon is the world’s largest river. Its catchment area covers about 40% of South America, and the discharge at its mouth represents 30% of total freshwater inflow to the Atlantic Ocean (Wisser et al. 2010).
As shown in Fig. 2a, the river source is located in the Peruvian Andes. The river flows through the Brazilian tropical rain forest and receives water from several important tributaries: the Japurá River, the Purus River, and the Negro River (16% of the total discharge) at Manaus. 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; Molinier et al. 1993).
From a geological perspective, the Amazon basin can be divided in three major morphostructural units: the shields at the eastern part of the basin (Guiana shield to the north and Brazilian shield to the south), the central Amazon trough, and the western Andean Cordillera. In terms of vegetation, the basin is covered by tropical rain forest (71%) and savannas (29%; Sioli 1984).
The whole Amazon basin is located in a humid tropical climatic zone. The central basin is under an equatorial climate zone, implying high surface temperatures, air humidity, and precipitation. Thus, a vast floodplain along the mainstream is filled every year, leading to the damping of discharge extremes. Northern and southern parts of the basin are under a tropical climate with a dry and a wet season. However, because of the seasonal shift of the intertropical convergence zone, the maximum rainfall season 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.
3. SA with the ANOVA method
This section presents the SA method used in this study.
a. SA generalities
Variance-based SA methods aim to analyze the model output variance against an ensemble of simultaneously modified parameters (Efron and Stein 1981; Sobol 1993) and are efficient even for nonlinear and nonmonotonic models (Saltelli and Bolado 1998). The goal is to estimate a set of sensitivity indices (Si; one for each parameter) that represent the contribution of the associated parameter in the model output unconditional variance. The two most usual measures of sensitivity are the main effect and the total effect. The main effect is used in factor privatization searching for the parameter that controls the most model output uncertainty. The total effect, on the contrary, is used in factor fixing searching for irrelevant parameters that can be constrained to an arbitrary value (Saltelli et al. 2008).
There are several methodological variants to calculate partial variances, for example, Sobol’s method (Sobol 1993) or the simple and extended Fourier amplitude sensitivity test (FAST and EFAST; Cukier et al. 1973; Fang et al. 2003). The main difference between FAST and the Sobol method lies in the approach used to calculate the variances multidimensional integrals. Whereas Sobol’s method uses a Monte Carlo integration procedure, FAST uses a pattern search based on sinusoidal functions. However, Reusser et al. (2011) tested the three sampling methods (FAST, EFAST, and Sobol) and found comparable results between each method. Even though these methods are quite popular, they still have a high computational cost by requiring a large number of model runs to estimate Si.
Currently, some studies focus on direct calculation of the double-loop multidimensional integrals underlying Si by trying to reduce computational cost (Jansen 1999; Saltelli et al. 2010; Sobol 1993). Alternatives have been found to overcome dependence on the problem of dimensionality (linked to the number of parameters) and thereby reducing computational cost, in particular through high-dimensional model representation (HDMR) and metamodeling (Rabitz et al. 1999; Ratto et al. 2007).
In addition, Sieber and Uhlenbrook (2005) and Reusser et al. (2011) introduced temporal dynamics of parameter sensitivity (TEDPAS) and time series of grouped error in order to perform temporal SA of any model output variable such as river discharge. This kind of analysis aims to identify dominant parameters during a given hydrological process and highlight parameter interactions. For example, Garambois et al. (2013) performed a temporal Sobol’s sample ANOVA to identify key parameters of the physically based and spatially distributed Model of Anticipation of Runoff and Inundations for Extreme Events (MARINE) during flash flood events. Guse et al. (2014) used a temporal FAST–ANOVA method on the Soil and Water Assessment Tool (SWAT) model to detect which parameter dominates in poor model performance cases.
In this study, the state-dependent parameter (SDP) metamodeling approach developed by Ratto et al. (2007) is used. Also, a temporal sensitivity analysis following Garambois et al. (2013) and the Reusser and Zehe (2011) TEDPAS methodology is applied. The ANOVA is based on Sobol’s sampling method as in Ratto et al. (2007).
b. The ANOVA formulation
The formulation of the ANOVA decomposition is based on Sobol (2001). The model output Y can be defined as a function of k independent uncertain input parameters :
where is the model operator and the parameters are defined over , the set of all possible values. It is ensured that the range of variation of these parameters corresponds to the assumed parameter uncertainties and, in the following SA studies, leads to realistic values for these parameters .
Now, let be a possible set of normalized parameters (i.e., scaled between 0 and 1). The relationship that links the normalized model parameters with Y is denoted by
where f is also the model operator but taking the normalized parameters as inputs.
It is assumed that Y can be decomposed in a sum of real functionals with increasing input dimensionality of the following form:
Under the hypothesis that f is integrable over the unit hypercube of size k, it is assumed that the following integral holds , ,
The variable can take any real values, so the above integral implies that the values have a zero mean. It follows from Eq. (14) that the functionals are orthogonal and can be expressed as an integral of . Indeed,
where signifies that the integration is performed over all parameters but . Successively,
Moreover, if f is also square integrable over the unit hypercube of size k, so are the summands and
Now, is considered as a random variable uniformly distributed over the unit hypercube of size k, meaning that each parameter follows a continuous uniform law between 0 and 1. The model output and all the summands , , and are random variables as well. This implies that integrals of f and represent expectations and integrals of the variances of and .
For instance, the model output expectation is given by
and the model output conditional expectations given first-order factors are defined as
meaning that the expectation of the model output Y is computed over all possible values of keeping fixed.
A similar definition holds for higher-order factors:
The total variance (Var) of the model output Y is then formulated as follows:
is called the main effect of parameter on Y representing the expected reduction of total variance if could be fixed and where
is the covariance caused by combined effects of parameter and on Y. Two parameters are said to interact when their effects on Y cannot be only explained by and , and so on for higher-order terms.
They represent the part of a given single parameter (or a given subset of parameters) variance in the unconditional variance. They could be seen as an estimate of the model sensitivity to this particular parameter. From Eq. (23), it appears that
In practice, calculating the sensitivity indices requires first the evaluation of the integrals in Eqs. (24) and (25) along with the double-loop integral in Eq. (22). This evaluation is achieved stochastically using an ensemble of model output realizations. First, the ensemble of input parameter sets , with the size of the ensemble, are generated using a quasi-random sequence generator. As each parameter is independent from others and normalized, each parameter follows a continuous uniform law between 0 and 1. Then, the model output ensemble is established by running the model f with each one of the input parameter sets.
Second, the are calculated using the “SDP metamodel” approach framed by Ratto et al. (2007). It is designed to estimate functionals until second order, that is, and in Eq. (13). This is done with the Smoothing Spline ANOVA (SS-ANOVA) routine (available at https://ec.europa.eu/jrc/en/econometric-statistical-software#SS-ANOVA-R). Once the terms are estimated for the entire ensemble, the associated are straightforwardly estimated using the estimator suggested in Doksum and Samarov (1995). The advantage of this method is its adaptability to any sampling method of the input ensemble, as it provides fast, accurate, and unbiased results. The SS-ANOVA procedure provides (besides an estimation of the coefficients) the standard errors of SDP and hence the relative significance of estimated HDMR terms. Finally, a last important routine output is the metamodel correlation coefficient measuring the likeness between the original and emulated model. In the present study, only first-order Si are estimated and studied. Model output variance explained only with first-order Si characterizes a purely additive model and will be confirmed by an close to 1 and low interaction effects (quantified by ). Otherwise, it will indicate nonnegligible interaction effects. For more details, refer to Ratto et al. (2007).
To sum up, the ANOVA method produces a set of sensitivity indices taking value between 0 and 1 and associated with a single or a subset of input parameters. These indices represent a percentage of the total variance.
Two SAs with different choices of parameters (but all with the same general configuration) have been carried out in this study and are detailed in the following subsection.
c. Simulation description
In this study, the SA is based on a sampling of the parameter space of size Ne = 1024 drawings [this value was chosen following Garambois et al. (2013) and has reasonable computational costs]. To generate the initial ensemble of normalized input parameter sets , the quasi-random sequence proposed by Sobol (1967) is used. These sequences use a base of two to successively form finer partitions of the unit interval and then reorder the coordinates in each dimension. They spawn the entire unit hypercube for any dimension (the computational routine can generate until 230 − 1 points of maximum dimension 52); and thereby the entire range of definition for each input parameter is explored.
A model output ensemble element consists of daily time series over the period 2008–10. The temporal SA is conducted over the three years. The years 2006 and 2007 are used as spinup to avoid effects of the initial condition and transitive states in the model output variance, but no SA is conducted over these two years. The year 2009 was a particularly wet year, as mid-June floods reached their highest level for 50 years. Conversely, in 2010, the Amazon River experienced an extreme drought and reached its lowest level for half a century. The year 2008 is considered an average year. Therefore, studying the model sensitivity for the full 3-yr period from 2008 to 2010 is relevant for all types of hydrological years.
ISBA–TRIP is run in offline mode and uses atmospheric forcing from the Global Soil Wetness Project phase 3 (GSWP3; http://hydro.iis.u-tokyo.ac.jp/GSWP3). This project consists of three global-scale experiments with the objective of investigating long-term changes of 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 Twentieth Century Reanalysis (Compo et al. 2011).
ISBA–TRIP estimates water levels and discharges at each of the 2028 cells contained in the Amazon basin mesh grid. For reasons of computational cost, the temporal SA is detailed for a selection of cells corresponding to known in situ stations (Fig. 2b). Table 1 briefly describes these evaluation cells.
To carry out this temporal SA, two different sets of parameters have been chosen and presented in the next section.
d. Choice of parameters
1) Global SA over the Amazon basin
The first step is to analyze TRIP model sensitivity to all parameters over the entire basin. All TRIP parameters are considered: the river and floodplain Manning coefficient n (unitless) and (unitless), the groundwater time constant (days), the riverbed slope s (unitless), the river width W (m), the river bankfull depth (m), and the river section length meandering coefficient μ (unitless). See Table 2 for an overview of these parameters’ nominal range.
Apart from the groundwater time constant and the meandering ratio, which are constant, all other parameters are spatially distributed over the entire basin. To avoid overparameterization, parameter fields are perturbed by applying a constant multiplying factor. Then, sensitivity to a given field is studied via this multiplying factor. Moreover, particular attention is paid to the perturbation of the riverbed slope field. Indeed, following Pappenberger et al. (2010), the riverbed slope is perturbed by first applying a powering constant and then a multiplying factor . Multiplying and powering factors have a nominal value of 1 and are defined in a ±ε% interval around it. The ε value is chosen so that the perturbation range is representative of the parameter’s uncertainty (Gatelli et al. 2009). Existing parameter value tables for the Manning coefficient (Maidment 1993), previous studies (Pappenberger et al. 2010; Decharme et al. 2010; Paris et al. 2016), and comparison to a remotely sensed optical-image-based database for river width (Yamazaki et al. 2014) were used to fix ε for each parameter, which led to the range of parameter listed in Table 3.
This first SA study is denoted by AMA8—“AMA” for Amazon (though it is planned to extend the ANOVA method to other basins in future work) and “8” since eight perturbed parameters are considered. Note that except for the riverbed slope (for which two coefficients are used), the sensitivity to any of the other parameters presented in Table 3 will be directly referred to as the sensitivity to the corresponding TRIP parameter (see fifth column in Table 3); for example, sensitivity to will be directly mentioned as the sensitivity to W. When discussing slope sensitivity, the use of and will be preserved.
2) Regional SA over the Amazon basin
In the first SA study, a preliminary insight into TRIP model sensitivities is drawn. Based on these results, a second SA study, referred to as AMA45 in the following, is carried out and aimed at regionalized sensitivities (disregarded in the first SA, since parameter perturbations are made at the entire basin scale for AMA8). A subset of input parameters considered as the most critical ones for understanding TRIP model behavior is then considered per geographical zone. For this purpose, the Amazon basin is divided into nine hydro-geomorphological zones. The interest of such a subdivision is to study the separated impact of the parameters of each zone. These zones were designed according to 1) hydrological arguments as the main course is separated from the tributaries that have their own zones and 2) geological arguments as three major morphostructural units are distinguishable.
The nine zones (Fig. 3) are the following: 1) the upstream Andean part of the basin until the city of Iquitos, Peru; 2) the main stream from Iquitos to Óbidos; 3) the main stream from Óbidos to the river mouth; 4) left-bank tributaries from the Napo River to the Japurá River, including 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.
The procedure of AMA45 is to select a reduced number of important parameters from AMA8 results and rerun the SA by considering each one of these parameters in any of the nine Amazon basin zones. The overall SA configuration will remain identical to the first SA in terms of number of members within the ensembles (Ne = 1024), using atmospheric forcing (GSWP3) and configuration of the ISBA–TRIP runs (three years from 2008 to 2010 with spinup from 2006 to 2007).
4. Results and discussion for AMA8
a. Water level sensitivity
Figures 4–6 display the resulting first-order sensitivity indices for TRIP water level [see Eq. (9)]. Table 4 shows the time average over the 3-yr study period of the metamodel and sum of all higher-order effects (i.e., the interaction effects) for each cell. The close to 1 indicates a good convergence of the method, and interaction effects close to 0 indicate that the main effects (first-order sensitivity indices) are sufficient to explain the water height variance, characterizing purely additive contribution of the parameters.
The simulated water height ensemble for each evaluation cell generally presents a large spread with amplitude up to 20 m for some cells such as Jatuarana and Óbidos (Fig. 4). The ensemble first and ninth deciles (plotted in black solid and dashed lines, respectively) are clearly distant all year long, even during low-flow season (when the ensemble dispersion is at its lowest). It can be already pointed out that water level presents an important sensitivity to TRIP parameters.
Overall, the sensitivity time series show no clear interannual variations, and four parameters are systematically activated: first, the river Manning coefficient (purple line), closely followed by the riverbed slope powering coefficient (green line), then the riverbed width (orange line), and the groundwater time constant (blue line). In some cases, the river section length meandering coefficient is slightly activated (red line in Figs. 5b, 6). In most figures, , , and explain 40%, 35%, and 20%, respectively, of the model output variance; and explain a small percentage of the sensitivity also during low-flow period. To summarize the sensitive parameters by order of importance are the following:
Three different behaviors are observed. First, evaluation cells along the Amazon mainstream—Itapéua, Jatuarana, and Óbidos (Fig. 4)—show very light interseasonal variations, implying that only three parameters’ Si are significant: and are quasi-constant and presents a smooth increase when the water level is at its lowest. According to Eq. (9), depends on the actual river reservoir storage and on the channel geometry (its width and length). The dominating parameters n, , and W are directly involved in the Manning formula [Eqs. (1) and (8)] that estimates the flow velocity. The flow velocity itself appears in Eq. (2) to determine the river outflow, thereby impacting the surface water storage. Therefore, these parameters represent the contribution of the S to .
Second, the left-bank-tributary cells—Acanaui, Caracarai, and Serrinha (Fig. 5)—present a similar behavior, but sensitivities to n, , and W are noisier, which can be related to precipitation events. Indeed, left-bank tributaries present smaller drainage areas and are therefore more sensitive to local events, namely precipitation. Besides, peaks are more numerous here, but, as for the general result cells (Fig. 4), they appear when the water level ensemble spread is very narrow and suddenly increases. It is worth noting that for the strongest peaks a smaller but synchronized peak is observed. This implies that during sudden regime changes, the reservoir geometry controls the water level dynamics.
Third, the right-bank-tributary cells—Canutama, Manicoré, Itaituba, and Belo Monte (Fig. 6)—present parameter sensitivity with clear interseasonal patterns: and mainly dominate as observed in previous results but decline slightly in low-flow season. During this period, distinctly appears and even overpasses other Si. Eventually, by the end of the low-flow season, and have peaks similar to those observed in Fig. 5, which occur when the water level rises again. The activates exclusively during low-flow season. This parameter represents the water inflow from the groundwater reservoir into the river. It is a continuous source of water in the river, but the ANOVA method, as used here, determines sensitivities as a fraction of the overall variance. Then, the groundwater storage contribution takes a higher proportion during the low-water season when other contributions (precipitation, surface runoff, etc.) impacting the water level are lowered. This is undoubtedly due to the very low water level during this season compared to the rest of the year. This seasonal pattern in and the high contribution of is more pronounced for right-bank tributaries than left-bank tributaries and is explained by the large-amplitude difference between the ensemble maximum and minimum values (almost 20 m).
b. Discharge sensitivity
Figures 7–9 display first-order sensitivity for TRIP discharges, and Table 5 presents the time average over the study period for and sum of all interaction effects. The metamodel coefficients are slightly lower than those obtained for water levels but generally remain above 0.85, showing again a good convergence of the method. However, discharge interaction effects account for 10%–15% of the unconditional variance (e.g., first-order effects explain 85%–90% of the variance).
Contrary to water level ensemble, the simulated discharge ensembles (gray area) present a quite narrow spread with tightened first and ninth deciles (plotted in black solid and dashed lines, respectively) close to ensemble extrema. The same ensemble of input parameter sets has been used to generate both the water height ensemble and the discharge ensemble (see section 3d). But, contrary to the quite large water height ensemble dispersion, the discharge ensemble dispersion is narrower. A first outcome is that discharge is much less sensitive to TRIP parameters.
In similarity with water height results, the discharge sensitivity displays no remarkable interannual sensitivity patterns. The dominating parameters are the same as those for water height but with a different ordering, that is,
The parameter is the most dominating parameter, and may exceed 0.9 (e.g., Fig. 9) and has its highest values during low- and high-flow seasons. Parameter is the parameter for groundwater flow exfiltration to the drainage network (Roux et al. 2011) and represents the mass transfer from the groundwater into the river. This high sensitivity is in accordance with the rainfall–runoff modeling literature and with the conclusions of Pappenberger et al. (2010). In a stable regime, the discharge is mainly driven by water mass transfer from the upstream river and groundwater.
During the transition period between high- and low-flow season, and also in response to precipitation events (Fig. 8 mostly), sensitivity drastically drops and other parameters are sensitive. Figures 7–9 display an evident anticorrelation between and all other parameter Si. For example, time- and station-averaged correlation is −0.89 and −0.87 between and and between and , respectively. On the contrary, those other parameter sensitivities are highly correlated, for example, 0.94 and 0.97 for the time- and station-averaged correlation between and and between and , respectively.
Between high- and low-flow seasons, or during precipitation events, the discharge is driven by regime changes, that is, by flow velocity changes. Thus, parameters involved in the Manning formula activate, respecting the ordering observed in water level sensitivities, that is,
In addition, the river section length, represented by μ, is also quite an important parameter for the discharge because 1) it is directly used to compile discharge [Eq. (10)] and 2) it is an important river geometry parameter during regime changes.
From this first study, where the parameter sensitivities are studied identically over the entire Amazon basin, five parameters appear to be relevant in explaining TRIP sensitivities: the groundwater time constant, the river Manning coefficient, the riverbed slope powering coefficient, the river width, and the river section length meandering ratio. These parameters translate the predominant impact of both the groundwater reservoir and the use of the Manning formula for flow velocity estimation. Also, it appears that water level is essentially sensitive to parameters describing the geomorphology of the reservoir, while discharge is more sensitive to mass transfer. As discharge seems to be weakly sensitive to TRIP parameters—deduced from the low dispersion of the simulated ensemble—it is rational to assume that discharge is sensitive to other parts of the ISBA–TRIP system such as precipitation forcing, ISBA parameters, and ISBA outputs.
It is of interest to conserve the five dominating parameters and to study the impact of regionalized parameters over TRIP water height and discharge. The next SA keeps the same configuration, but a different set of five parameters will be taken for each subbasin zone introduced in section 3d (Fig. 3). Giving a total number of 45 parameters (five parameters for nine zones), this simulation will therefore be denoted as AMA45.
5. Results and discussion for AMA45
Concerning the AMA45 experiment, results are shown in Figs. 10–13. The same color code as AMA8 is used to plot parameter sensitivity time series (e.g., purple line for n and orange line for W) and to differentiate two parameters of the same nature but from two different regions; different panels are used for each region.
Results are then discussed for a fewer number of cells than AMA8: Caracarai, Belo Monte, Serrinha, Itaituba, and Óbidos (stations 6, 10, 7, 9, and 3, respectively, in Fig. 2b). These five locations are representative of the results obtained for all cells in this experiment. For Caracarai (in zone 5) and Belo Monte (in zone 8), the pixels are situated in the same zone as their contributing upstream pixels. Serrinha and Itaituba are located in zone 2 and zone 3 (respectively) while their upstream contributing pixels are only located in one different zone–zone 5 and zone 8 (respectively). Finally, Óbidos’ cell is the first cell of zone 3 along the Amazon mainstream and receives flows from zones 1, 2, 4, 5, 6, and 7.
a. Water level sensitivity
The results at Caracarai and Belo Monte (not shown here) are similar to those in the AMA8 experiment (Figs. 5, 6). This was expected as they are located in the same zone as their drainage area and therefore have the same configuration as AMA8.
Contributions from other zones are observed at Serrinha and Itaituba (Fig. 10). At these cells, the activated parameters are the local (from the cell’s zone) geometrical (W and μ) and morphological (n and ) parameters, and the parameter from the upstream zone (5 for Serrinha and 8 for Itaituba) is activated. The overall tendencies are the same as in AMA8. Observations are identical at Óbidos (Fig. 11), but with from all upstream zones being activated. Therefore, materializes the water inflow in the cell.
With the new perspective of regionalized parameters, it appears that the local reservoir geomorphology monitors the river water level by controlling the amount of mass that leaves the cell. Meanwhile, during low-flow season, water level is sensitive to water mass inflow from upstream areas.
b. Discharge sensitivity
Similar to water level discussion, discharge sensitivities for AMA45 at Caracarai and Belo Monte (not shown) are equivalent to those for AMA8 (section 4b; Figs. 8, 9). Then, results at Serrinha and Itaituba (Fig. 12) show the integrator behavior of discharge as only upstream zones are activated. Similarly, at Óbidos (Fig. 13), all upstream zones are activated.
However, a particular behavior is observed at Óbidos. Water from the drainage area necessarily flows through zone 2 before reaching Óbidos. Therefore, all farthest upstream zones (1, 4, 5, 6, and 7) are activated through while geomorphological parameters are activated in the closest zone (zone 2). In zone 2, is slightly activated as well, but because the drainage area of the zone is smaller than others, its impact is less important. Focusing on , two types of peaks are observed. These peaks correspond to low-flow season for the highest peak (high-flow season for smaller peaks), but with some delay due to the transfer time between the actual high- and low-flow season and the time when the water reaches the observed pixel.
In conclusion, concerning discharge sensitivity, represents the continuous water inflow while the other geomorphological parameters represent regime changes (cf. AMA8 in section 4b).
From the results above, it seems important to have a better knowledge of geomorphological parameter values to improve water height and discharge diagnostics in ISBA–TRIP. Currently, the river width is defined using an empirical relationship, but it could be directly estimated from existing recent databases such as North American River Width Data Set (NARWidth; Allen and Pavelsky 2015) or Global Width Database for Large Rivers (GWD-LR; Yamazaki et al. 2014). Similarly, other DEMs could replace the currently used DEM in ISBA–TRIP (GTOPO30) to get a better estimate of the riverbed slope, such as Shuttle Radar Topography Mission (SRTM) edited in Hydrological Data and Maps Based on Shuttle Elevation Derivatives at Multiple Scales (HydroSHEDS; Lehner and Grill 2013) for latitudes under 60°. Additionally, the upcoming Surface Water and Ocean Topography (SWOT; Alsdorf et al. 2007; Fjortoft et al. 2014; Biancamaria et al. 2015) mission will also provide such kinds of information on the river width and the river surface slope, which could be used to set TRIP parameters. Optical images could also be used to spatialize the meandering ratio. A better definition of those parameters would limit the model parameter uncertainty. The most uncertain parameter would therefore be the Manning coefficient and the groundwater time constant, which still remain difficult to estimate.
These SAs are preliminary works in preparation for model reanalysis and dynamic parameter estimation through data assimilation (Moradkhani 2008; Reichle 2008; Pedinotti et al. 2014) that use in situ or remotely sensed water level and discharge data to correct model parameters and/or state. Among them, the incoming SWOT mission offers a valuable potential to improve global-scale hydrological modeling. The idea is to use SWOT satellite products in a parameter estimation configuration to correct the main TRIP parameters highlighted during the SAs.
Additionally, this study focused on the water height exclusively, which is a physical variable of interest, for example, for water depth data assimilation. However, current satellites provide water elevation of top water body distance to a reference geoid or ellipsoid. Therefore, we quickly investigate the contribution of model parameters with respect to water depth anomaly (i.e., obtained by subtracting the averaged water height over the study time period to the water height). Sensitivities were calculated as previously with the same perturbation ranges. Preliminary results for anomalies are presented at Óbidos in Fig. 14 with the associated for the time-averaged water height. While the overall behavior is quite the same for the Manning coefficient, the riverbed slope, and even the river width (anticorrelated behavior with n and s), clearly activates for short periods when the anomaly ensemble is near to zero. It appears that the time-averaged water height (used to calculate the anomalies) takes all variance from the parameters n, s, and W, letting dominate when water heights are close to the averaged water height (i.e., when anomalies are close to zero). Those results will be further developed and studied in future works.
The present study focused also only on river-routing model parameters and was applied over a unique river basin—the Amazon—but it offers numerous perspectives. In a global perspective, as SWOT will observe all rivers wider than 100 m between 78°S and 78°N, the present platform will be extended to other river basins situated in other climatic zones, such as the Mississippi or the Niger. In addition, it will be of great interest to study the impact of the LSM (here the ISBA model), the atmospheric forcing (more precisely precipitation), and even the initial reservoir states on TRIP outputs to improve our understanding of the continental part of the water cycle.
One important final remark about this SA is that the following results are quite dependent on the chosen parameters along with their perturbation range, but also on the chosen model itself. Indeed, the TRIP model considers only the kinematic wave propagation equation for the river reservoir. Other studies include diffusive wave propagation equation (Yamazaki et al. 2011; Winsemius et al. 2013) in their routing models along with a finer description of the topography and the flood dynamics. Even 2D-type finescale hydrodynamic model has been applied at continental scale (Sampson et al. 2015). The ANOVA mathematical formalism can be easily exported to these models and will probably give different results because of the different model physics (e.g., kinematic wave against diffusive wave) and may lead to potentially different sensitivity patterns in space and time.
This study aims to analyze the ISBA–TRIP large-scale hydrological model sensitivity over the Amazon River basin. An output model variance decomposition method was used to identify key river-routing model parameters during a 3-yr period (2008–10).
Two analyses were carried out to evaluate the sensitivity of model parameters at different spatial scales. The first study (AMA8) considered parameters whose uncertainty was defined at the entire catchment scale. The second study (AMA45) used the same parameters but with regionalized uncertainty according to a geological and hydrological division of the Amazon River basin. The objective was to separate the local and upstream impacts of the parameters on both water height and discharge.
For AMA8, the results showed no interannual sensitivity. For both water height and discharge, the importance of river Manning coefficient and riverbed slope was highlighted. These observations are consistent with the use of the Manning formula to estimate the flow velocity. Also, river width and length had a nonnegligible influence. Contrary to water height, the groundwater time constant significantly dominated discharge sensitivity.
The second study used the same parameters but with regionalized uncertainty according to a geological and hydrological division of the Amazon River basin. The aim was to separate contribution from local parameters and upstream parameters.
Water level sensitivity is relatively constant through time, with the Manning coefficient multiplicative constant and the riverbed slope powering coefficient explaining 40% and 35%, respectively, of the unconditional variance. The river width is less present by taking only 20% of the variance. However, when the water level is very low, a nonnegligible impact of the groundwater time constant is observed. Also, in response to sudden regime change and precipitation events, the river width and length may be important. Overall, water level is essentially sensitive to geomorphological parameters. The second regionalized SA indicates that the local geomorphology of the reservoir drives the water level most of the year. In a low-flow period, the inflow from upstream zones is materialized by the activated upstream and turns into an important contributor to water height as local phenomena (e.g., triggered by precipitation events) are minimized.
Discharge sensitivity presents several temporal patterns: seasonal patterns linked to the alternation of high- and low-flow seasons and short-term patterns associated with precipitation events in upstream regions. The groundwater time constant dominates, sometimes at over 90% of the variance. In a permanent regime, this parameter represents the main water inflow into the river. Geomorphological parameters activate during transitional regimes. It is worth noting that, because of the narrow dispersion of simulated discharge ensemble, discharge is actually weakly sensitive to TRIP parameters. A reasonable assumption would be to consider precipitation and even the ISBA configuration (e.g., soil conductivity distribution and vegetation cover) as the main drivers of discharge. The second regionalized study confirms that discharge is mainly driven by mass transfers, in particular upstream mass transfers.
Ultimately, this SA method offers an extensive variety of experiments. The SA results depend on both the studies model and the parameters range. Therefore, the same SA formalism applied with a different set of parameters or another model may give a different behavior of the sensitivity indices. Keeping the ISBA–TRIP model, one could study the impact of the precipitation forcing or the ISBA parameters on TRIP outputs or preserve the same formalism on other river basin. The method could also be applied to other models with a different physics and different parameters.
This work was performed using high performance computing resources from Calcul en Midi-Pyrénées (CALMIP; Grant 2015-P1408). The authors are very thankful to M. Ratto and I. M. Sobol for supplying the open-source SS-ANOVA and LPτ routines, respectively. The GSWP3 team is acknowledged for letting the authors use their different forcing fields. The authors acknowledge funding from the CNES Terre–Océan–Surfaces Continentales–Atmosphère (TOSCA) committee. C. M. Emery is supported by a CNES/Conseil Régional de Midi-Pyrénées Ph.D. grant.