2.1 Study area, river discharge observations, and validation metrics

The simulated domain includes 12 important rivers which flow to the Mediterranean Sea (Fig. 1) and correspond to contrasting climates, from the mountainous climate leading to pluvio-nival hydrological regimes (alpine rivers) to the semi-arid climate in northern Africa or southern Turkey (Ceyhan River). These rivers contribute significant amounts of freshwater to the Mediterranean Sea. The Nile River was excluded as it is strongly affected by the association of irrigation and dam operation; thus, it is not very representative for this validation. This validation on the Mediterranean contributes directly to the HyMeX program (Drobinski et al., 2014). For each selected river, the simulated river discharge is compared to observations at the closest available station to the river mouth (Table 1). These stations have a wide range of upstream areas which vary from 2000 to 95 000 km². The corresponding observed discharge was gathered from 3 sources: (1) the daily and monthly river discharge time series from the Global Runoff Data Centre (GRDC, 56002 Koblenz, Germany); (2) the daily discharge from the French national hydrometry portal (La Banque Hydro, http://www.hydro.eaufrance.fr, last access: 3 December 2018); and (3) the daily discharge for Italian rivers (Po and Tiber) from the corresponding Italian database (Luca Brocca, CNR-IRPI, Italy, personal communication, 2016). Only seven stations provide daily discharge values. For comparison with these observations, which exhibit some data gaps (less than 15 % of the observational record but for three stations, cf. Table 1), the simulated discharge values corresponding to the missing values are eliminated. This is done to construct the hydrographs studied, as well as the derived validation metrics. Following Moriasi et al. (2007), we used the following metrics: the Pearson correlation coefficient (CC), the Nash–Sutcliffe coefficient (NS), and the normalized standard deviation (NSD, normalized by the observed standard deviation), which all indicate good performance when reaching 1; and the percent bias (PCBIAS, normalized by the observed mean discharge), the ratio of the root mean square error to the observation standard deviation (RSR), and the cross-correlation lag time (CLT), which all indicate good performance when approaching 0. The CLT is defined as the lag in days needed to achieve the maximum correlation between the simulated and observed series.

Figure 1Extended simulation domain. The main watersheds are colourized as a function of maximum upstream area (km²). They were extracted for an ORCHIDEE resolution of 1∕4^∘, with a threshold HTU number of 50. The 12 river basins studied are coloured in red. The river network is plotted in blue based on the data set from the Generic Mapping Tools (http://gmt.soest.hawaii.edu, last access: 3 December 2018). The river names are Rhône, Ebro, Po, Chelif, Maritsa, Moulouya, Ceyhan, Tiber, Adige, Shkumbinit, Devollit, and Var corresponding to the numbers from 1 to 12, respectively.

Download

Table 1Information on stations used for validation.

^* Daily data are available. Area: catchment area upstream of gauging station (km²); MQ: monthly mean discharge (m³ s⁻¹); LQ: minimum value of monthly discharge (m³ s⁻¹); HQ: maximum value of monthly discharge (m³ s⁻¹); and SD: standard deviation of monthly discharge series (m³ s⁻¹). ${}^{\mathrm{a},\mathrm{b},\mathrm{c}}$ Observation data are sourced from La Banque Hydro, the GRDC, or Luca Brocca (CNR-IRPI), respectively.

Download Print Version | Download XLSX

2.2 Numerical design

The RRS is integrated into the ORCHIDEE land surface model (Krinner et al., 2005), where the surface water budget and the resulting surface and subsurface run-off fluxes are estimated by a multi-layer soil hydrology scheme (De Rosnay et al., 2002). In this study, the model is forced by atmospheric forcing data, and three data sets are available, all based on the WATCH Forcing ERA-Interim data set (Weedon et al., 2014). They differ by the precipitation data set they are combined with to obtain bias-corrected precipitation: (1) the CRU data set (Harris et al., 2014); (2) the GPCCv5 data set (Schneider et al., 2014); and (3) the MSWEP data set (Beck et al., 2017). The first two forcing data sets are at a 1∕2^∘ spatial resolution, and the latter is at 1∕4^∘ resolution. The simulation period extends from 1979 to 2013 after a 10-year spin-up. The time step of ORCHIDEE is 30 min, whilst the time step of the RRS is 1 h, which is compatible with the abovementioned spatial resolutions based on previous theoretical experiments (not shown).

3.1 General framework

The original routing scheme in the ORCHIDEE LSM implements the linear reservoir routing method with river network building on a 1∕2^∘ global map of the main watersheds (Oki et al., 1999; Vörösmarty et al., 2000). A brief description of this scheme can be found in Ngo-Duc et al. (2007a) and Guimberteau et al. (2012a). Each river basin is constructed by connecting a number of sub-basins which are defined inside the ORCHIDEE grid boxes, following eight outflow directions. In other words, each sub-basin represents the section of the river basin within the grid box. A sub-basin can be smaller than the ORCHIDEE grid box if the ORCHIDEE model runs at a coarser resolution than 1∕2^∘ and each grid box might encompass sub-basins of different rivers. Hereafter, we will use the term “hydrological transfer unit” (HTU) to designate these sub-grid basins. Run-off from one grid box can flow into a neighbouring one or stay in the same grid box, depending on the downstream HTUs. Thus, we propose to call our scheme a “unit-to-unit routing” to distinguish it from the classical grid-to-grid routing.

To describe the transformation of run-off into river discharge along the river network, each particular HTU consists of three linear reservoirs with decreasing water residence times; these reservoirs describe the lags imposed onto groundwater flow, overland flow, and streamflow (Hagemann and Dümenil, 1997; Miller et al., 1994). The lag times are the product of a slope index (k), characterizing the HTU, and a constant specific to the type of reservoir (g), calibrated by Ngo-Duc et al. (2007a). The values of g are 0.24, 3, and 25 (day km⁻¹) for the stream, overland, and groundwater reservoirs, respectively. Following Ducharne et al. (2003), the slope index is given by $k=d/\sqrt{S}$, where d and S are the respective distance and slope between a pixel and its downstream pixel. It is first defined at the 1∕2^∘ resolution and then averaged across all 1∕2^∘ pixels composing a HTU. Surface run-off and drainage, which are computed by the soil moisture module of the ORCHIDEE model, are first lagged locally by the overland and groundwater reservoirs, then they are routed along the river network through a linear cascade of stream reservoirs. Applications of this scheme not only enhance the study of large-scale water balance (Ngo-Duc et al., 2005a, c), but they also allow for the simulation of various interactions between rivers and their watersheds (De Rosnay et al., 2003; Guimberteau et al., 2012a, b).

The availability of a high-resolution DEM (digital elevation model), namely HydroSHEDS (Lehner et al., 2008), has presented the opportunity to modify this routing scheme, with the ability to constructing more adequate HTUs in each ORCHIDEE grid box. This new RRS with a higher-resolution river graph can operate at a fine spatial resolution (e.g. finer than 10 km) and is expected to better represent the complexity of river basins and respond better to the inhomogeneity of precipitation patterns. The next sub-section briefly presents the HydroSHEDS data set, and the following sub-sections elaborate on the modifications made to the original routing scheme of ORCHIDEE to create this new version.

3.2 HydroSHEDS data

As a seamless near-global hydrological data set, HydroSHEDS is a suitable database for improving the river routing scheme in the ORCHIDEE model. It is available at resolutions from 3 to 30 arcsec, and provides all of the required information including hydrologically conditioned elevation, drainage directions, and watershed boundaries. Derived from the Shuttle Radar Topography Mission, the quality of HydroSHEDS has the limitations of this radar product, returning a complex mix of terrain elevation and vegetation height, and only covering the land areas from 56^∘ S to 60^∘ N due to the shuttle's orbit. Despite tremendous effort to void-fill and properly condition the drainage directions some errors remain, such as spurious inland sinks, whilst the assumption of single flow direction prevents river bifurcations, including deltas, from being properly described (Lehner et al., 2008). Nevertheless, this database is widely considered as the best available DEM for hydrological applications, and has shown its advantages for large-scale high-resolution river routing (Gong et al., 2011; Yamazaki et al., 2011). The stream network in HydroSHEDS is comparable with other global hydrographic data such as HYDRO1k, ArcWorld, and the Digital Chart of the World (Lehner et al., 2006).

In this study, we use a resolution of 30 arcsec (1∕120^∘ with a pixel size of 0.86 km² – ca. 1 km at the Equator). It is believed to be sufficient for large-scale applications, and offers the possibility to complete the region north of 60^∘ N with the HYDRO1k database at 1 km to achieve full global information (Lehner et al., 2008; Marthews et al., 2015; Wu et al., 2012). A preliminary quality control of HydroSHEDS was performed to ensure that all pixels drain to ocean or to a recognized endorheic basin. This was based on a non-exhaustive comparison with the global ESA-CCI land cover classification (Bontemps et al., 2013) to verify the existence of lakes or wetlands at inland outflow points. This procedure also provided an identifier for all of the river basins with an identified outlet, in particular for the small coastal basins with a zero ID in HydroSHEDS. Based on the elevation, flow directions, and basin/outlet identifier, the following information could be calculated for each HydroSHEDS pixel for further use: slope index, flow accumulation (quantifying the amount of upstream pixels), and total downstream distance to the outlet (ocean or lake).

3.3 HTU and river basin construction

As mentioned above, the major improvement of the new RRS is the possibility to better describe the geometry of the river basins by defining high-resolution HTUs beneath the ORCHIDEE grid-mesh. Figure 2a shows the major steps carried out in the river basin construction. There are two steps that control the number of HTUs in an ORCHIDEE grid-mesh which is an important issue for unit-to-unit routing with respect to the limitation of computational resources. First, the model user needs to specify the maximum area of a HTU in an ORCHIDEE grid cell which allows for the preservation of as much of the hydrological information of HydroSHEDS as required. Second, the user specifies the maximum number of HTUs. If this threshold is exceeded, in the last step of the river network construction, the number of HTUs is reduced by merging smaller units and thus increasing the average HTU's area.

The construction of the HTUs is illustrated in Fig. 2b for the Rhône River (98 000 km² in Switzerland and south-eastern France, with an outlet in the Mediterranean Sea near Marseille) and an ORCHIDEE grid-mesh resolution of 1∕8^∘. For clarity, we impose a maximum number of nine HTUs per ORCHIDEE grid box in this study, but the optimal HTU number is discussed in Sect. 4. In this framework, the number of upstream HTUs contributing water to each ORCHIDEE grid cell is depicted in Fig. 2b. Logically, the ORCHIDEE grid cells which cover the main stream of the Rhône River have a large number of upstream HTUs; therefore the main stream is denoted using dark blue, while smaller tributaries are identified in yellow. This river basin design depends on the arrangement of HTUs inside each ORCHIDEE grid box, as shown by the simple example in Fig. 2c–d. These two figures outline the definition of HTUs in the ORCHIDEE grid cell that is marked by the orange rectangle in Fig. 2b.

All of the outlets from the grid box are identified, and then all of the HydroSHEDS pixels sharing the same grid box outlet are combined to build preliminary HTUs (Fig. 2c). Following this step, a preliminary HTU can be larger than the aforementioned user-defined size (e.g. a user can set the maximum area of a HTU to 2 % of the ORCHIDEE grid box). For instance, the largest preliminary HTU in Fig. 2c covers about 80 % of the grid box area. Hence, a procedure was developed to partition these large preliminary HTUs to the user-defined size. The partitioning process relies on the Pfafstetter topological coding system for streams and basins (Verdin and Verdin, 1999): the flow accumulation is used to identify the main stream of the HTU to partition, and its four main tributaries; this results in the division of the large HTU into nine smaller HTUs which comprises the basins of the four tributaries and five inter-HTUs. All of the coloured HTUs in Fig. 2d initially belong to the preliminary HTU flowing out of the grid box at the outlet point marked by the orange circle on the eastern edge. The division process will continue if the HTU size is larger than the user-defined size, which can obviously not be smaller than the HydroSHEDS pixels (ca. 0.86 km²).

While the delineation based on the Pfafstetter codification informs on the connectivity between the HTUs deriving from the same preliminary HTU (inside one ORCHIDEE grid box), the linkage between HTUs belonging to different grid boxes requires a supplementary procedure. After a HTU is constructed it gets a unique identifier, and the corresponding HTU outlet is located (as the pixel with the biggest flow accumulation). We further calculate the total area upstream from the HTU outlet, and the distance from the HTU outlet to the river basin outlet (ocean or endorheic lake). A complex procedure uses this information to identify the downstream HTU in neighbouring grid boxes, and re-establish the connectivity and coherence of the river network. The advantage of this procedure is the possibility to function on both regular and non-regular grids. Finally, the last step not only ensures the final number of HTUs does not exceed the user-defined threshold (e.g. nine HTUs per grid box in the case of Fig. 2b) but also balances the HTU areas to avoid excessive asymmetry of the HTUs' area distribution. This is done by combining the smallest HTUs that flow out of the ORCHIDEE grid box. Note that the number of tiny HTUs increases remarkably as the ORCHIDEE resolution decreases (e.g. coarser than 1∕4^∘). For example, there are tiny preliminary HTUs which only include one HydroSHEDS pixel, as marked by the orange squares in Fig. 2d.

Figure 2A schematic diagram of catchment construction in the new river routing scheme of the ORCHIDEE land surface model (a), and an illustration of the hydrological transfer unit (HTU) idea (b–d). (b) The number of upstream HTUs contributing water to each grid box for the Rhône River at a resolution of 1∕8^∘ with a maximum number of nine HTUs per grid box (note the logarithm scale for the coloured legend). Red rivers come from the Generic Mapping Tools data set (http://gmt.soest.hawaii.edu, last access: 3 December 2018). The small orange box shows the grid box displayed in (c) and (d). (c) The orange ORCHIDEE grid box in (b) with information derived from HydroSHEDS at a 30 arcsec resolution: arrows show the flow direction of each HydroSHEDS pixel, the coloured legend shows flow accumulation, bold red arrows show the flow direction of the grid box outlets, the orange lines delineate the boundary of the preliminary HTUs, and the orange circle highlights the outlet points of the largest preliminary HTU in this grid box. (d) Partitioning of HTUs that share the same grid box outlet (marked by the orange circle). Hexagons denote the outlets of inter-HTUs based on the Pfafstetter codification, and the coloured shading indicates the different HTUs. Flow direction arrows are displayed in black and white. Orange squares highlight the approximate 1 km² HTUs.

Download

Figure 3Comparison of the modelled area (using the old and new RRS) and the reference area (from Wu et al., 2012) for the 12 river basins in this study (a). Representation of the Tiber River basin in the ORCHIDEE model with the old (b) and new (c) RRS at a regular latitude–longitude grid of 1∕2^∘, and with the new RRS at a grid of 1∕16^∘ (c). In (b–d), the colours denote the upstream area (km²) contributing water to each grid box, the orange circle is the outlet point, and blue rivers come from the Generic Mapping Tools data set (http://gmt.soest.hawaii.edu).

Download

3.4 Limitation of the old RRS

The river basin construction described above is the solution for the main limitation of the old RRS – its poor representation of small river catchments, i.e. catchments with an area under 2500 km² ($\mathrm{1}/\mathrm{2}{}^{\circ }\times \mathrm{1}/\mathrm{2}{}^{\circ }$). The reason for this shortfall is that the old RRS implements the river network with a basin description at a 1∕2^∘ resolution (Oki et al., 1999; Vörösmarty et al., 2000). Figure 3a compares the simulated areas of 12 river basins in this study to a reference area, which is from the global river network data at 1∕16^∘ spatial resolution (Wu et al., 2012). The blue triangles and green circles denote the modelled area in the new RRS with the ORCHIDEE resolutions of 1∕2 and 1∕4^∘, respectively. The modelled areas in the old RRS with the 1∕2^∘ resolution are shown by the orange diamonds. The numbers correspond to the river names given in Table 1. Significantly, the new RRS accurately constructs the river basin area of all 12 rivers (i.e. the CC values are 0.99 and 1.0). The old RRS only represents the area of the Rhône (1), the Ebro (2) and the Maritsa River (5) well, while the errors are high for other rivers such as the Po (3), the Chelif (4), and the Ceyhan River (7). In addition, the old RRS can not represent four small river basins (i.e. the Adige (9), the Shkumbinit (10), the Devollit (11), and the Var rivers (12)). Figure 3b shows the description of the Tiber River basin (in Italy, with an outlet to the Tyrrhenian Sea at Roma station) with the old RRS and a regular latitude–longitude grid at a resolution of 1∕2^∘. This river basin covers only four grid boxes with a maximum upstream area of about 9100 km², whereas the upstream area of the Roma station is about 16 545 km² according to the GRDC database. In this case, the old RRS misses the southern part of the Tiber River basin. The simulated outlet point (shown by orange circle) of the Tiber River is located about 70 km from the coastline. At the same resolution (1∕2^∘) the new RRS more satisfactorily depicts the southern part of the Tiber River (Fig. 3c). The total area is about 13 600 km² with the outlet positioned about 30 km from the coast. As the ORCHIDEE resolution increases, the representation of small rivers become more realistic using the new RRS. An example is shown on a regular latitude–longitude grid at 1∕16^∘ in Fig. 3d. The total area of the Tiber River reaches nearly 16 600 km² and the river mouth is properly placed at the coast. The position of the mouth can vary between one resolution and another because of the changing land-sea mask in ORCHIDEE. Due to the limitation of the old RRS with respect to small catchments (Fig. 3), a comparison of the simulation quality between the old and new RRS is difficult in this region. Therefore, the following sections only focus on understanding the new RRS.

3.5 Water routing process

Surface and subsurface run-off from an ORCHIDEE grid box are distributed to the overland and groundwater reservoir of the embedded HTUs in a fashion proportional to their areas. As in the original routing scheme, run-off is routed downstream with a delay time that is controlled by the number of HTUs along the stream and the properties of each HTU, namely their slope index (k) and reservoir parameter (g); the product of these two metrics defines the time lag of each HTU. The slope index is first calculated at a 1 km resolution based on the slope and length of the HydroSHEDS pixels (i.e. the aforementioned formula in Sect. 3.1), and it needs to be properly aggregated at the HTU scale. When used with the new topographic information based on HydroSHEDS, the simple averaging performed in the original version of the routing scheme leads to the consideration of water travelling a distance of 1 km, which is underestimated for most HTU sizes. We developed a new algorithm that uses the drainage directions and the resulting distance of each pixel from the HTU outlet. In practice, for each pixel, we define K as the sum of all of the 1 km values of k along the corresponding downstream line. The upscaled value of k for the HTU is then given by the product of the sum of K across all the pixels composing the HTU and the fractional area of the HTU. As a result, the slope index of the HTUs changes with the area and length of the stream lines in the HTUs (D), so that the streamflow velocity (given by $Kg/D=g/(\sqrt{S}$) does not depend, or only weakly depends, on the HTU scale. In addition, the three reservoir parameters (g) are recalibrated, leading to values of 0.01, 0.5, and 7.0 (day km⁻¹) for the stream, fast, and slow groundwater reservoirs, respectively, i.e. smaller than with the former 1∕2^∘ topographic data. These values were estimated empirically for the Rhône River basin and were implemented over the entire simulated domain. As the objective of our study is to explore the value of the new information brought by the high-resolution watershed descriptions, the parameters were determined so that, on a 1∕2^∘ grid and using the finest HTU decomposition, the routing scheme reproduces the quality of the inter-annual variability and annual cycle of the coarse resolution version. This provides a baseline against which the impact of the degradation of the HTU resolution can be evaluated on various grids.

Figure 4Sensitivity of simulation results (Nash–Sutcliffe coefficients) to different HTU sizes (HTU) in the ORCHIDEE mesh, with spatial resolutions varying from 1∕16 to 1∕2^∘. The Nash–Sutcliffe coefficient is calculated with respect to the reference case (see the text for details) at daily (a, c, e) and monthly (b, d, f) timescales. Note that the range of the Nash–Sutcliffe coefficient is different for each panel. The x axis gives the ratio of the average HTU size to the ORCHIDEE grid box area (in %). Panels (a) and (b) are Beaucaire station (the Rhône River), (c) and (d) are Pontelagoscuro station (the Po River), and (e) and (f) are Roma station (the Tiber River).

Download

4.1 Experiment design

As previously mentioned, the number of HTUs in each ORCHIDEE grid box increases the computing requirements, but it is also expected to increase the simulation quality, owing to a better description of the river flow directions and basin boundaries. Therefore, here we investigate the impact of the average size of HTUs on the simulated hydrographs, in order to better understand the numerical aspects of the new RRS, and find the best compromise between computational needs and simulation quality. To this end, the ORCHIDEE model is first used without the RRS to generate surface and subsurface run-off at a 1∕4^∘ spatial resolution, using the WFDEI_MSWEP atmospheric forcing from 1979 to 2013 (Sect. 2). These fields are then interpolated to four other horizontal resolutions of approximately 1∕16, 1∕12, 1∕8, and 1∕2^∘, using the first-order conservative remapping module of the Climate Data Operators (Schulzweida, 2018). These five sets of run-off data are then used to simulate river discharge with the new RRS. For each resolution, a reference case is chosen to preserve most of HydroSHEDS information. This would ideally be achieved with an average HTU size of 0.86 km², but we used higher values at the coarsest resolutions (1∕4 and 1∕2^∘) to limit the computing requirements. As a result, the average HTU areas of the reference cases are approximately 0.86, 0.86, 0.97, 2.60, and 8.20 km² for resolutions of 1∕16, 1∕12, 1∕8, 1∕4, and 1∕2^∘, respectively. These reference HTU areas vary between 2.0 % and 0.3 % of the grid box areas. Finally, for each resolution, several alternative RRS simulations are performed by varying the average HTU size from 0.3 % to about 80 % of the ORCHIDEE grid box areas.

4.2 Results

The analysis of these experiments focuses on three stations: Beaucaire (on the Rhône River, with an upstream area of 95 590 km²), Pontelagoscuro (on the Po River, with an upstream area of 70 091 km²), and Roma (on the Tiber River, with an upstream area of 16 545 km²).

Figure 5Same sensitivity analysis as in Fig. 4, but for the skewness of the HTUs' area distribution (a, c, e) and the total number of HTUs (b, d, f).

Download

Figure 4 shows that, for each resolution, the simulation results degrade when the ratio of the average HTU size to the area of the ORCHIDEE grid box (further abbreviated as the area-ratio) exceeds a certain value. The 1∕2^∘ simulations start to deviate remarkably from the reference case when the area-ratio exceeds 0.7 %. For resolutions of 1∕4, 1∕8, 1∕12, and 1∕16^∘, the downgrade points are at area-ratios of 1.5 %, 2.0 %, 3.0 %, and 3.0 %, respectively. Another remarkable point is that the degradation is weaker at the monthly than at the daily timescale. At the monthly timescale, the Nash–Sutcliffe coefficient with respect to the reference simulation remains above 0.7, while it drops to a negative value or near zero at the daily timescale. This highlights that the simulation results are more sensitive to the resolution and variation of the HTUs' arrangement at the daily timescale. Nevertheless, the degradation occurs at the same area-ratio for both timescales.

These area-ratio thresholds can be explained by analysing the distribution of HTUs' areas for each resolution (i.e. revealed by the skewness in Fig. 5a, c, e). For all three rivers (the Rhône, the Po, and the Tiber), the skewness peaks at these area-ratio thresholds. Skewness indicates a lack of symmetry produced by the existence of a few much larger HTUs among plenty of smaller ones. The large HTUs appear due to the combination procedure involved to control the number of HTUs. As we move to larger area-ratios, the large and more equal HTUs start to dominate, reducing the skewness of the distribution, and degrading the simulated flow due to a lack of detail in the river graphs and basin characteristics at small area-ratios.

Figure 5b, d, f shows a decrease of the total number of HTUs as the average size ratio of the HTUs increases. The number of HTUs in the case of the 1∕16^∘ resolution shows the steepest decrease. For the Rhône River, it decreases from around 120 000 to 3000 HTUs. For a small river such as the Tiber, it also drops from about 18 000 to 500. In contrast, the decreased range for the case of the 1∕2^∘ resolution is only about 11 000 and 2000 HTUs for the Rhône River and the Tiber River, respectively. The change of total catchment area and the HTUs' slope factor are also investigated, but their impacts on the behaviour of the new RRS are small (not shown). The simulated basin area only varies strongly in the case of the 1∕2^∘ resolution grid. The catchment construction with a resolution higher than 1∕4^∘ gives more steady river areas. The examination with other metrics (i.e. correlation coefficient, root mean square error, standard deviation, and kurtosis of the HTU size distribution) supports the following analysis with similar signals, thus it is not shown here.

Figure 6 highlights another aspect of the impact of HTUs' size on the simulated river discharge, by focusing on the performance of the simulations against observed discharge. With the average HTU size maintained at approximately 13 km² (corresponding to an area-ratio ranging from 0.5 % at 1∕2^∘ to 30 % at 1∕16^∘), this figure compares monthly simulated river discharge series with the observation data as the ORCHIDEE resolution increases. Therefore, this figure provides a simple investigation on the dependency of the new RRS on the ORCHIDEE resolution. The comparison is focused on two metrics, the Pearson correlation coefficient (CC) and the Nash–Sutcliffe coefficient (NS), calculated at a monthly timescale. These metrics indicate satisfactory performance, with good stability at most resolutions. The only decrease in performance is found at Beaucaire at a 1∕2^∘ resolution. As the run-off fluxes were interpolated from the 1∕4^∘ to other resolutions, this test does not include the impact of the resolution change on the other components of ORCHIDEE (water and energy budgets). If the full ORCHIDEE model was run at these resolutions, we would expect the changes in run-off and drainage to be larger than the impact of the resolution on the routing demonstrated here. The main point is the stable performance of the RRS at a fixed HTU size (i.e. 13 km²) as the ORCHIDEE resolution varies. Although the new RRS inherently depends on the ORCHIDEE resolution as it implements a unit-to-unit routing method, it can be seen that with a certain average HTU size a stable simulation quality can be expected for a wide range of ORCHIDEE resolutions.

4.3 Practical HTU size

At a given resolution of the ORCHIDEE model, a higher area-ratio (ratio of average HTU size to the ORCHIDEE grid box area) means less HTUs and thus less computational requirements. Conversely, the abovementioned analysis regarding the numerical behaviour of the new RRS indicates that the quality of the simulations deteriorates at large area-ratios, and that this degradation starts at the area-ratio corresponding to the highest skewness of the HTU size distribution. Thus, it can be assumed that the best compromise between simulation quality and computational cost is achieved at this practical HTU size, which seems rather constant across the basins studied. Here, for the 1∕2^∘ resolution, we find a practical area-ratio of 0.7 %, corresponding to a practical average HTU size of about 22 km²; for the 1∕4^∘ resolution, the practical area-ratio is 1.5 %, corresponding to a practical average HTU size of about 11 km². This allows a strong computational gain compared to RRS simulations that would be carried out at the highest possible resolution (ca. 0.86 km²), and a gain by a factor of 10 is obtained for the recommended resolution. For the domain discussed here at a resolution of 1∕2^∘ and a 33-year simulation period, the ORCHIDEE model without routing takes 5 h (wall-clock time) on a 64 core computational node of the IPSL MesoCentre (https://mesocentre.ipsl.fr/, last access: 3 December 2018). At full HTU resolution the routing adds 30 h to this simulation, while in the optimized configuration the simulation only adds 3 h. Thus, the high-resolution routing with a practical area-ratio of 0.7 % only increases computing time for ORCHIDEE by 60 %.

Table 2Evaluation metrics for monthly river discharge simulations by the ORCHIDEE model with the new RRS. NS: Nash–Sutcliffe efficiency; PCBIAS: percent bias; and RSR: ratio of root mean square error with observation standard deviation.

(1), (2), and (3) are simulation results with forcing data from WFDEI_CRU, WFDEI_GPCC, and WFDEI_MSWEP, respectively.

Download Print Version | Download XLSX

Figure 6Performance of the new RRS as the resolution of the ORCHIDEE mesh changes from 1∕2 to 1∕16^∘. All simulations are performed with the same average HTU area of about 13 km². RRS performance is assessed at a monthly timescale against observed monthly time series at Beaucaire (solid lines) and Pontelagoscuro (dashed lines) stations. Red lines denote the correlation coefficient and blue lines represent the Nash–Sutcliffe coefficient.

Download

5.1 Experiment design

The goal of this section is to evaluate the performance of full ORCHIDEE simulations against river discharge observations. To this end, the new RRS is fully coupled to the ORCHIDEE LSM, which is driven by the three atmospheric forcing data sets presented in Sect. 2: the WFDEI_CRU and WFDEI_GPCC with a spatial resolution of 1∕2^∘ and the WFDEI_MSWEP with a spatial resolution of 1∕4^∘. The three corresponding simulations are performed with the practical average HTU sizes identified above, viz. 22 km² at 1∕2^∘ and 11 km² at 1∕4^∘. For comparison, similar simulations are also performed with different average HTU sizes (from 0.3 % to 50 % of the grid box area), and also with the old RRS. In the latter case, the ORCHIDEE simulations are only run with the two 1∕2^∘ forcing data sets, owing to the limitations of the OLD RRS at finer resolutions. All the simulations extend from 1979 to 2013 after a 10-year spin-up.

5.2 Monthly timescale

Figure 7 shows that the new RRS satisfactorily captures the seasonal cycle of observed discharge at the Beaucaire (the Rhône river) and Pontelagoscuro (the Po river) stations. For Pontelagoscuro, the simulations adequately reproduce two characteristic peaks of the pluvio-nival regime, in October and May. For Beaucaire, the new RRS not only captures the high peak in January but also reproduces the gradual decrease to low flow in August, except with the WFDEI_CRU data. Most simulations also display a positive bias, which can be attributed to excessive run-off fluxes (as the RRS is conservative and does not change the long-term mean discharge) and explained by systematic errors in the water budget parameterization of ORCHIDEE or biases in the three forcing data sets. The purple box plots present the sensitivity of simulated discharge to the average HTU size (from about 0.3 % to about 50 % of the grid box area) and quantify the uncertainty coming from the numerical choices of the scheme. As the HTUs average size varies, the simulated discharge fluctuates in a range which is much smaller (about 40 % smaller) than the magnitude of the bias when compared to observation. In other words, it can be said that the numerical uncertainty is small compared to the uncertainty in the forcing data. According to Beck et al. (2017), the quality of the precipitation data from MSWEP is generally better than WFDEI_CRU when compared with observations from 125 FLUXNET stations. This is confirmed by the lower magnitude of the bias in the case which uses WFDEI_MSWEP (Fig. 7e, f), while this error, when using WFDEI_CRU, is larger at both stations (Fig. 7a, b). Thus, we can confirm that accurate atmospheric forcing is an important factor which determines the performance of a RRS, as already reported by many studies (Guimberteau et al., 2012a; Ngo-Duc et al., 2005c; Pappenberger et al., 2010). In addition, the simulation results using the old RRS (shown by the orange lines in Fig. 7) confirm that the new RRS better matches the observed discharge, which is likely related to its better estimation of the river basin area and structure.

Figure 7Mean annual cycle of river discharge at the Beaucaire (a, c, e) and Pontelagoscuro (b, d, f) stations simulated with three different forcing data sets (WFDEI_CRU: a, b, WFDEI_GPCC: c, d, and WFDEI_MSWEP: e, f). The solid black line shows the observed discharge, and the orange line shows the simulation results with the old RRS. The other lines show simulation results with the new RRS and the practical average HTU size (see text for detail). The small coloured ticks along the y axis give the average values, and the purple box plots show the range of results from simulations with the new RRS and different average HTU sizes (the limits of the boxes correspond to the first and third quartiles).

Download

Validation at the 12 stations highlights that the new RRS simulates acceptable river discharge at a monthly timescale (Fig. 8 and Table 2). The black star shown in the Taylor diagram (Fig. 8) indicates where simulated results would have the same amount of variation and perfect linear correlation with observations. For this reason, simulations that agree well with the observations will lie close to this reference point. Over the 12 stations, the CC are all above 0.6 and the NSDs are in a range from 0.5 to 1.5. The new RRS achieves the best performance at the Beaucaire (1), Pontelagoscuro (3), Yakapinar Misis (7), Kokel (11), Paper (10), Malaussene (12), and Tortosa (2) stations. It is interesting that these stations display a wide range of upstream areas, with a monthly mean discharge from about 30 to 1500 m³ s⁻¹ (Table 1). The worst overall results are found at the Boara Pisani (9), Dar El Caid (6), Roma (8), Sidi Belattar (4), and Meric Koep (5) stations. For Boara Pisani (Adige River, northern Italy), the NSD values from the WFDEI_CRU and WFDEI_MSWEP experiments are both higher than 2.2, which is linked to the high positive bias with these forcing data sets (i.e. the estimated discharge is about twice as high as observations, as shown in Table 2). As the overestimation in the WFDEI_GPCC experiment is smaller (i.e. the PCBIAS is about 22 %), the NSD value stays around 1.10. However, it should be noted that the CC values for this station range from 0.77 to 0.94, so the monthly variability is quite well captured. In particular, the simulated monthly series reproduce the observed hydrological regime well with rather weak seasonal variations and two moderate peaks (in June and October); poor performance is generally found in the driest part of the Mediterranean Basin. The negative bias at Roma station can be attributed to the underestimated run-off during summertime (July to September) or water management practices. The underestimation of discharge is even worse at Sidi Belattar station, which belongs to the longest river in Algeria and rises from the Saharan Atlas. The observed annual mean discharge in this area is lower than 20 m³ s⁻¹, with very low values in summertime (close to zero); none of the simulations can reproduce this characteristic, regardless of the forcing data set utilized. The time series of the anomaly of monthly discharge (with respect to the mean seasonal cycle) are also analysed in Fig. 8b to assess the inter-annual variability of discharge. It shows about the same perspective as the monthly series with lower CC values, although good CC (about 0.9) are still found at Beaucaire (1), Pontelagoscuro (3), and Yakapinar Misis (7) stations. The errors at the Boara Pisani (9), Sidi Belattar (4), and Dar El Caid (6) stations are more clearly demonstrated. Regarding the effect of different forcing data, no clear hierarchy is found, although the above findings confirm that the biases of surface and subsurface run-off by the ORCHIDEE LSM, which are at least partly due to the biases of the forcing data sets, have a strong impact on the simulated streamflow. This is probably the main reason for the low NS values and the high absolute PCBIAS values (shown in Table 2).

Table 3Evaluation metrics for daily river discharge simulations by the ORCHIDEE model using the new RRS: CC – Pearson correlation coefficient; NS – Nash–Sutcliffe efficiency; NSD – normalized standard deviation; CLT – cross-correlation lag time; and RSR – ratio of root mean square error with observation standard deviation.

(1), (2), and (3) are simulation results with forcing data from WFDEI_CRU, WFDEI_GPCC, and WFDEI_MSWEP, respectively.

Download Print Version | Download XLSX

Figure 8Taylor diagrams comparing simulated and observed river discharge at the 12 stations for (a) monthly series and (b) monthly anomalies with respect to the mean seasonal cycle. Blue triangles, red squares, and green circles correspond to simulations with forcing data from WFDEI_CRU, WFDEI_GPCC, and WFDEI_MSWEP, respectively. The number of each station is given in the figure legend and in Table 1. The black star shows the observed value at which the normalized standard deviation and correlation coefficient are 1.0.

Download

It is worthy to note that the impact of human regulation on natural streamflow is neglected in the current version of the new RRS. As a result, irrigation is not represented at all in these simulations, although it is known to play an important role in the Mediterranean region (Margat and Treyer, 2004). According to Montanari (2012), the annual water withdrawal for irrigation in the Po River basin is about 17 km³, i.e. a third of the mean annual discharge (47 km³). Not all withdrawals are transformed into evaporation; thus, a part of these abstractions will return to the river. Nevertheless, one can assume that the observations for Pontelagoscuro (displayed in Fig. 7) probably underestimate the natural river discharge. Snoussi et al. (2002) also underline that the water discharge from the Moulouya River (Dar El Caid station) has been reduced by almost 50 % due to the construction of the Mohamed-V reservoir (in 1967), which could explain the large positive biases of our simulations. Figure 9a shows that, in the Ebro Basin (Spain), the discharge peaks in May and June are largely overestimated by the simulations, regardless of the forcing data set used. However, this is not the case when referring to the observed record from 1920 to 1930. During this period, human impacts on the natural Ebro river flow were not significant and the discharge peak from May to June was stronger than the peak associated with the winter rains. This analysis strongly supports the fact that the overestimation of river discharge over the last decades could be alleviated by a proper representation of human water abstractions for irrigation in the new RRS. For the Maritsa River, Artinyan et al. (2007) noted that about 7 % of the river flow is used for irrigation during summertime (June to August). Despite the fact that this study only covers the year 1996, their findings nevertheless underline the role of anthropogenic influence on river processes, which probably contribute to the positive bias of our simulations. However, in contrast to the Ebro, the spread in the river discharge bias caused by the forcing data uncertainty is larger than the bias itself, suggesting that the role of human processes is not as large in the Maritsa River as it is for the Spanish basin.

Figure 9Mean annual cycle of river discharge at the Tortosa station on the Ebro River (a) and Meric Koep station on the Maritsa river (b) simulated with forcing data from WFDEI_CRU (blue), WFDEI_GPCC (red), and WFDEI_MSWEP (green). The solid black line is the observation data after 1979, while the orange line is the observation data for the period from 1920 to 1930 (at the Tortosa station). The small coloured ticks along the y axis give the average values.

Download

5.3 Daily timescale

Simulation quality at a daily timescale is validated at seven stations, listed in Table 1. Table 3 presents five metrics considering the variability of daily time series (i.e. CC, NS, NSD, and CLT) and the magnitude of the error (i.e. RSR). The daily CC values are slightly smaller than the monthly CC values, but remain higher than 0.6 for five stations, where the short-term variability of streamflow is correctly reproduced. The best result can be found for the Pontelagoscuro station with forcing data from WFDEI_GPCC (Table 3). The two exceptions are Roma station (Tiber) and Sidi Belattar station (Chelif), which displayed weak CC at the monthly timescale. Table 3 also presents CLT values in a range from −25 days to 1 day, which suggests that some CC values can be improved by accelerating the transfer of water in the RSS. We also find that the daily NS values are much weaker than the monthly values. At the daily timescale, acceptable NS values, which are usually taken as higher than 0.5 (Moriasi et al., 2007; Tavakoly et al., 2017), are only found at the Pontelagoscuro (Po) and Yakapinar Misis (Ceyhan) stations, which are the stations maximizing the monthly NS. The reason for this is that the NS is very sensitive to mass balance errors, i.e. river discharge bias (Gupta and Kling, 2011; Krause et al., 2005), and to daily time series of small dry rivers (Schaefli and Gupta, 2007). Although the NS is very commonly used to evaluate hydrological simulations, this criteria still raises many questions regarding its application (Criss and Winston, 2008; Gupta et al., 2009). The largest biases also explain why all RSR values in Table 3 are larger than 0.5. The large fluctuations of all the metrics between the simulations using the three different forcings again highlight the importance of the accuracy of the atmospheric conditions imposed on the LSM.

For the four stations with the best daily NS values, Fig. 10 visualizes the simulation results at the daily time scale using flow duration curves, which represent the percentage of time over which discharge exceeds the discharge value indicated on the vertical axis. This analysis shows that discharge at Beaucaire station (Fig. 10a) is overestimated over the full range of frequencies, except with WFDEI_MSWEP, which induces lower flows than observed for the 20 % lowest flows (exceedance frequency above 80 %). Except when using WFDEI_CRU forcing data, realistic statistics are found at almost all frequencies at Pontelagoscuro station (Fig. 10b), explaining the good results presented in Table 3. In contrast, at Yakapinar Misis station (Fig. 10c) the RRS captures the low flows under 500 m³ s⁻¹ well, but the highest flows which occur less than 10 % of time are strongly underestimated. The poor performance at the Sidi Belattar station is confirmed by Fig. 10d. At this station, streamflow exceeding 100 m³ s⁻¹ only occurs 5 % of time, whilst more than 50 % of the time the flow is close to zero. It is difficult to capture the daily discharge at this station, and the large agricultural water demand over the Chelif Basin (Mahe et al., 2013) probably contributes to the errors. However, the main error source seems to be the forcing data sets, as the observations fall in their wide uncertainty range.

Streamflow fluctuations possess alternative frequencies in addition to the daily and monthly oscillation. This can be seen in the power spectrum of daily river discharge at Beaucaire in Fig. 11a. The power spectrum corresponds to the squared amplitude at each frequency and was extracted using discrete Fourier transform spectroscopy and was smoothed with a Savitzky–Golay filter (Savitzky and Golay, 1964). It is a robust method for analysing the multi-scale temporal variations in the various variables contributing to the river discharge. As has been shown by Weedon et al. (2015), it facilitates the diagnosis and development of hydrological models. In the power-spectra trend lines for high frequencies (variation faster than 30 days) and low frequencies (slower than 30 days) are also plotted.

Based on the previous diagnostics in Sect. 5.2 and 5.3 we expect the new RRS to replicate the variation of river flow at various frequencies well. Figure 11a shows the good match of the two power spectrum patterns. The β values for the simulation are 0.91 and 2.39 at frequencies above and below 30 days, respectively, very close to the corresponding slopes for the observations, i.e. 0.81 and 2.28. At high frequency, there is a mismatch of the power density around the 3-day period. It can be traced back to the simulated subsurface run-off, which shows a similar peak during the same period (Fig. 11b). The peak at 3 days is probably characteristic of the soil moisture diffusion scheme of the ORCHIDEE model, as it does not have a signature in precipitation and evaporation averaged over the entire Rhône catchment (Fig. 11). The spectra of these two fluxes which characterize the exchange with the atmosphere are more noisy at the synoptic scales (below 15 days) and display a stronger slope for low periods. These results are indicative of soil moisture diffusion in ORCHIDEE that is too fast or of an insufficient buffering of the resulting subsurface flow (drainage) by the routing reservoir representing groundwater flow. It also further suggests that the link between the soil hydrology and the routing scheme is too simple in ORCHIDEE and lacks an appropriate representation of the aquifers. Human activities and the regulation of river flows also affect the river discharge variability at many frequencies (e.g. hydropeaking as a result of the storage for hydropower plants; Meile et al., 2011); thus, the validation of this variability either requires the analysis of pristine catchments or the representation of water management infrastructures in the model.