Modelling river runoff and estimating its weather-related uncertainty for 11 large-scale rivers located in different regions of the globe

Simulating the dynamics of water balance components of river basins accounting for the processes of their formation is one of the main issues of modern hydrologic science. Observed and projected global climate change are motivating the scientific community towards the solution of this problem, forcing improvement of the existing methods and to look for new approaches for assessment and projection of river runoff and the other water balance components (Peel & Blöschl 2011). In spite of the fact that during the last 40 years there have been many attempts to create new hydrological models, which should be unique and universal, it can be stated that this problem has not been solved at the moment (Singh & Woolhiser 2000; Beven 2011).

This is also confirmed by the appearance of various international projects devoted to comparison of different hydrological models. One such project is the Inter-Sectoral Impact Model Intercomparison Project (ISI-MIP), launched in 2013 (Warszawski et al. 2013), the water sector of which is aimed at validation and improvement of hydrological models on regional and global scales and application of the models for hydrological projections of possible changes in water resources due to expected climate change. The authors of the present paper are participating in this project with their Soil Water–Atmosphere–Plants (SWAP) model (Gusev & Nasonova 2000; Gusev et al. 2006) and the results presented here were partly obtained within the framework of the second phase of the ISI-MIP (ISI-MIP2).

The SWAP model has been widely validated during the last two decades. As a rule, the validation was carried out using meteorological forcing data derived from observations at meteorological stations. However, such observations are not always available, do not cover the whole land surface of the Earth and represent point data. At the same time, there are many global data sets containing meteorological data at different temporal and spatial resolutions (e.g. Xu et al. 2016). That is why the first goal of the present work was investigating the ability of SWAP to reproduce runoff of rivers located in different regions of the globe under a wide variety of natural conditions (from areas with permafrost and severe arctic climate to humid areas of equatorial monsoons and arid regions of tropical and subtropical latitudes) using meteorological forcing data and land surface parameters derived from the global data sets suggested in the framework of the ISI-MIP2 project.

Similar large-scale multi-basin modelling studies were performed by nine hydrological models participating in the ISI-MIP. The summarized results can be found in the overview paper (Krysanova & Hattermann 2017).

Another important problem associated with hydrological simulation and projections is their inevitable uncertainties resulting from application of different climate and hydrological models, different greenhouse-gas emission scenarios, model parameters and so on. There are many publications devoted to estimation of these uncertainties, which should be taken into account at the stage of decision-making in water resource management and planning. Some of these uncertainties (epistemic ones) can be reduced, for example, due to obtaining new knowledge about climatic and hydrological systems, and improving the quality of models and data. However, there are also structural uncertainties which do not depend on our knowledge or data and are an inherent property of these systems. Such uncertainties cannot be reduced and should be treated as a lower limit of predictability in climatic or hydrological modelling (Gelfan et al. 2015). Considering river runoff, structural uncertainty of its simulations and projections is associated, in particular, with natural variability of runoff.

This variability, as shown in Gelfan et al. (2015) and Gusev et al. (2017), is caused by the chaotic nature of characteristics of atmospheric processes (so-called weather noise), their instability with respect to small changes in the initial conditions. The latter is confirmed by numerical experiments with atmosphere-ocean general circulation models (AOGCMs) when small errors in the initial conditions resulted in large scatter in simulations of the evolution of meteorological variables (see, for example, Gelfan et al. (2015)). In nature, differences in the initial state of the atmosphere (for example, at the beginning of each year), result from different (even very small) positions of the Earth on its path around the Sun, as well as the difference in the state of the atmosphere which has been forming during a previous period. Small changes in the initial state of the atmosphere (a nonlinear system with dissipation) lead to quite large changes in its further dynamics. This is one of the reasons that it is impossible to predict an actual dynamics of the climate system, because the initial conditions of the state of the atmosphere can never be prescribed with absolute accuracy.

The weather noise of atmospheric processes causes a corresponding noise of hydrological processes, resulting in natural uncertainty of hydrological variables (in particular, river runoff). This uncertainty can be treated as the lowest limit of predictability of river runoff, which means that further improvement of climatic, hydrological and land surface models (LSMs) cannot help one to obtain long-term hydrological projections with a lower uncertainty than this limit. Such estimates are of great interest because they provide an opportunity to compare changes in runoff due to climate change or man-made impacts with natural variability of river runoff (i.e., its variability in the case of stationary climate). In this context, the second goal of the present work was estimating and analysing the natural variability of river runoff caused by the weather noise.

