Can multi-objective calibration of streamflow guarantee better hydrological model accuracy?

In order to understand hydrological processes and thus predict the likely influence of global climate change on hydrology, a means of extrapolating into the past and the future has been deemed necessary, mainly due to the limitations in space and time of hydrological measurement techniques and scope (Beven 2012). There are many types of models that can enable us to quantitatively extrapolate, predicting flooding, runoff, the availability of various water resources and nutrient transport (Mimikou et al. 2015; Molina-Navarro et al. 2017). However, these models are only simplified representations of the real world, and their accuracy cannot be taken for granted (Muleta 2012). In addition, some model parameters may not physically manifest and therefore cannot be determined through direct measurements. Thus, model calibration is always necessary (Kim & Lee 2014). Model calibration is the process of identifying model parameter values from the available parameter ranges to enable the maximization or minimization of the objective function (a function representing difference between simulated and observed values) (Muleta 2012).

There has been considerable research into various optimization algorithms used to identify optimal model parameter values. In the past, this particular research has mainly focused on the calibration of a single objective (Vrugt et al. 2003a). For example, a widely used algorithm for identifying optimal parameter values with a single-objective function, the shuffled complex evolution, was developed at the University of Arizona (SCE-UA) (Kan et al. 2016). This type of calibration algorithm is primarily concentrated with matching one component of the hydrograph produced by the model to the observed data (Boyle et al. 2000). However, in regards to the observed data, it has been argued that any single-objective function, no matter how elaborately chosen, is usually inadequate to properly represent all of the characteristics (Vrugt et al. 2003a; Sadeghi-Tabas et al. 2017). To circumvent this problem, multiple-objective calibration, which represents more than one aspect of a hydrological system and its behaviour, has been proposed as an alternative as it can incorporate more information from the available data (Zhang et al. 2015).

In the past, assumptions have been made that multi-objective calibration improves model accuracy (Kim & Lee 2014). Nevertheless, little research has gone into testing this supposition prior to now, which this study undertakes. To carry out this test, the Nash–Sutcliffe efficiency and its logarithmic form were selected as the objective functions for both multi- and single-objective calibration for the FLEX-Topo (Topography-driven FLux EXchange model) used in the upper stream of the Heihe River basin (UHRB) in northwest China. At the basin outlet, streamflow discharge was used for model calibration. Model accuracy was assessed by comparing the simulated streamflow with the observed streamflow data.

The Upper Heihe River Basin (UHRB) is the upper stream of the second largest inland river in Northwest China, the Heihe River (Figure 1(a)). The river originates in the Qilian Mountains. The UHRB has an area of 10,009 km² and stretches for 303 km. The elevation of the UHRB ranges from 1,700 m to 5,000 m (Figure 1(b)). This mountainous region surrounding the UHRB is characterized by a cold desert climate, with an annual temperature of UHRB 2–3 °C (Liu et al. 2008). The long-term average annual precipitation is ∼430 mm y⁻¹ and potential evaporation is ∼520 mm y⁻¹. More than 80% of the annual precipitation occurs from May through September. Soil types in this region are predominantly mountain straw and grassland soil, chestnut-coloured soil, chernozemic soil and desert. The land cover includes forest, grassland, bare rock or bare soil, wetland and permanent snow (Figure 1(c)). On average, the UHRB produces ∼70% of the total river-borne runoff of the whole Heihe River Basin (Chen et al. 2005). As it is the main runoff-producing region in the Heihe River basin, the UHRB is essential for an integrated water management practice, and so the hydrodynamics of the UHRB have garnered a great deal of research interest.

