The Surface Water and Ocean Topography (SWOT) mission will provide free water surface elevations, slopes, and river widths for rivers wider than 50 m. Models must be prepared to use this new finescale information by explicitly simulating the link between runoff and the river channel hydraulics. This study assesses one regional hydrometeorological model’s ability to simulate river depths. The Garonne catchment in southwestern France (56 000 km2) has been chosen for the availability of operational gauges in the river network and finescale hydraulic models over two reaches of the river. Several routing schemes, ranging from the simple Muskingum method to time-variable parameter kinematic and diffusive waves schemes, are tested. The results show that the variable flow velocity schemes are advantageous for discharge computations when compared to the original Muskingum routing method. Additionally, comparisons between river depth computations and in situ observations in the downstream Garonne River led to root-mean-square errors of 50–60 cm in the improved Muskingum method and 40–50 cm in the kinematic–diffusive wave method. The results also highlight SWOT’s potential to improve the characterization of hydrological processes for subbasins larger than 10 000 km2, the importance of an accurate digital elevation model, and the need for spatially varying hydraulic parameters.
Remote sensing from spaceborne platforms is increasingly used for the monitoring of components of the hydrological cycle, including river discharge (Santos da Silva et al. 2010). The surface soil moisture can be observed by the Soil Moisture Ocean Salinity (SMOS), Advanced Scatterometer (ASCAT), and Soil Moisture Active Passive (SMAP) satellites (Pierdicca et al. 2013; Kerr et al. 2010; Flores et al. 2012). The Gravity Recovery and Climate Experiment (GRACE) satellite provides terrestrial water storage variations by measuring large-scale gravity fluctuations over time (Syed et al. 2008; Landerer and Swenson 2012). Several altimetric satellites have been launched in the past to measure water surface elevations. The missions that have observed river free surfaces are ERS-1 (1991), TOPEX/Poseidon (1992), ERS-2 (1995), Jason-1 (2001), Envisat (2002), Jason-2 (2008), and the Satellite with Argos and Ka-band altimeter (SARAL; 2013) (Biancamaria et al. 2010; Santos da Silva et al. 2010). They provide information at the global scale, even over ungauged basins. However, they do have limitations, principally, their long revisit time (between 10 and 30 days; Biancamaria et al. 2010) and their coarse spatial resolution: the first nadir altimeters were unable to observe rivers less than 1 km wide. Additionally, their track spacing was poor for hydrology. Recent missions have been able to observe rivers 100 m wide (Santos da Silva et al. 2010), but with the same long revisiting time.
The Surface Water and Ocean Topography (SWOT) mission (launch planned for 2020)1 is a swath mapping radar interferometer designed to measure spatial and temporal water elevation changes in lakes, reservoirs, and large river channels over continental surfaces between 78°S and 78°N (Durand et al. 2014; see also https://swot.jpl.nasa.gov). With regard to rivers, the fundamental SWOT measurements will consist of river water surface elevation, slope, and width. When it comes to water surface elevation, the accuracy will depend on the average area. Errors are expected to be 10 cm (within one standard deviation σ) over an area of 1 km2 inside the river mask (e.g., for a 10-km reach of a 100-m-wide river). Errors will be lower if the average area is larger, and higher for narrow rivers, because of additional difficulties in the determination of the water mask. The revisit times will depend on latitude and will be around four revisits per 21-day-orbit repeat period at midlatitude. In addition to the above-mentioned products, a global estimate of discharge (at the time of the observation) will be produced for all rivers wider than 50 m (SWOT Project 2014), and algorithms for this purpose are currently being developed and tested (Durand et al. 2010, 2014; Gleason and Smith 2014). This data will complement the existing gauge network by providing data between gauges and over ungauged rivers.
In terms of its contribution to the understanding of the continental water cycle, SWOT data are anticipated to be used in conjunction with the above-mentioned remotely sensed data, in situ data, and models of varying complexity. The reach-averaged SWOT products are likely to be used by the vast majority of these models. Raw pixel-by-pixel data may be used only for some specific local studies because of the large amount of data needed, the complexity, and the large associated error.
In addition to classical model evaluation, assimilation is a promising way to use these data and fosters new research on assimilation techniques and model development. Andreadis et al. (2007) used the LISFLOOD-FP hydrodynamic model (Bates and De Roo 2000) and synthetic observations of water elevation to estimate river discharge over a 50-km reach of the Ohio River. Durand et al. (2008) estimated bathymetric depth elevation at five points over a 240-km reach of the Amazon River using synthetic water surface elevations and an ensemble-based data assimilation. Pedinotti et al. (2014) assimilated synthetic data over the entire Niger basin using the large-scale hydrometeorological model Interactions between Soil, Biosphere, and Atmosphere (ISBA)–Total Runoff Integrating Pathways (TRIP). In the first study, the model states were updated, while in the two others, invariant parameters were estimated. Both types of application may be used in the future. The possibility of assimilating the discharge product instead of the surface elevation product must be further investigated depending on the application and the model complexity. While the use of discharge can be envisaged for simple applications, uncertainties in the discharge algorithm and the model are likely to interact in advanced applications (floods, inundations, and low flows). None of these studies has focused on the production of the initial runoff, except by perturbing or imposing it, but the question will obviously be raised for applications modifying the model state (a bias coming from rain or the water surface budget must be corrected in conjunction with the river water depending on the basin time characteristics).
The many models that simulate water elevations can be divided into two main categories. At a large scale, the water surface elevation simulation is mainly motivated by the simulation of large flood plains, the dynamics of which are essential in order to estimate the discharge of some major rivers of the world correctly. Hence, the LISFLOOD-FP model has been used over the lower Ob (Biancamaria et al. 2009) using observed discharge as incoming flow and a 1° × 1° ISBA–TRIP run for lateral transfers, and Hydrological Modeling and Analysis Platform (HyMAP) has been used to simulate the Amazonian basin with a 0.25° × 0.25° resolution (Getirana et al. 2012). At the global scale, models such as Catchment-based Macroscale Floodplain (CaMa-Flood; Yamazaki et al. 2013) and ISBA–TRIP (Decharme et al. 2012) use a relatively low resolution for the land surface scheme that produces the runoff and finer-scale DEM information to characterize the floodplain areas within a grid. Second, at a very fine scale, numerous hydrodynamic models (over several tens of kilometers) were developed for various applications, mainly in relation to flooding.
SWOT will be the first altimetric mission to document intermediate (or regional)-scale basins with a relatively high frequency (e.g., for temperate regions such as western Europe: 50 000–200 000 km2, every 3–10 days; Pavelsky et al. 2014). This offers a new opportunity to link water surface elevation, land surface processes, and meteorology more closely at this scale.
Several issues must be addressed in order to evaluate SWOT’s potential ability to provide information on the water cycle at the regional scale. On the model side, it is important to be prepared to evaluate models based on instantaneous outputs (including discharge and free water surface elevations) rather than the usual daily averages and to evaluate them using proxies of the future SWOT products (simulators or finescale hydrodynamic models). It is also crucial to validate the routing schemes at this intermediate or regional scale and verify whether the relevant parameters for hydrodynamics can be estimated accordingly. For SWOT to be used to improve simulations of hydrology (and the surface water budget), SWOT errors on river elevation, width, and slope must also be compared with errors of current hydrological models at this scale. Finally, research is needed to evaluate the scale of hydrologic or meteorological phenomena that can be documented by SWOT using hydrological models and well-instrumented test beds. The methodologies developed in such test beds could be transferable to regions of the world with less in situ data coverage.
The Garonne catchment in France (56 000 km2) is among the basins that could be used as a test bed for SWOT studies, particularly for evaluating the coherence between the hydrological and meteorological scales. In this basin, the river width is less than 200 m and the inundations are limited in extent, occurring in conjunction with heavy large-scale precipitation events or snowmelt. A regional hydrometeorological model, ISBA–Modèle Couplé (MODCOU), was established using the finescale Système d’Analyse Fournissant des Renseignements Atmosphériques à la Neige (SAFRAN; Habets et al. 2008) and the surface parameter database for land surface models ECOCLIMAP2 (Faroux et al. 2013). Two reaches of the rivers are covered by two high-resolution hydraulic models, which would allow discussion of the scaling issue, and 97 gauging stations are available. However, ISBA–MODCOU, which uses the Routing Application for Parallel Computation of Discharge (RAPID) scheme (David et al. 2011a,b), is only validated with daily discharge.
This paper seeks to contribute to the establishment of a test bed in the Garonne basin within the SWOT framework. Its objectives are to further improve the routing scheme of the model to produce output that can be compared to SWOT; to discuss the level of detail needed for hydrodynamics, including the spatialization of the hydrodynamic data; and to evaluate the results at the basin scale through comparison with high-resolution hydraulic models. More specifically, the routing scheme of ISBA–MODCOU is further developed in order to validate the model on the basis of both discharge and river depth. Then, the application of such a model within the framework of SWOT is discussed.
Section 2 describes the hydrometeorological model together with the input and validation data. Section 3 describes the experimental setup, including the various routing methods tested in order to simulate river depths. The results for both discharge and water depth are presented in section 4. A discussion of the results within the SWOT framework, including issues related to the choice of river routing, the implications for river depth and water elevation simulation accuracy, and the relevance of SWOT to the documentation of hydrologic processes is proposed in section 5.
The hydrometeorological model used to produce results in this paper is a distributed regional-scale model composed of the land surface model ISBA and the hydrological MODCOU. In this version, the original river routing scheme of MODCOU has been replaced by RAPID, as in David et al. (2011a,b), in order to allow a simulation of discharge in all the river reaches of the model. The various components of the model are described below and summarized in Fig. 1.
The ISBA land surface model (Noilhan and Planton 1989) within the Surface Externalisée (SURFEX) platform (Masson et al. 2013) is used to simulate the physical variables in the upper soil, soil surface, and vegetation and to simulate water and energy exchanges within the soil–surface–atmosphere continuum. Its parameters are derived from the ECOCLIMAP2 ecosystems and surface parameters database (Faroux et al. 2013) at 1-km resolution. ISBA is forced in our study by the SAFRAN meteorological analysis (see section 3b). In this study, the ISBA–diffusion (DIF) configuration is used: a multilayer approach is employed to solve the one-dimensional Fourier law and the mixed form of the Richards equation explicitly to calculate the time evolution of the soil energy and water budgets (Boone et al. 2000; Decharme et al. 2011, 2013). This version describes the soil using 14 layers; the hydrological active depth depends on the vegetation. A subgrid runoff scheme (Habets et al. 1999a) is employed to account for subgrid heterogeneities for precipitation, surface parameters, and soil wetness. The surface runoff over saturated areas, or Dunne runoff, is computed using the Variable Infiltration Capacity model (VIC; Dümenil and Todini 1992; Wood et al. 1992; Zhao 1992; Habets et al. 1999a) in which the saturated fraction of the grid cell depends on soil moisture, precipitation intensity, and a shape parameter B fixed at 0.5 (Habets et al. 2008). Orography is derived from the Shuttle Radar Topography Mission (SRTM) 90-m model (Farr et al. 2007). ISBA is run at the 8 km × 8 km regular gridcell resolution.
The hydrological and hydrogeological model platform MODCOU computes the spatial and temporal evolution of the piezometric level of multilayer aquifers using the diffusivity equation (Ledoux et al. 1989) and routes the continental surface water into the river. The surface runoff simulated by ISBA is transferred to the river by the isochrone transfer (ISO) module (Ledoux et al. 1984) and then routed within the river by the RAPID module (David et al. 2011a,b). The bottom runoff (or drainage) is passed to the aquifer module, which exchanges with the river. In the original SAFRAN–ISBA–MODCOU (SIM)-France model (Habets et al. 2008), the river flow is computed every 3 h, and the evolution of the aquifer is computed daily. However, in the Garonne basin, the alluvial aquifers are not simulated by MODCOU. The drainage is routed similarly to the surface runoff. Figure 1 is an illustration of the ISBA–MODCOU hydrometeorological model.
The subsections below describe the river network component, the transfer to the river, and the routing in the river.
1) Transfer to the river
The ISO module (Ledoux et al. 1984; Habets et al. 2008) transfers the surface (overland) and bottom runoff (deep soil drainage) to the river. In this way, the runoff partitioned by the production function (ISBA) is routed to the river network. Each drainage area is divided into a number of isochronal zones equal to the number of time steps necessary for the flow to reach the nearest river cell. The transfer times depend on the topography and concentration time, which is a parameter to be fitted (Habets et al. 2008). The transfer time ttra (s) from one cell to the neighboring downstream cell is given by Eq. (1):
where Δl (m) is the distance between the centers of two cells; Sl (m m−1) is the slope between the two cells; SDA (m2) is the accumulated drainage area; and β is a calibration parameter, usually taken to be equal to 0.25 (Habets et al. 1999b). This routing was run at the daily time step to compare daily simulated discharges with daily observations and at the 3-h time step to compare simulated 3-h discharges with 3-h observations, or 3-h discharges averaged over the day with daily observations.
2) River routing
River discharges are found by the parallel-computing-based RAPID (David et al. 2011a,b). The spatial resolution of river grid cells considered by RAPID is 1 or 2 km. Unlike the initial routing scheme of MODCOU used in Habets et al. (2008), RAPID computes the routing in the entire river network using a matrix-based Muskingum routing scheme instead of calculating flows in a small, predetermined number of cells, usually corresponding to the gauging stations. Hence, RAPID increased the possibilities for scientific evaluation and improvements at a similar computing cost.
Even if the Muskingum method can lead to good-quality results, as shown by David et al. (2011a,b), in our case it has severe limitations. The flow velocity is constant whatever the regime, levels are not simulated, and backwater and floodplain storage effects are not taken into account. River models, which have been further improved to use more detailed routing schemes based on the kinematic or diffusive wave, develop floodplain inundation schemes (Getirana et al. 2012; Yamazaki et al. 2013) and use the full Saint-Venant momentum equations (Paiva et al. 2011) or approximations that are stable enough to allow longer time steps (Yamazaki et al. 2013), such as the local inertial equation (Bates et al. 2010).
This study focuses primarily on the simulation of river depths and flow velocity. David et al. (2011a) proposed a method to calibrate the Muskingum parameters based on observations. In addition, Saleh et al. (2011) used rating curves distributed on a 188-km river network in northern France to estimate river depth within the framework of a river aquifer exchange study. The latter method relies on distributed hydraulic parameters that must be estimated. Based on the same parameters, a kinematic wave routing scheme [based on Decharme et al. (2010)] and a diffusive wave scheme [based on an improved Muskingum scheme (Todini 2007)] are also tested.
The derivation of the original Muskingum approach (Gill 1978) is based on Eqs. (2) and (3) written for a river reach without lateral flow. In each reach, the water flow is computed with a time-invariant velocity. The main advantage of this method is that it is not necessary to know the hydraulic parameters of the river network. The principal disadvantage is that the water velocity is time invariant, a hypothesis that is not true in nature because the flow velocity in a river channel increases with discharge:
where S (m3) is the volume of water stored in the reach; k (s) is the time for water to be transferred between two reaches; ε (dimensionless) is the space-weighting factor ranging from 0 to 0.5, as a function if the storage in the river is controlled by the downstream conditions (ε = 0), or if the inflow and the outflow have exactly the same weight (ε = 0.5); and I (m3 s−1) and O (m3 s−1) are the upstream (inflow) and downstream (outflow) discharges of the reach.
Finally, the outflow of the reach at time step t + Δt, where Δt is the time step of the model (1800 s in the present study), can be written as
Knowing the hydromorphology of the river channel and the equations that relate the discharge to the river depth (see next section), the river depth in Muskingum routing is calculated from the discharge at each time step using the Newton–Raphson method (Todini 2007).
(ii) Manning–Strickler kinematic wave
The derivation of the original Manning–Strickler kinematic wave (MS) method is based on Eqs. (2) and (7), which were both written for a river reach without lateral flow (Decharme et al. 2010). This method is equivalent to the one used in some global hydrological models (e.g., Alkama et al. 2010; Pedinotti et al. 2012):
where V (m s−1) is the mean flow velocity in the river channel and L (m) is the length of the reach.
The water storage is calculated as a function of using a Runge–Kutta fourth-order scheme to prevent numerical bias given by the nonlinearity of Manning’s formula (Decharme et al. 2010). The Manning equations are given by Eq. (8) for rectangular channels and Eq. (9) for trapezoidal channels:
where Kstr is the Manning–Strickler factor, which quantifies the roughness of the riverbed; Rh or Rh* is the hydraulic radius [the ratio between the wetted area A or A* (m2) and the wetted perimeter P or P* (m) of the river channel] for rectangular and trapezoidal channels, respectively; So (m m−1) is the slope of the riverbed; and W (m) is the river width for rectangular channels. For trapezoidal channels, Bo (m) is the horizontal length of the bed and α is the angle between the banks of the bed and the vertical plane. The variables υ and υ* (m s−1) are the flow velocities related to Rh and Rh*, respectively, and h and h* (m) are the river depths for rectangular and trapezoidal channels, respectively, between the river free surface and the bed. They are calculated as a function of the water storage in the reach and the channel geometry as shown in Eqs. (10) and (11):
Once has been calculated, is deduced as a function of the difference between and , as shown in Eq. (12):
To avoid numerical instabilities, the flow velocity in all grid cells of the river network must be lower than a maximum velocity Vmax given in Eq. (13):
In this study, a time step of 300 s was chosen. This time step corresponds to a maximum flow velocity of 6.7 m s−1 for 2-km reaches located in the major plain rivers, which is always higher than the flow velocity in the case of high discharges.
Ponce and Yevjevich (1978) extended the original Muskingum method to time-variable parameters: the time transfer and the space-weighting factor vary with time. This routing was inspired by the Muskingum–Cunge method (Cunge 1969). Unlike the original Muskingum method, the variable parameter Muskingum–Cunge suffered from a loss of mass, which increased with the flatness of the bed slope. Todini (2007) resolved these inconsistencies and proposed a modified routing scheme, hereafter referred to as the Muskingum–Cunge–Todini (MCT) scheme, by introducing a celerity c (m s−1), distinct from the flow velocity. The expression of the outflow calculation written in Eq. (14) is close to the Muskingum formulation, and Eq. (5) is still valid. However, the difference is that C1, C2, and C3 are not time invariant (Todini 2007):
where is a first guess of . The following equations can be repeated two or more times to obtain an optimal expression of .
At the end of the guess loop, Eq. (15) is computed:
The three time-variant coefficients C1, C2, and C3 are calculated in Eq. (16):
The coefficient β at times t and t + Δt is
The celerity is calculated as a function of the average discharge equal to , and of the hydraulic parameters of the river channel (Todini 2007).
The Courant number C* at times t and t + Δt is
The cell Reynolds number D* at times t and t + Δt is
As in the Muskingum formulation, the river depth in MCT is calculated using the Newton–Raphson method (Todini 2007) to extract the river depth from the discharge. The calculation of the river depth is computed by from the Manning–Strickler equations, given that the discharge and the hydraulic parameter values in every river reach of the catchment are known.
To avoid numerical instabilities, the flow velocity is limited as a function of the time step, as shown in Eq. (13). For the same reason as in the kinematic wave routing (see previous section), the value of Δt in the MCT formulation is 300 s.
c. Hydraulic models
The 1D model developed by the Institut de Mécanique des Fluides de Toulouse (IMFT; Larnier 2010) is used for hydraulic simulation of the upstream part of the Garonne catchment (an 80-km reach centered on the gauging station Verdun-sur-Garonne, see Fig. 2). It uses a 1D finite-difference hydrodynamic scheme to compute water depth, velocity, and discharge in a channel flow. The model is forced by river discharges observed at Toulouse, 30 km upstream of the gauging station of Verdun-sur-Garonne. The spatial step of the river transects described in the 1D shallow water model was between 500 and 1000 m, and 153 measured river transects were used to build the model geometry.
MASCARET is a one-dimensional, free-surface hydraulic model based on the Saint-Venant equations used for modeling flood events, submersion waves resulting from the failure of hydraulic infrastructures, the regulation of river infrastructures, and the propagation of canal waves. It was developed by the Laboratoire National d’Hydraulique et Environnement (LNHE) at Électricité de France (Goutal and Maurel 2002). MASCARET is used for the downstream part of the Garonne (a 60-km reach around Marmande). The 1D hydraulic model is forced by river discharge observed at the gauge of Tonneins (near Le Mas d’Agenais), and 83 measured river transects were used to build the model geometry. By comparing simulations of the model and observations at Marmande, the efficiency for discharge and RMSE for river free-surface elevation are 0.98 and 0.26 m, respectively.
3. Model setup and experimental design
a. The Garonne catchment
The Garonne River basin (56 000 km2) is located in southwestern France and drains the northern slopes of the Pyrenean chain (along the French border with Spain). The Pyrenees and the Massif Central mountains border the basin to the south and the east, respectively. Its main tributaries are the Ariège, Dordogne, Tarn, and Lot Rivers. The Garonne and the Ariège flow from the Pyrenees, while the Dordogne, Tarn, and Lot flow from the Massif Central. The climate over the basin is influenced by oceanic conditions over the western part of the domain. It is also characterized by heavy rainfall events during winter and relatively warm weather during summer. Hydrological data (observations of mean annual discharge, mean winter discharge, and mean summer discharge) of the main river gauges located along the Garonne (Fig. 2) are given in Table 1. The annual average discharge in the river gauges of Marmande and Le Mas d’Agenais (downstream Garonne) is about 500 m3 s−1, but values of more than 2000 m3 s−1 can be observed during flood episodes, and less than 100 m3 s−1 during severe droughts. The Garonne basin is highly impacted by human activity: hydropower dams are present in the mountainous areas, while a high number of small farm dams and some reservoirs of intermediate capacity were built in the plain area for irrigation purposes and to sustain low water flows. However, the impact is limited in the downstream Garonne, except during low-flow periods. The impact is higher for the tributaries (in summer and during the snowmelt period). Weirs for navigation or regulation are also present in the basin (especially along the Dordogne River) and can locally have an influence on the hydraulics of the river. It must be noted that RAPID is not able to explicitly take into account these weirs.
b. The meteorological data
The meteorological data were provided by SAFRAN (Quintana-Seguí et al. 2008). The following eight physical variables are analyzed: 2-m air temperature, 2-m relative humidity, 10-m wind velocity, cloudiness, incoming solar radiation (short waves), atmospheric incoming radiation (long waves), rainfall, and snowfall. SAFRAN computes vertical profiles of temperature, humidity, wind speed, and cloudiness every 6 h for 615 climatically homogeneous zones across France. In our study, the first estimates for these profiles usually come from the ECMWF operational archives, at resolutions decreasing from 25 to 20 km over the period, and are refined with surface observations through an optimal interpolation method. A precipitation analysis is performed daily based on a first estimate derived from climatological fields. All analyzed values are then interpolated at the hourly time step, and solar (visible) and infrared radiation are calculated using a radiative transfer scheme (Ritter and Geleyn 1992) that uses vertical profiles of temperature, humidity, and cloudiness. The hourly distribution of precipitation is inferred from the analyzed hourly specific humidity and further constraints from the snow–rain transition elevation (Quintana-Seguí et al. 2008; Vidal et al. 2010). Atmospheric variables are ultimately projected on the regular 8-km grid used by ISBA. In the analysis, 400 stations are used for precipitation and around 100 are used for the other surface parameters.
c. Evaluation datasets
The river depth and discharge observations (the time measurements of which are variable) used in this study were obtained from the flood forecasting service [Service de Prévision des Crues (SPC)] for 97 hydrological stations in the Garonne catchment (see Fig. 3). At these stations, discharges were calculated as a function of the river depth using the operational rating curves of the SPC.
Data could be obtained for only 19 of these stations for river depth validations. The rating curves were reprocessed to determine the water depth above the riverbed, assuming a rectangular bed. The roughness coefficient Kstr and the difference between the riverbed and the origin of the operational scale were deduced by fitting the rating curve using the Manning–Strickler equations. Hence, the river depths used in the following were calculated as the difference between the water elevation and the bed elevation obtained from the rating curve processing.
d. Determination of the hydrological model parameters
Special attention was paid to the determination of the RAPID parameters for the routing methods tested in this study (Table 2). The detailed description and the method of determining the parameters for each method are given below.
1) Muskingum parameters: Transfer time and space-weighting factor
(i) Muskingum original formulation
In the original Muskingum method (MD11), the parameters k and ε determined by David et al. (2011a) were used. David et al. (2011a) computed k in all MODCOU reaches of the Garonne catchment by using topographic information: the slope of the riverbed and the upstream catchment area of the relevant reach.
(ii) Muskingum based on the lagged cross correlation
A detailed analysis of the simulations based on the original parameters proposed by David et al. (2011a) showed that discharges were not well phased with observations in time. The lagged cross-correlation method (David et al. 2011b) led to the calculation of an optimal value of k, corresponding to the flow velocity during floods. The goal was to maximize the correlation factor between two hydrographs (discharge observations) of two stations located in the same portion of the river. Once the transfer time maximizing the lagged cross correlation was determined, the flow celerity and k could be calculated (assuming that the distance between the two stations was known).
A total of sixteen river portions were used to calculate optimized transfer times in the Garonne catchment. The lagged cross-correlation method was applied using hourly observed discharges over the full study period (1995–2006), but the results have a low sensitivity on the period used, provided the period is longer than 6 months.
2) Hydraulic parameters
The following three hydraulic parameters: bed slope, roughness factor, and river width must be carefully evaluated in order to simulate accurate river depths. The methods for determining them are described below.
(i) Bed slope
The bed slope (Fig. 4) was estimated by using the topography database SRTM 90 m (Farr et al. 2007). Given the 1- or 2-km river gridcell resolution of RAPID, which is lower than that of SRTM, a specific method was developed. To obtain positive downward slopes over the river (needed by the routing schemes), the representative elevation of the RAPID grid cell was chosen to be equal to the minimum elevation of SRTM over this cell. Then, the original slopes were smoothed over large river portions to obtain downward slopes in all reaches of the Garonne catchment. The slopes were averaged over reaches bounded by the confluences of the river network in the model. The fact that SRTM gives surface elevation and not riverbed elevation is discussed further in section 5.
(ii) River width
Following Leopold’s original method [Leopold and Maddock (1953); see also Arora and Boer (1999)], the river width was obtained in every reach through the relationship between the average discharge over a chosen period and the river width. We associated the width with the average discharge (during the 1995–2006 period) on 20 reaches where we had obtained the values of measured width and observed discharges. A logarithmic regression (Fig. 5) was performed over these 20 points. Assuming that the values of the average discharges simulated by the original Muskingum routing (that we know in all reaches of MODCOU) were not too highly biased, it was possible, after some preliminary tests, to determine an acceptable river width in every reach of the study area. The 20 observed river widths were not used in the model to avoid spatial discontinuities of simulated river depths (and thus surface slope between two grid cells) over a river reach, given that there is a strong dependence of the simulated river depth on the river width value. The final relation between the river width and the average annual discharge is given by Eq. (20). Its coefficient of determination is 0.84:
(iii) Manning–Strickler factor
Because the roughness of the riverbed has a greater effect on small river flow than on large mouth flow, a linear relationship was taken between Kstr and the river width (Arora and Boer 1999). The two selected minimum and maximum values of the Manning–Strickler coefficient, and , are equal to 10 (representative of mountain rivers) and 40 (representative of plain rivers). The minimum and maximum river widths simulated in the catchment are and [see Eq. (20)], and i is the index of the reach considered.
The value of Kstr was corrected for reaches where rating curves were available. The rating curves were obtained from the SPC. Operational river discharges were estimated in situ by relating river depth measurements taken nearly continuously to periodic measurements of flow velocity and channel cross-sectional area, from which instantaneous river discharges were derived.
In this case, the rating curve was approximated by the Manning–Strickler equation by supposing that W and So were known in every reach. For large rectangular channels, the Manning–Strickler equation (discharge as a function of the river depth) is
The roughness coefficient Kstr was optimized so that the root-mean-square error (RMSE) between the O(h) curve and a rating curve was minimized. As we did not have a rating curve for every grid cell of the Garonne catchment, it was necessary to carry out spatial interpolation between the points where we had rating curves. This means that the value of Kstr evolved linearly along the river channel between two reaches where Kstr was calculated using the rating curve.
e. Experimental design
The period from 1 August 1995 to 31 July 2006 was chosen for the validation of the hydrological MODCOU over the Garonne catchment. The simulations of MODCOU were forced by surface and bottom runoffs produced by ISBA–DIF (Decharme et al. 2013). RAPID transferred the water in the river channels by four different routings (see Table 2). RAPID was initialized with no water in the river channels at the beginning of the simulation period; river reaches received water after 1–5 days of simulation, depending on their location in the Garonne catchment. RAPID was forced by ISBA-simulated runoff at daily or 3-hourly time steps. Discharge and river depth results could be averaged over the day or every 3 h. Validations were performed with daily or 3-hourly discharge and river depth observations, using mainly the Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe 1970) and the discharge ratio for discharge evaluation and RMSE for river depths evaluation. The time step of the original data was variable. The three-hourly and daily discharges and river depth data were calculated by averaging instantaneous 1-hourly data obtained by the interpolation of the variable data. Georeferenced simulations of RAPID were compared with river gauge observations at locations that should correspond to that of the grid cell under consideration in the model. A rectangular shape was used for the river channels in the simulations (see section 5 for a discussion of the impact of the use of a trapezoidal channel). Special attention was paid to the six main stations listed in Fig. 2 and Table 1: these stations are located in the plain area of the basin, where comparison with the hydraulic models is possible and the SWOT observations will be of good quality.
a. Discharge validations: Daily time step
In this section, RAPID is forced by ISBA at a daily time step, and discharge outputs are averaged over a period of 24 h. Figure 7a shows the performance of the original version (MD11) of David et al. (2011a). Over the basin, the NSE varies from −14 to 0.76. At the outlet (Le Mas d’Agenais), the NSE is equal to 0.70. The NSE is higher over some river portions (mainly in the Dordogne basin, in the north of the domain), while an NSE lower than zero can be seen over very small and highly human-impacted rivers (reservoirs and water uptake for irrigations) or in the upper Garonne (in the Pyrenees). In the latter case, the poor results were attributed to difficulties in snowmelt simulations and the presence of dams. Other results of intermediate quality (NSE < 0.5) in some other tributaries in the plain area (mainly south of Le Mas d’Agenais) were attributed to the presence of dams and water uptake for irrigation.
1) Impact of improved transfer times for the Muskingum routing method
The calibration of the celerity in 16 river portions [section 3c(1)] led to improved efficiency at six stations, as shown in Fig. 7. The improvement could be mainly observed in the plain area of the Garonne catchment, where the transfer times were underestimated by David et al. (2011a). Given that the lagged cross correlation gave a celerity for high discharges, the best phasing was observed when discharge values were high. For medium and low discharges, the simulated celerity was slightly overestimated and the water arrived at the downstream stations too early. As a result (considering low, medium, and high discharges), the efficiency was improved because the NSE is mainly controlled by high discharge values. At Le Mas d’Agenais, the daily NSE over the 1995–2006 period improved from 0.70 (MD11) to 0.83 in the Muskingum method based on the lagged cross-correlation (MLCC) run.
At a few river gauges of the Dordogne, Isle, Tarn, and Aveyron Rivers, the lagged cross-correlation method degraded the scores in comparison with MD11. These poor results were attributed to celerity variations within the reaches used to calculate the model parameter k. Usually, the celerity increased from upstream to downstream areas. Thus, the typical consequence was an overestimation of the celerity in the upstream part of the reach, so the discharge phasing was degraded at a gauge located in the middle of the reach. Second, we showed that the lagged cross-correlation ρ in the Dordogne and Isle Rivers was not very sensitive to the lagged time. This meant that the accuracy of the maximum value of the lagged cross correlation was not as high as it was for the portions in the Garonne River. The related hypothesis was that the presence of dams in the river network could have an influence on the flows, disturbing the natural discharge. Hence, the propagation of flows in the Dordogne could not be accurately simulated by the routing schemes used in this study.
2) Manning–Strickler kinematic wave method
The application of the parameters Kstr, So, and W determined in section 3c(2) improved discharge simulations in comparison with MD11 on the main rivers of the basin (Table 3). In the total study area, the scores of 44 river gauges were improved and 53 degraded in comparison with MD11, but many of these scores were only slightly degraded. When considering bigger changes, 17 gauges were improved by 0.05 or more, while the NSE decreased by more than 0.05 on four gauges. The best improvement of the NSE was observed in the plain area of the Garonne River (Fig. 7), where hydraulic parameters were relevant. When considering the results by river flow classes in the plain area (Table 4), MS consistently improved the results for high discharge over both MD11 and MLCC. For low flows, the phasing between MS simulations and observations was poor because discharges were underestimated (underestimation of water produced by ISBA or human impacts such as water uptake or release). Underestimated discharges led to underestimated flow velocity and poor phasing.
However, Table 4 shows that in Portet-sur-Garonne and Verdun-sur-Garonne, medium discharges were better simulated in MLCC than in MS. As both MLCC and MS underestimated the value of the maximum discharge (not shown), the poor phasing of the MLCC maximum discharge at Portet-sur-Garonne and Verdun-sur-Garonne resulted in an improvement in the scores of the middle quantile range, which corresponds to the recession period after the peak (Fig. 8a).
3) Muskingum–Cunge–Todini method
Like MS, the MCT diffusive wave routing scheme improved discharge scores in the Garonne catchment in comparison with the original formulation, MD11 (Table 3). MCT was run in order to test whether the formulation of the diffusive wave was able to improve the scores in comparison with the MS formulation. In general, MCT did not improve the scores (Fig. 7), especially for small slopes. The diffusion resulted in a diminution of the vertical extension of the hydrograph, so the maximal discharge value of a flood was decreased. As the flow velocity is a function of the discharge, the reduction of the maximum discharge value tended to decrease the flow velocity. Thus, the water in the channel flowed too slowly, and the phasing between simulations and observations was degraded. Here, it was decided, after some preliminary tests, to limit the minimum value of the slope to 0.001 in the calculation of the diffusion factor D* [see Eq. (19)]. Note that all the results presented here include this limitation. This point is further discussed in section 5.
b. River depth validations: Daily time step
As the river depth in a river channel is a function of the discharge, the improvement of discharge scores resulted in an improvement in the water depth scores (Fig. 9). In MD11, the RMSE between simulations and observations in the Garonne catchment generally ranged between 30 and 70 cm. In the three new routing methods (MLCC, MS, and MCT), the RMSE was reduced by 10–20 cm in the Garonne River channel (Table 5). The differences between MCT and MS were zero or around 1–3 cm in most cases. Similarly to the discharge scores [section 4a(3)], MCT was slightly worse than MS because the diffusion in MCT degraded the phasing between simulations and observations in downstream areas. The NSE calculated on river depths (not shown) in the downstream Garonne River were about 0.7–0.8 for MD11, and they are improved by 0.10–0.15 with the other routings. The NSE values of MCT were slightly degraded in comparison with the MS routing (by 0.01–0.03). The NSE values for small rivers were often negative, confirming the model’s poor performance at simulating the river depth dynamics in these areas of the basin. The apparently good RMSE results over these rivers must be compared to SWOT errors that can be significantly higher for rivers in mountains than for rivers in plains because of difficulties in the determination of the river mask for smaller rivers and additional layover errors due to the surrounding topography. Hence, the practical interest of SWOT to improve our knowledge on these rivers will be highly limited.
c. Discharge and river depth validations: 3-h time step
In this section, RAPID is forced by ISBA at a 3-h time step (instead of a 24-h time step) and the simulated river discharges are averaged at a 3-h time step. This is the first attempt to compare ISBA–MODCOU–RAPID simulations at a time step smaller than 24 h. This comparison is very important in the context of the SWOT mission, which will provide instantaneous measurements.
First, the 3-h simulated discharges were averaged at a daily time step and compared to daily observations. Whatever the routing scheme, the NSE for discharges was in most cases about 0.01–0.02 worse, and the RMSE for river depths was about 1–2 cm higher. The scores were, in fact, very close to the scores of the previous section. Given the values, the degradation can be considered insignificant. There are several possible causes for the slight degradation, which must be further examined: the SAFRAN hourly data present some discrepancies (Quintana Seguí et al. 2008), it is very difficult to validate the transfer time within the soil of ISBA, and the transfer time in ISO is very simple as it depends only on the slope and area.
When 3-h simulation results were compared with 3-h observations, the results showed that the NSE for discharges was about 0.05–0.10 lower in comparison with daily simulations because of the intraday discharge variability. Apart from that, 3-h simulations showed that the maximum discharge value (peak) during a flood was higher than the peak value of daily simulations (Figs. 8a,b). For river depths, the RMSE was degraded by about 3–5 cm on average. The results are consistent over all the routing schemes. The results of the kinematic wave run (MS) are shown in Table 6.
a. Routing schemes
The various routing methods tested in this study were a significant improvement over the scores of the original method MD11. This result is not surprising since the MD11 parameters were fitted not only on the Garonne catchment but over the whole of France, and the optimization function was not based on the NSE but on a square error cost function (David et al. 2011a). The MCT method appeared to be inaccurate in the case of low slopes, because of a nonlinearity in diffusion. This problem has been partially solved by imposing a minimum value for the bed slopes in the calculation of the diffusion factor D* [see Eq. (19)]. Another solution could have been to change the calibration of the Manning coefficients specifically for this method, in order to have a better agreement of the discharge phase (the drawback would have been to use different Manning coefficients for the routing schemes used in this paper). In addition to these hydraulic considerations, the problem may also come from an underestimation of the runoff by the land surface scheme during low-flow periods (leading to less water in the channel and a lower velocity) and should be further investigated. The Muskingum approach is probably more suitable when computational efficiency is particularly important or when hydraulic parameters of the river channels are difficult to determine. In our case, the Manning–Strickler kinematic wave approach seems to be a good compromise as it uses a variable velocity scheme and is less sensitive than MCT to errors in slope. The mean NSE over the six downstream stations presented in Fig. 1 is increased by 0.13 in MLCC, 0.17 in MS, and 0.14 in MCT. Additional tests using only the routing part of RAPID confirmed the results obtained with the full hydrometeorological model. When forced at Tonneins with observed discharge and compared to the Marmande gauge (25 km downstream), the MS routing scheme obtained an NSE of 0.97, while the reference hydraulic model obtained 0.98. The MCT routing scheme obtained the same value, while both versions of the Muskingum scheme obtained lower scores (0.93). The evaluation of the model results at a 3-h time step shows a slight degradation of the scores, but the peak discharge is more realistic and should better compare with the instantaneous measurements of SWOT.
b. Reach averaging of river depths
River depths, which are highly variable in space, are similarly highly dependent on the values of hydraulic parameters. River depths are well simulated in river reaches where the operational rating curves are available, allowing relevant values of the roughness coefficient to be fitted. In the river reaches without operational rating curves, the simulated river depth is less accurate: as the relation between the discharge and the river depth is not known in these reaches, the relationship of Arora and Boer (1999) [see Eq. (21)] is used to determine Kstr. This usually gives good results, but may be false locally because of the high spatial variability of the river width, and to the possible bias of the average simulated discharges in the Garonne catchment.
The spatial variability of the river free-surface and riverbed elevation for MODCOU and MASCARET along the Garonne River, between the gauges of Tonneins and La Réole, is shown in Fig. 10. The riverbed altitude is more variable in MASCARET than in MODCOU, because MASCARET uses 83 observed river cross sections, while MODCOU relies on smoothed SRTM data alone. Hence, the river depth variability is underestimated in MODCOU. At the local scale (the 2-km grid cell of MODCOU), the results are erratic: considering the hydraulic models as a reference, the river depth RMSE of MODCOU is 0.53 m at Verdun-sur-Garonne (a good value by chance) and 1.95 m at Marmande. Comparisons over longer reaches lead to better and more consistent scores: 0.83 and 0.91 m, respectively, for a 10-km reach and 0.66 and 0.76 m, respectively, for a 20-km reach. It must be noted that these values are impacted by a negative bias, as the river depths of MODCOU are consistently lower than those of the hydraulic models, as can be seen in Fig. 10 for the MASCARET model: this fact can be attributed to the underestimation of the ISBA runoff in case of low discharge regimes and to the complex shape of the riverbed in MASCARET, which induces more variable flow velocities and higher river depths. Note that the bed and water elevations of Fig. 10 are not directly derived from SRTM, but from the smoothed slopes derived from SRTM and from a fit of the MODCOU water elevations on the water elevations simulated by MASCARET. An offset of 8 m has been introduced in the water elevation simulated by MODCOU in order to maintain the same water elevation in MASCARET and MODCOU in the uppermost river transect at the beginning of the simulation (1 September 1995). See section 5d for further discussion of water surface elevation.
c. Riverbed geometry
One limitation of this approach is the estimation of river widths: they are estimated as a function of the averaged simulated discharge in every grid cell of the Garonne catchment. Consequently, the river width regularly increases along the river channel. In reality, the width variability constrains the flow and influences the hydrodynamic of the flow. In the most downstream reach of the Garonne (Fig. 10), the average width computed in MODCOU is 193 m, while the estimation based on 83 observed river transects is 169 m. The standard deviation of the river width is 0.5 m for MODCOU against 8.7 m for the fine-resolution hydraulic model MASCARET. SWOT is expected to improve the estimates of river width and will complement the present efforts to estimate the widths of large rivers using satellites (Yamazaki et al. 2014).
Second, in this study we considered a rectangular geometry in every river cell of the Garonne catchment. In reality, the shape can be very heterogeneous, as shown by the high temporal variability of river widths simulated by the hydraulic models. To verify the possible impact of a modification in the shape of the bed, we tested the impact of a trapezoidal channel (not shown) on the discharge and river depths simulations at Le Mas d’Agenais. The results depended on the angle of the riverbed [as shown in Eq. (11)]. For α = 30°, the impact was very low; for α = 60°, the impact was higher. For example at Le Mas d’Agenais, the difference between the river depth simulations for a rectangular channel and the trapezoidal channel (α = 60°) was greater than 10 cm when the discharge was higher than 1500 m3 s−1 (a value higher than Q95). This means that the rectangular bed approximation is usually valid in our case, except for high angles (60° and more) and high discharge. Finally, it must be noted that because of the SRTM 90-m inaccuracies, we used smoothed bed slopes for MODCOU. Hence, the slope variability was highly underestimated in this study.
d. Water surface elevation simulations
The initial choice of this study to use the SRTM DEM in MODCOU, in order to allow an easy application of the methodology to other basin of the world, led to inaccuracies in slope estimations and water surface elevation. For example, the comparison of the simulated water surface elevation simulated by MODCOU and the hydraulic models for the 10-km reach centered on Marmande and Verdun-sur-Garonne gave an RMSE of 5.16 m for Verdun-sur-Garonne and 0.90 m for Marmande. The bias component is particularly important in the first case because of divergences in the slope variability in the area. Figure 11 also shows that the water surface elevation is less sensitive to bathymetry during high discharge periods than during low discharge periods, as found by previous studies (e.g., Trigg et al. 2009). The limitation of our approach can be explained by the following:
SRTM is representative of the elevation as of February 2000, when the mission occurred. The subtraction of a representative depth of the river may partly correct this bias.
SRTM uncertainties, combined with the fact that the MS scheme needs only downward slopes, imposed a smoothing of the riverbed leading to significant local inaccuracy of the model (e.g., in Verdun-sur-Garonne). In our case, the use of the minimum value of the cell improved the results compared to the mean value over a 2-km cell, but this choice may not be valid with a better DEM.
To compare measured and observed water surface elevation accurately, an improved DEM and additional data to account for riverbed elevation must be used. In the case of the Garonne River, national sources of data can be used, but this is not the case for other parts of the world.
Comparing water elevation (or level) variations between two SWOT observations is a potential way to lower the direct effect of water bed elevation errors on the scores. Figure 11 presents the RMSE of the river depth differences for three stations as a function of the time difference (data are not reach averaged). The RMSE increases along the river and with the time difference. However, DEM errors will indirectly influence the results by perturbing the river routing simulations.
e. Relevance for hydrometeorology
Figure 12 shows how a simulated flood propagates from upstream to downstream within the catchment. The propagation time is about 2 days in the case of a typical flood. Some river depths are poorly simulated on a limited number of river reaches because of the very simple interpolation of hydraulic parameters over the basin. As SWOT will rarely observe the Garonne catchment twice during a period of 48 h, the SWOT products must be used in conjunction with hydrometeorological models and other data in order to inform phenomena at the relevant scale in the basin. At a short time scale (1–3 days), the RMSE of the water elevation (or level) differences for three stations in the downstream Garonne strongly increases with time (Fig. 11). The errors come primarily from river routing and are progressively mixed with errors from the land surface modeling (runoff production and routing to the river) and meteorological forcing (position, chronology amount, and rate of precipitation). After 4 days, the RMSE increases at a lower rate, probably influenced by weekly to monthly and seasonal errors and biases in the land surface model and the forcing data (Fig. 11).
With four revisits irregularly distributed in time and space over a 21-day orbit cycle, further studies are needed to evaluate the phenomena that can be informed by SWOT, especially for ungauged basins. Sensitivity studies with perturbed meteorological forcing and land surface characteristics should allow the determination of the spatial and temporal scale of the hydrometeorological processes that can be observed by SWOT. These studies could benefit from an assimilation of water elevation within MODCOU, using synthetic observations with a realistic time and space distribution and realistic errors. A raw estimate of the spatial scale can be done by comparing model errors (Table 5) to the anticipated SWOT errors (10 cm over an area of 1 km2 inside the river mask, for example, for a 10-km reach of a river that is 100 m wide). Model errors at Portet-sur-Garonne and Villemur-sur-Tarn (estimated mean width 100 m, basin size 10 000 km2) are on the order of 0.28 and 0.53 m. It therefore seems reasonable to anticipate that SWOT will be relevant at this scale in basins similar to the Garonne and at a smaller scale for ungauged basins or basins with meteorological data of a lower quality.
6. Conclusions and perspectives
This study is a contribution to the building of a test bed over the Garonne basin within the framework of the SWOT mission. More precisely, the objective of this study was to evaluate the ability of a regional hydrometeorological model to simulate both discharge and river depths in the Garonne catchment over the period from 1995 to 2006.
The introduction of transfer times on some river portions significantly improved the results of David et al. (2011a). The introduction, in the model, of a flow velocity close to those observed during high floods led to an improvement of the efficiency from 0.70 in MD11 to 0.83 in MLCC for the Garonne at Le Mas d’Agenais.
The introduction of a variable flow velocity was made possible through prescribed hydraulic parameters. The hydraulic parameters were deduced from assumptions that are transferable to other comparable basins and were improved using observations where available. The three routing schemes tested further improved the results on average, especially in the downstream Garonne River.
However, this work has limitations, mainly related to the determination of the hydraulic parameters at the basin scale: riverbed elevation, width, and slope. The choice of a method based on a DEM easily available over most of the world led to some inaccurate results when compared to detailed data and hydraulic model results over two reaches of the river.
Validation with water elevation differences between two successive observations instead of absolute water elevation may be a way to reduce the bias introduced by the method. In addition, the influence of the water management on the results remains to be evaluated.
In the short term, there are two possible extensions of this work. The first is to further improve the model by using a more precise local DEM to better determine elevation and slope. The determination of width must be improved, at least over the main river. The airborne campaign AirSWOT, the projected calibration–validation and science support instrument for the SWOT mission over the Garonne, should constitute an opportunity to improve the model. With finer resolution than the SWOT mission, it will help advance the link between the finescale processes that can be simulated by hydraulic models and those at the scale of a regional hydrological model.
The second perspective is to evaluate the added value of the SWOT products for hydrometeorology in the Garonne basin and ungauged basins of similar size (especially in terms of spatial and temporal scales). This work will be achieved by establishing synthetic assimilation experiments and by perturbing the meteorological forcing and other land surface variables or parameters. The transfer of the system to other basins of the world might then be considered.
V. Häfliger is funded by Météo-France and Centre National d’Études Spatiales (CNES), as the project is partly funded by the TOSCA Programme of CNES. Cédric H. David is supported by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (NASA). The authors thank Florent Lobligeois (IRSTEA), Béatrice Vincendon (CNRM), Noël Watrin, Gérard Rauzy, and Didier Narbais (SPC/Garonne) for their help with accessing the hydrological data. Discharge observations were provided by the French Hydro database (Ministère de l’Écologie, du Développement Durable et de l’Énergie; http://www.eaufrance.fr), which gathers data from many sources. The authors also thank S. Faroux and S. Donier (CNRM) for their help in installing and using computer tools and three anonymous reviewers for their valuable comments on the manuscript.
NASA’s decision to proceed with the SWOT mission will not occur until completion of the National Environmental Policy Act (NEPA) compliance process. SWOT is a proposed NASA mission at this time and the information in this paper is predecisional, to be used for planning and discussion purposes only.