Eleven large-scale river basins selected in the ISI-MIP were used in the present study. They are situated in different geographical zones in different continents (Figure 1) and include the Rhine and Tagus located in Europe, the Lena, Ganges, Upper Huanghe (Yellow), and Upper Yangtze in Asia, the Niger in Africa, the MacKenzie and Upper Mississippi in North America, the Upper Amazon in South America and the Darling in Australia. The areas of the basins range from 67,490 km² to 2,460,000 km² (Table 1). A description of the basins can be found in Krysanova & Hattermann (2017). Here, it should be noted that the basins substantially differ from one another by natural conditions (climatic, soil, vegetation, geomorphological, etc.) which is very important for the solution of the stated problems.

A spatial schematization of the river basins was derived from Wu et al. (2012). The basins are presented by a set of regular grid cells (with a spatial resolution 0.5° × 0.5° in latitude and longitude) connected by a river network. The schemes of the selected river basins are shown in Figure A1 (available with the online version of this paper) and the number of computational grid cells for each basin is given in Table 1.

Meteorological forcing data needed to drive the SWAP model are the air temperature and humidity, precipitation, incoming fluxes of longwave and shortwave radiation, wind speed and air pressure. These data for each calculation grid cell of the selected basins for the period of 1970–2001 were derived from the Water and Global Change (WATCH) data set with 0.5° × 0.5° spatial resolution (Weedon et al. 2011). WATCH was created on the basis of combining monthly values of ground-based measurements (including air temperature, the number of days with precipitation, cloudiness and precipitation) taken from the Global Precipitation Climatology Center (GPCC) and the Climatic Research Unit (CRU of University of East England) data sets, and the ERA-40 reanalysis data (with corrections to seasonal changes of atmospheric aerosols improving incoming solar radiation).

Most vegetation parameters for the selected river basins were derived on the basis of the global ECOCLIMAP data set which has 1-km resolution and is based on classification including 215 ecosystems (Champeaux et al. 2005). The obtained parameters were aggregated by means of mosaic (tile) approach for each calculational grid cell.

As well, SWAP needs the following soil hydrophysical parameters: porosity, field capacity, wilting point, hydraulic conductivity at saturation, as well as soil matric potential at saturation and B-parameter in parameterization by Clapp & Hornberger (1978). These parameters were calculated by making use of equations given in Cosby et al. (1984) and data on soil texture (CLAY and SAND) taken from the ECOCLIMAP.

The obtained parameters were checked for consistency and some unrealistic values were corrected as proposed in Gusev et al. (2006).

Streamflow data for the gauges listed in Table 1 were provided with a daily time step (with the exception of the Ganges, for which only monthly data were available) within the framework of the ISI-MIP2 project (see Huang et al. (2017) for details). Additional daily runoff data for subbasins of the Lena and MacKenzie rivers were taken from the the Global Runoff Data Centre (GRDC).

The LSM SWAP is a physically based model describing heat and water exchange processes in a groundwater–soil water–vegetation/snow cover–near-surface atmosphere system (SVAS). The model treats the following processes: interception of liquid and solid precipitation by vegetation; evaporation, melting and freezing of intercepted precipitation, including refreezing of melt water; formation of snow cover at the forest floor and at the open site during the cold season; partitioning of non-intercepted precipitation or water yield of snow cover between surface runoff and infiltration into the soil; formation of the water balance of the aeration zone including transpiration, soil evaporation, water exchange with lower layers and dynamics of soil water storage; water table dynamics; formation of the heat balance and thermal regime of SVAS; soil freezing and thawing.

The SWAP model makes it possible to simulate the dynamics of heat and water balance components for terrestrial ecosystems and river basins on different scales and under different natural conditions. Spatial heterogeneity of the territory is taken into account by application of the function of distribution of hydraulic conductivity at saturation (Gusev & Nasonova 1998). Different versions of SWAP and their validation have been detailed in a number of publications (e.g., Gusev & Nasonova 2000; Gusev et al. 2008, 2016a; Nasonova et al. 2009).