One hydrological station and four meteorological stations are located in and around the UHRB (Figure 1(c)). For this study, these meteorological stations provided input data for the hydrological model. In this model, Thiessen polygons were used to define the monitored area of these meteorological stations. Since the meteorological stations are all located at relatively low elevations, empirical formulas were used to adjust precipitation and temperature data. The streamflow discharge data for model calibration were obtained from the hydrological station (Yingluoxia Station) at the outlet of the UHRB to calibrate the model. All these input data were obtained from the Cold and Arid Regions Science Data Center (http://westdc.westgis.ac.cn/).

The FLEX-Topo model, a semi-distributed hydrological model, was used as the basis for different calibration strategies. Input data included daily precipitation and daily temperature, and the outputs included streamflow discharge and evapotranspiration. The FLEX-Topo model consists of four identically structured, parallel components representing four hydrological landscapes (hydrological upland, sunny-hillslope, shady-hillslope and hydrological lowland) which are classifed through topographic information that included height above the nearest drainage (HAND), elevation, slope and aspect (Gao et al. 2014). There are 25 parameters, shown in red in Figure 2. Table 1 shows the definition of the parameters and the ranges of their values determined by experience or references.

The four landscapes are each characterized by various runoff-producing mechanisms, which are embodied in the parallel components of the FLEX-Topo model. The hydrological processes in each landscape are briefly described next. More details on the model parameters, the model structure and the water balance equations have been described by Gao et al. (2014).

In the hydrological upland, the dominant hydrological processes are Hortonian overland flows (HOF) and deep percolation (DP). Precipitation (P) can be stored as snow cover (S_wB) if the daily temperature is below the threshold value T_t. It is assumed that there is no interception due to a lack of significant vegetation. The sum (P_eB) of rainfall and the snowmelt, which is calculated by a degree-day model, flows towards the unsaturated zone reservoir (S_uB). When P_eB is greater than the infiltration capacity (P_t), Hortonian overland flow (R_HB) occurs. Infiltration (R_uB) recharges the unsaturated reservoir (S_uB). Percolatation (R_pB) from S_uB that flows into the slow-response reservoir (S_s) is calculated from the relative soil moisture (S_uB/S_umaxB) and maximum percolation (P_maxB). Actual evapotranspiration (E_aB) from the unsaturated reservoir is estimated by the relative soil moisture (S_uB/S_umaxB) and the potential evapotranspiration that is calculated by the Hamon equation (Hamon 1960). Saturation excess overland flow (R_eB) occurs if the amount of water in the unsaturated reservoir exceeds the maximum storage capacity (S_umaxB). R_eB plus R_HB is partitioned by a splitter parameter (D_B) into flow (R_fB) going into the fast-response reservoir (S_fB) and flow (R_sB) re-infiltrating into the slow-response reservoir (Ss). R_fB turns into R_LfB after convolution using the time lag parameter T_lag, which represents the time interval between storm and fast runoff generation. Flow (Q_fB) from S_fB is routed to the stream channel after time K_fB.

The major runoff-producing mechanism of shady-hillslope is the storage excess subsurface flow (SSF). The existence of vegetation indicates the necessity for the interception reservoir (S_iFH). Evapotranspiration (E_iFH) from the interception reservoir is assumed to be equal to potential evapotranspiration if the storage of the interception reservoir is nonzero. The sum (P_eFH) of the remainder of rainfall after interception and the snowmelt that is calculated from a degree-day model flow towards the unsaturated reservoir (S_uFH). P_eFH is partitioned into runoff and flow that is routed to S_uFH through the use of the runoff coefficient that is calculated from the parameter β_FH representing the heterogeneity of soil properties. Runoff from P_eFH and S_uFH is separated by a splitter parameter (D) into flow (R_fFH) going into the fast-response reservoir (S_fFH) and flow (R_sFH) re-infiltrating into S_s. Actual evapotranspiration (E_aFH) from the S_uFH is estimated by parameter C_e and potential evapotranspiration, the latter of which is calculated by the Hamon equation. C_e is a threshold value. If S_uFH/S_umaxFH is larger than C_e, actual evapotranspiration is assumed to be equal to potential evapotranspiration. R_fFH becomes R_LfB after convolution using the time lag parameter T_lag and flows into the fast-response reservoir (S_fFH). Water (Q_fFH) from S_fFH is routed to the stream channel with timescale (K_fH).

The major runoff-producing mechanism of sunny-hillslope is assumed to be the same as that of shady-hillslope, but the parameters for this landscape take different values.

In hydrological lowland, the dominant hydrological process is saturation excess overland flow (SOF). Interception and snowmelt are the same as they are in shady-hillslope, but soil routine is different. Runoff (R_fW) produced from P_eW and S_uW is directly routed to the fast-response reservoir (S_fW) without delay because of its proximity to the stream channel. Due to the shallow groundwater level and limited storage capacity, capillary rise (C_R) occurs. C_R is represented by a parameter (C_Rmax) that indicates a constant amount of capillary rise.

The MOSCEM-UA (Multi-Objective Shuffled Complex Evolution Metropolis-University of Arizona) algorithm that was developed by Vrugt et al. (2003a) was proposed for multi-objective calibration because it maintains a uniform sampling density within the Pareto set of solutions and includes the single-criterion end points. The evolution strategy employed in the MOSCEM-UA algorithm is similar to that of the SCEM-UA algorithm (Vrugt et al. 2003b), but the probability ratio concept in the SCEM-UA algorithm was replaced by a multi-objective fitness assignment concept in order to develop the initial population of points toward a set of solutions resulting from a stable distribution. The MOSCEM-UA algorithm has been successfully applied in a wide range of hydrological and environmental models, demonstrating its reliability and effectiveness (Efstratiadis & Koutsoyiannis 2010). In this study, the MOSCEM-UA algorithm (Vrugt et al. 2003a) was adopted to identify the parameters in the model (Figure 3). With these parameters, simulated streamflow is generated. Using the simulated and observed streamflow, objective function values were calculated by a module within the FLEX-Topo model.

Three basic parameters had to be set in order to control the operation of the algorithm: the number of iterations, the number of parallel sequences, and the number of random samples. Kuczera (1997) suggests that the number of parallel sequences be set to be equal to the number of the calibration parameters. The other two algorithmic parameters were set according to the recommended values in Vrugt et al. (2003b). In line with these suggestions, these three parameters were set as 100,000, 25 and 2,500, respectively. The 20-year-period from 1959 to 1978 was used for model calibration. The first year was disregarded in order to minimize any influence from the initial conditions. A computer having a 3 GHz Intel Xeon CPU model E5-2687Wv4 was used to run the algorithms.

There are many commonly used objective functions, such as the Nash–Sutcliffe efficiency (Nash & Sutcliffe 1970) and its logarithmic form (Legates & Mccabe 1999; Kim & Lee 2014). The Nash–Sutcliffe efficiency and its logarithmic form highlight high flow and low flow, respectively. Therefore, both indexes were included for multi-objective calibration in order to simultaneously consider both high flow and low flow. Meanwhile, these two indexes were also separately used for single-objective calibration. Since the MOSCEM-UA algorithm generally deals with a minimization procedure, the Nash–Sutcliffe efficiency and its logarithmic form were subtracted from one. Table 2 shows the mathematical formulation of the objective functions. Comparison between different calibration strategies was made to assess the simulation accuracy in modelling observed streamflow, focusing on high flow and low flow.

For single-objective calibration, one set of parameters was obtained, and for multi-objective calibration the Pareto front consisting of 33 sets of parameters was obtained. Figure 4 shows the results of different calibration strategies.

It is clear that several tradeoffs exist in multi-objective calibration, indicating that no single solution can optimize both objective functions simultaneously. The values of NSE and lnNSE corresponding to solutions in the Pareto front (dots marked in blue) vary widely. The NSE values range from 0.154 to 0.199, while the lnNSE values range from 0.100 to 0.151, indicating good performance of the model.

In regards to the NSE calibrated model (dot marked in red), the NSE value is 0.157, which is very close to the value (0.154) of the corresponding extreme end of the Pareto front. As regards the lnNSE calibrated model (dot marked in green), the lnNSE value is 0.097 and it is close in value to the other extreme end of the Pareto front (0.100). This proximity in value confirms that the Pareto front includes a single-criterion end point, which is consistent with the results of Vrugt et al. (2003a). In addition, it can be found that the extreme end of the Pareto front on lnNSE is closer to the single-objective calibration, yet the other extreme end is not proximal to the NSE calibrated solution.

Since all Pareto optimal solutions are equally important, it may be difficult to prefer one solution over another. Therefore, every solution in the Pareto front was used to run the model, obtaining 33 sets of simulated streamflow. The average of these simulated flows was calculated and was denoted as the ‘Average’. The hydrographs of the different objective calibrated simulations are also shown in Figure 5. The simulated flows are respectively marked as Q_NSE, Q_lnNSE and Q_Average.

It can be seen that the simulated flows fluctuated synchronously with the observed flow, which indicates that both strategies could capture the dynamics of the observed discharge. It seems, however, that the resulting hydrograph based on NSE appeared slightly truer to the observed high flow than the other two resulting hydrographs.

Figure 6 shows a scatterplot of the residuals in regards to observed flow. For low flows, these three simulations all tend to show a more centred scatterplot and the distributions of the dots were similar. This pattern indicated that these three simulations performed similarly in modelling low flow, although it was supposed that the lnNSE calibrated simulation would behave better. For higher flows the dots show a rising trend, demonstrating that all three simulations underestimated the high flows but to various degrees. The multi-objective calibrated simulation performed no better than the simulation based on NSE, the latter which presented a relatively low degree of underestimation. Taking the values of NSE into consideration, it is interesting to find that the lower the NSE value, the better the modelling of the high flow. In general, the NSE tends to perform better in modelling high flow and in reproducing the observed streamflow dynamics. There were no apparent differences among these three simulations in simulating low flow.

The MOSCEM-UA algorithm is a random sampling approach, thus initial samples may have an impact on the solutions. In addition, three parameters control the behaviour of MOSCEM-UA: the number of iterations (NI), the number of random samples (NRS) and the number of parallel sequences (NPS). To examine the influence of samples and these three parameters on the simulations, various combinations of these algorithm parameters were carried out. The results of these trials are shown in Figure 7. It seems that trials of the changing parameter combinations of the optimization algorithm almost made no difference in modelling high flow and low flow. The results of these trials also prove the simulations based on NSE outperformed the simulations based on multi-objective.

This study aimed at testing whether multi-objective calibration would guarantee a better simulation of observed streamflow. Two widely used objective functions, NSE and lnNSE, were selected for model calibration. Since all the settings, including the input data, the hydrological model, the parameter ranges and the optimization algorithm were kept the same (except for the calibration strategies), it can reasonably be assumed that it was the objective functions that influenced model performance.

For high flow, the results show that simulation based on NSE performed slightly better. This may be because NSE weighs error on high flow more than error on low flow. Because of the use of the logarithmic form, lnNSE weighs more on the error on low flow (Fenicia et al. 2007; Kim & Lee 2014). Therefore, the simulation calibrated by NSE outperformed the simulation based on lnNSE in reproducing high flow. Since the multi-objective calibration takes both NSE and lnNSE into consideration, it is not surprising that simulation based on NSE outperformed simulation based on multi-objective calibration in modelling high flow. By contrast, there was no significant difference among the three types of simulations based on different objective function(s), although simulation based on lnNSE was supposed to improve the low flow modelling. These results may be due to two factors. The first is that NSE still reacts to low flows although it highlights high flow. The second factor is that there generally tends to be more measurement errors in the observed low flow values. (Laaha et al. 2006). The greater level of error in observed low flow can be confirmed by Figure 6, in which the centred scatter illustrates a nearly symmetrical distribution around the line whose ordinate value is equal to zero.

It can be argued that a great deal of attention should be paid to objective function selection for multi-objective calibration as objective function has impact on model performance (Muleta 2012). It is crucial that performance of objective function be assessed before application in multi-objective calibration.

It has been suggested that the inclusion of more variables in model calibration might provide additional information about the catchment characteristics (Zhang et al. 2016). In this study, the variable was streamflow. However, it is worth mentioning that there are other observed variables that could be used for multi-objective calibration, such as time shift between observed and simulated peak flow (Fenicia et al. 2007), energy (Scheerlinck et al. 2009) and chemical constituent (Haas et al. 2016; Zhang et al. 2016). For example, Zhang et al. (2016) calibrated both water quantity and water quality simultaneously. This type of multi-objective calibration was not included in this particular study, but it may prove to be an area of interest in the future.

This hydrological modelling study aimed at testing the supposition that a multi-objective calibration guaranteed better model accuracy in modelling high flow and low flow. To capture the dynamics of high and low flow, the Nash–Sutcliffe efficiency and its logarithmic form were selected as the objective functions. These two criteria were both included in the multi-objective calibration, while in the single-objective calibration they were separately used. The FLEX-Topo, a semi-distributed hydrological model, was used to simulate the Upper Heihe River basin hydrological processes, and the MOSCEM-UA was used for calibration. Assessment of different calibration strategies was carried out by comparing the simulated streamflow with the observed streamflow, focusing on the high flow and low flow.

It seems that multi-objective calibration, in fact, did not significantly improve model performance compared to single-objective calibration. Our results show that the Nash–Sutcliffe efficiency performed slightly better in modelling high flow while the logarithmic form of the Nash–Sutcliffe efficiency did not significantly improve the low flow simulation and at the same time more seriously underestimated the high flow. Multi-objective calibration led to no better accuracy in modelling high flow and similar simulation in low flow compared to single-objective calibration which was based on the Nash–Sutcliffe efficiency. Thus, it cannot be taken for granted that multi-objective calibration of streamflow will necessarily guarantee better model accuracy. Indeed, model performance depends on the selected objective function(s). With this in mind, in order to design better and more accurate hydrological simulations of hydrological processes, a serious look at the selection of the objective functions should be undertaken.

This work was supported by Natural Science Foundation of China (91325302; 91425303), the National Natural Science Foundation of China (41571022, 41625001), the National Science Fund for Distinguished Youth Scholars (41625001), the Beijing Natural Science Foundation Grant (8151002), and the National Science and Technology Major Project (2015ZX07203-005). We thank Li Wang for her support during this study. We also thank Kun Ma for his constructive advice and Yvonne Cranmer for her suggestions on language issues.