In hydrological applications, in the case of a grid box or a small catchment, a kinematic wave equation is used to simulate streamflow at the box/catchment outlet. Application of the SWAP model for a large-scale river basin requires schematization of the basin as a number of computational grid cells connected by a river network. Heat and water regimes are simulated independently for each grid cell. Runoff simulated for each cell is then transformed by a river routing model to obtain streamflow at the basin outlet. Here, a simple linear transfer model in river channels (Oki & Sud 1998) is used.

Participation of the SWAP model in numerous international projects, such as PILPS (phases 2a, 2c, 2d and 2e), Rhona-AGGregation, MOPEX, ALMIP-2 devoted to the comparative analysis of the quality of existing LSMs and hydrological models has shown that SWAP is able to reproduce the dynamics of water balance components and other hydrothermal characteristics of river basins quite adequately provided that input data of high quality are available (e.g., Wood et al. 1998; Gusev & Nasonova 2000; Boone et al. 2004; Nasonova et al. 2009; Getirana et al. 2017; Grippa et al. 2017).

If NS = 1, the simulation is ideal. If NS < 0, temporal variability of variable х is reproduced poorly (in this case, a simple averaging of observations is better than model simulation). As to the Bias, it should be taken into account that it is desirable to reduce the error of runoff measurements to ∼5% (WMO 1994) (this value may be much higher for flood periods). As such, the quality of runoff modelling can be considered as ‘good’ when |Bias| ≤5–10%.

Parameters' calibration was carried out automatically using the algorithm for searching for the global extremum SCE-UA (Duan et al. 1992). The Nash and Sutcliffe efficiency was used as an objective function. The search for the maximum of the objective function was performed, as it was suggested in Nasonova et al. (2009), under the condition that an absolute value of the systematic error of simulation does not exceed 5%.

For the snow-dominated river basins (here, the Lena and MacKenzie), seven model parameters were calibrated: hydraulic conductivity at saturation, soil column thickness, rooting depth, albedo of vegetation in the warm period, albedo of snow cover, Manning's roughness coefficient, and the effective velocity of water movement in the channels (Gusev et al. 2008). For the rain-dominated rivers, the correcting factor for mean monthly leaf area index LAI was calibrated instead of albedo.

As previously mentioned, the forcing data were derived from the WATCH based on the reanalysis ERA-40 which may contain systematic errors as any product of reanalysis. To reduce the errors, the procedure of hybridization with monthly observations is often applied. In the WATCH, monthly measurements of precipitation were used for bias correction, while another important variable influencing runoff formation − incoming radiation − was not corrected. That is why in our study, correcting factors to incoming shortwave and longwave radiation were applied and calibrated (together with the model parameters) for all river basins (as in Gusev et al. 2008; Nasonova et al. 2009).

As well, to improve the results of calibration, the Lena and MacKenzie river basins were divided into three parts because of their large areas and optimal values of parameters were obtained for each part. As shown in Gusev et al. (2016a, 2016b), such a division allows the improvement of the quality of simulations.

Calibration was performed for eight years (different for different rivers depending on available measured runoff); the rest of the years were used for model validation with the exception of the Tagus and Ganges rivers for which all available measured runoff (during eight and four years, respectively) was used for calibration.

In order to analyse the differences in runoff uncertainties obtained for different rivers, the spectral analysis which is widely used both in technical fields and in atmospheric sciences (Coley & Tukey 1965; Panofsky & Brier 1968; Belotserkovskii 1993) was applied. In hydrology, the spectral decomposition of the water balance components (precipitation, evapotranspiration and river runoff) was performed and analysed in Gusev et al. (2017).

Since natural uncertainty of river runoff is associated with its stochastic component (weather noise), one can try to determine, as the first approximation, the ratio of a regular component of runoff to a stochastic one. For this purpose, the expansion of a temporal evolution of river runoff into the Fourier series can be used. Since any variable (both deterministic and random) can be represented as the sum of harmonic oscillations (Fourier series), this procedure was carried out for an annual course of runoff of the rivers under study and spectral densities S of daily river runoff were calculated for harmonics with different frequencies f. For this purpose, the measured values of daily runoff were used for all rivers except for the Ganges River, for which daily measurements were not available. In this case, the simulated daily runoff was used.

Figure 2 shows a typical spectral density S of river runoff for harmonics with different frequencies f (Gusev et al. 2017).

As can be seen in Figure 2, the spectral density of river runoff has a clearly pronounced maximum at the frequency f_аn = 2.738 × 10⁻³ day⁻¹ (marked with the number 1), which corresponds to a temporal period of one year. This situation is obvious due to the rotation of the Earth around the Sun with an annual period, and the component S(f_an) as well as components which are close to it in frequency (as Hamming window was used for the calculation of S ‘to involve a trade-off between resolving comparable strength components with similar frequencies and resolving disparate strength components with dissimilar frequencies’ (Harris 1978)) are regular deterministic components of S(f). As well, there may be other peaks corresponding to periods T/n, where T = 1 year, n = 2, 3, … (marked with the numbers 2 and 3 in Figure 2). In particular, the second peak corresponds to a semiannual period and so on. All these peaks are caused by deterministic reasons, for example, by influence of spring flood, monsoon rain period, etc. However, the most important thing is that the function S(f) has a maximum f_an. The value of S(f_an) was chosen as an indicator of the contribution of deterministic components into the spectral density S(f).

Random fluctuations of river runoff caused by the weather noise also contribute to fluctuations of the spectral functions. An index reflecting the contribution of random components of runoff was created as follows. The largest peaks (usually 1–3, see Figure 2) of S(f) were ‘cut off’ to remove the most significant part of the regular component. It should be noted that the number of peaks can be increased; however, it was found that cutting off further peaks, as a rule, has little effect on the final result. Then, for the transformed function S(f), the standard deviation σ_s caused primarily by the weather noise was calculated. Note that values of S(f) at frequencies larger than 0.03–0.1 (when S(f) sharply decreases) do not make a significant contribution to σ_s-values.

It should be noted that the suggested criterion reflects this relationship approximately, because it is based on the data of only one random realization (occurred in nature) of the dynamics of river runoff (out of a large number of potentially possible realizations). Nevertheless, it can be used as the first approximation that is confirmed by the obtained good agreement between the uncertainty of river runoff and the proposed criterion (see below).

At the first stage SWAP was run using meteorological forcing data and land surface parameters derived, respectively, from the WATCH and ECOCLIMAP data sets and called a priori input data. Simulations were carried out with a daily time step for the period of 1970–2001. The simulated runoff was compared with measurements for different periods depending on the availability of measurements.

Application of a priori model parameters resulted in a poor agreement between measured and simulated monthly river runoff: for most basins NS < 0, Bias ranged from 3% (for the Upper Mississippi) to 340% (for the Darling) (Table 2).

Calibration of the parameters resulted in substantial improvement of model performance. Figure 3 illustrates the ability of SWAP to reproduce monthly hydrographs of rivers located under contrasting natural conditions: for the Upper Amazon (with the largest runoff R and runoff ratio R/P equalled to 1,460 mm/year and 0.69, respectively), Darling (with the lowest R and R/P equalled, respectively, to 8 mm/year and 0.01) and MacKenzie (with the negative mean annual air temperature and R/P = 0.39, the latter is close to the average value among the 11 rivers).

Figure 4 provides mean monthly hydrographs for all the rivers simulated by SWAP with a priori and optimized parameters. Comparison of both model runs with measured runoff shows that calibration of model parameters substantially improved the quality of simulations. This is also confirmed by goodness-of-fit statistics given in Table 2 for the calibration and the whole simulation periods. As can be seen, for the calibration period, absolute value of Bias does not exceed 5.1% (mean value is 2.4%), NS coefficient ranges from 0.68 to 0.95 (mean value is 0.85). For the whole period, absolute Bias is lower than 13.8% with a mean value of 4.4%, NS-values are within the interval of 0.58–0.90 with a mean value of 0.78. The obtained results can be treated as good for all the rivers with the exception of the Lena, for which the quality of runoff simulations is the lowest and can be treated only as satisfactory. This can be explained by the lower quality of meteorological forcing data for the Lena River basin. In Gusev et al. (2016a), where forcing data were derived from measurements at meteorological stations located within the Lena basin, the quality of simulation of monthly river runoff was much higher in terms of NS (NS = 0.93), in spite of a coarser spatial resolution (1° × 1°).

Nevertheless, in general, for all the rivers the specific features of river runoff are reproduced quite adequately, although the rivers are situated in different natural conditions with respect to both climate (from arctic to equatorial and from arid to humid) and soil/vegetation (from unfrozen and seasonally frozen soils to permafrost, as well as from tundra to subtropical forests).

Figure 5 shows averaged over the simulation period annual runoff for the 11 rivers and the boundaries of its natural variation (absolute uncertainty Un_abs), as well as annual values of measured runoff for the same period. The values of uncertainty of annual runoff were estimated under the assumption that the probability distribution of annual runoff can be approximated by a log-normal distribution function. As can be seen in Figure 5, the measurements are within the range of uncertainty with 95% probability for all rivers. This is also confirmed by the results obtained for monthly runoff. Figure 6 shows the modelled hydrographs averaged over the simulation periods, the ranges of their absolute uncertainties and monthly runoff measured in different years.

Relative uncertainties of annual and monthly runoff are given in Table 2. As can be seen, monthly uncertainties for all rivers are higher than the annual ones. The reason is that for the greater averaging intervals, the multidirectional deviations of a variable from its mean value are mutually compensated, resulting in a decrease in the relative scatter of the values.

As well, analysis of the obtained results shows that uncertainties of river runoff largely depend on natural conditions of a river basin. In particular, they depend on the degree of determinism of the annual course of river runoff. The lowest uncertainty of annual runoff corresponds to the rivers with a clearly pronounced annual course of their runoff. Thus, the northern rivers (the Lena and Mackenzie) with a snow-dominated feeding and clearly pronounced spring flood have relatively low values of runoff uncertainty. The Upper Amazon is also characterized by a relatively small uncertainty because snow feeding also plays an important role in this part of the river basin. In addition, some smoothing of runoff hydrograph of the Upper Amazon occurs due to the fact that its right tributaries are located in the southern hemisphere, while the left tributaries are in the northern one. The Niger, Ganges, Yangtze and Huanghe rivers with more or less pronounced periodical impact of summer monsoon have slightly higher uncertainty of runoff compared to the rivers with a distinct spring flood peak.

Further increase in runoff uncertainty was found for the Rhine and Upper Mississippi rivers which are characterized by a mixed feeding (both snow and rain). The annual course of precipitation is smoothed in the Mississippi basin and absent in the Rhine basin (see Figure 2 in Krysanova & Hattermann (2017)). In addition, these rivers (especially the Rhine) do not have a clearly pronounced spring flood every year. This means that the contribution of regular components to runoff formation is lower compared to the above seven rivers.

Finally, the Tagus and Darling, being the rain-dominated rivers, have the highest uncertainty of river runoff. The climate in the Tagus basin is Mediterranean with a strong continental influence. Not only monthly, but also annual precipitation varies greatly in space and time. As to the Darling, the climate in the basin can be treated as tropical desert, and the river is supplied with water very irregularly. The Darling's numerous tributaries have their sources in the arid part of Australia, and many of them reach the main river only during occasional periods of flood. In the dry season, the river becomes much shallower or dries up. Nevertheless, in summer (although rarely), the river can cause damaging floods due to the climate and orography of the area. This results in the highest uncertainty in the Darling's runoff (Table 2).

The above considered qualitative analysis of the relationship between the uncertainty of river runoff and the ratio of its regular and stochastic components can be illustrated quantitatively using the above described technique for spectral decomposition of the dynamics of river runoff. Figure 7 shows examples of the spectral densities S(f) of daily runoff for the Upper Yangtze, Lena and Darling rivers characterized by different climatic conditions. The spectral densities for the other rivers are given in Figure A2 (available with the online version of this paper).

As can be seen in Figure 7, the spectral density of river runoff of the Lena and Yangtze has, as was expected, a clearly pronounced maximum at a frequency f_an = 2.738 × 10⁻³ day⁻¹ (marked with the number 1) corresponding to a temporal period of one year. Such a maximum is absent for the Darling River. This will be explained below.

When analysing the spectral densities of runoff for the 11 rivers (including those listed in Figure A2), it can be seen that the spectral density does not have an annual peak only for the Darling and Tagus rivers. This means that in these cases a regular component has a small effect on runoff formation. For the remaining rivers, the highest peaks are mainly at the frequencies corresponding to the annual and semiannual periods. The rivers with a clearly pronounced spring flood or the monsoon character of precipitation (the Lena, Mackenzie, Upper Amazon, Upper Yangtze, Upper Huanghe, Niger and Ganges) have, as a rule, the largest amplitude of the main peaks. The Upper Mississippi and Rhine rivers also have peaks 1 and 2, but their amplitude is smaller than for the above-mentioned group of rivers and these two rivers have an intermediate position between the group of the above seven rivers and the Darling and Tagus rivers. All that has been said is vividly confirmed by Figure 7.

Figure 8 shows the relationship between RSR and uncertainty of annual runoff Un_an (Figure 8(a)) or mean uncertainty of monthly runoff Un_mn,aver (Figure 8(b)). The obtained dependencies show that when the regular component in runoff formation increases the river runoff uncertainty becomes lower both on annual and monthly scales. Thus, Figure 8 is in good agreement with the above qualitative analysis of the dependence of runoff uncertainties for different rivers on the annual periodicity of runoff formation.

The most outstanding river is the Darling. As mentioned above, the spectral density of runoff S(f) for the Darling River does not have a maximum at the frequency f_аn (nevertheless, for the numerator in the RSR criterion, the value of S(f_an) was used). Such a situation is extremely rare and occurs only in the case of very strong stochasticity. In Figure 8, it is clearly seen that the Darling has the lowest value of RSR, i.e., its flow is characterized by an extremely high degree of stochasticity. In fact, the Darling River flow is known to be very irregular. Its numerous tributaries have sources in the arid inner part of Australia, and many of them reach the main watercourse only in occasional periods of flood. In the dry season, the river in its lower course dries up and breaks up into separate river parts. During a drought, which sometimes lasts more than a year, the Darling becomes much shallower, transforming into a stunted brook. However, although rare, in summer the river can become damaging.

The simulations of streamflow of 11 rivers with a priori input data including meteorological variables and land surface parameters derived from the global WATCH and ECOCLIMAP data sets were performed using the SWAP model. The simulated river runoff showed poor agreement with measurements, motivating calibration of model parameters influencing runoff formation to the most extent. Calibration of selected parameters (including correcting factors to incoming radiation) for all river basins was performed against measured values of monthly runoff for several years. Application of the calibrated parameters substantially improved the results. Thus, for the validation period, the coefficient of efficiency for nine rivers ranged from 0.58 to 0.9 (averaged to 0.78), while mean absolute bias was within the range from 0.5% to 13.8% (averaged to 4.4%).

Analysis of the obtained results shows that simulation of river runoff by the SWAP model in combination with input data derived from the global WATCH and ECOCLIMAP data sets can be considered (after calibration of the most influencing parameters) as fairly good for river basins located under different natural conditions. This conclusion is of great importance because meteorological observations and land surface parameters are not available for many regions of the globe and in such cases global data sets may be the only source of information for hydrological calculations.

Natural variability of river runoff caused by the impact of weather noise (associated with a chaotic character of atmospheric processes resulting in stochastic instability of meteorological characteristics of near-surface layer of the atmosphere and determining the lowest limit of uncertainty in the estimates of meteorological variables simulated by climatic models) was estimated on two time scales (annual and monthly) for the selected rivers. This uncertainty can be treated as the lowest limit of predictability of river runoff.

Analysis of the estimated uncertainties has revealed that monthly uncertainties of river runoff are larger than annual ones. In addition, it was shown that differences in runoff uncertainties obtained for different rivers depend greatly on natural conditions of a river basin, in particular, on the ratio of deterministic and random components of the river runoff formation.

The obtained values of natural uncertainty of both annual river runoff and climatic hydrographs, as well as the dependencies of the uncertainties on regularity-to-stochasticity ratio can be used for comparison with possible values of projected changes in annual and monthly runoff due to climate change and increasing impact of human activity. This study creates a sound basis for application of climate scenarios and projecting the dynamics of river runoff for the selected river basins by the end of the 21st century.

The work was supported by the Russian Science Foundation, project no. 16-17-10039. River runoff data were kindly provided by ISI-MIP organizers and the Global Runoff Data Centre (GRDC), D-56068 Koblenz, Germany.