Evaluating the Effect of Numerical Schemes on Hydrological Simulations: HYMOD as A Case Study

The importance of numerical schemes in hydrological models has been increasingly recognized in the hydrological community. However, the relationship between model performance and the properties of numerical schemes remains unclear. In this study, we employed two types of numerical schemes (i.e., explicit Runge-Kutta schemes with different orders of accuracy and partially implicit Euler schemes with different implicit factors) in the hydrological model (HYMOD) to simulate the flow hydrograph of the Leaf River basin from 1948 to 1988. Results computed by different numerical schemes were compared and the relationships between model performance and two scheme properties (i.e., the order of accuracy and the implicit factor) were discussed. Results showed that the more explicit schemes generally lead to the overestimation of flow hydrographs, whereas the more implicit schemes lead to underestimation. In addition, the numerical error tended to decrease with increasing orders of accuracy. As a result, the optimal parameter sets found by low-order schemes significantly deviated from those found by the analytical solution. The findings of this study can provide useful implications for designing suitable numerical schemes for hydrological models.


Introduction
Hydrological models are widely used in catchment studies and water resource management because they provide a reasonable representation of hydrological processes [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. Over the past few decades, attempts have been made to establish a theoretical framework for developing hydrological models [17][18][19][20]. In general, the model building process can be divided into three steps including the development of: (1) a conceptual model, (2) a mathematical model and (3) a computational model. The first two steps formulate the mathematical governing equations for the hydrological system, while the last step generates numerical solutions to the governing equations with numerical schemes. Although numerical schemes have been extensively explored in many fields, such as aerodynamic and hydrodynamic modeling [21,22], their importance has not been well recognized in hydrological sciences.
A key problem in developing hydrological models is whether the structure of the model can adequately reflect the types of hydrological characteristics [2]. Errors in the modeling results may arise from model uncertainty, which can be identified as the sum of three types of uncertainties including: (1) the theoretical model uncertainty, which refers to the uncertainty due to simplification and approximation in conceptualization and mathematical representation, (2) the parameter uncertainty, which refers to the uncertainty in model parameterization and (3) the numerical scheme and solution uncertainty, which refers to the uncertainty associated with the numerical schemes and solution algorithms [23,24]. The error induced by numerical scheme and solution uncertainty is termed numerical error, as opposed to structural error or sometimes termed structural inadequacy [19], that is associated with the theoretical model uncertainty and parameter uncertainty. Inadequate numerical schemes may induce excessive numerical errors, which complicates the structural identification process and prevents us from getting the right answer for the right reasons [18]. In recent decades, the importance of numerical schemes has been increasingly recognized. Newly emerged hydrological models, such as the SUPERFLEX modeling framework [13,14], the Framework for Understanding Structural Errors (FUSE) [25], the Structure for Unifying Multiple Modeling Alternatives (SUMMA) [26] and the Weather Research and Forecasting (WRF) modeling systems [27], have identified the implementation of numerical schemes as a fundamental step in their modeling frameworks. The performances of different numerical schemes have also been explored. For example, Clark and Kavetski [25] and Kavetski and Clark [28] analyzed the performances (e.g., fidelity, computational cost and time step dependence) of the fixed-and adaptive-step explicit/implicit schemes, respectively, in six conceptual rainfall-runoff models and assessed their impacts on sensitivity analyses, model optimizations and parameter inferences. They found that the fixed-step explicit schemes tended to produce severe numerical errors that misleadingly compensated for the structural error and proposed the use of implicit schemes for smoother solutions. Other numerical schemes have also been proposed based on the comparison between a number of selected numerical schemes [29]. However, a generalized relationship between the model performance and the properties of numerical schemes is still not available.
Numerical schemes have two important properties: the order of accuracy and the implicit factor. The order of accuracy is closely associated with the magnitude of the numerical error. With a fixed time step and grid size, high-order numerical schemes generally produce smaller numerical errors and are, therefore, less likely to degrade the model performance. However, Clark and Kavetski [25] and Kavetski and Clark [28] suggested that high-order explicit schemes could not completely prevent the problems of spurious multimodality and inconsistent parameter inference. On the other hand, the implicit factor, which is not directly related to the magnitude of the numerical error, also strongly affects the model performance. Clark and Kavetski [25] compared several explicit and implicit schemes that were of the same orders of accuracy but with different implicit factors and found that their performances were very different. Therefore, the relationship between model performance and the properties of the numerical schemes should be carefully studied in order to develop reliable hydrological models.
Given these considerations, this study aims to evaluate the effects of numerical schemes on hydrological models using the hydrologic model (HYMOD) as an example. The original HYMOD employs a low-order explicit scheme to calculate the streamflow hydrograph, which may also suffer from serious numerical error that leads to poor prediction and misleading diagnostic information [25]. By evaluating the performances of different numerical schemes to be implemented in HYMOD through extensive numerical tests, we reveal how the model predictions can be improved by controlling numerical error via the manipulation of solution algorithms. Instead of merely comparing the performances of a number of selected numerical schemes, we focus on finding the generalized relationship between the model performance and the two common properties of numerical schemes (i.e., the order of accuracy and implicit factor). The results of this study are expected to demonstrate the importance of numerical schemes in reducing the uncertainty in hydrological models and quantify the effects of scheme properties on model output, which helps to improve our understanding of numerical schemes and provides a reference for developing suitable numerical schemes for hydrological models. This study is organized as follows. Section 2 describes the data, the mathematical formulation of HYMOD and the numerical schemes used in this study. Sections 3 and 4 present the results of the HYMOD simulations, sensitivity analysis and parameter optimization using two classes of numerical schemes with different orders of accuracy and implicit factors and discuss their effect on model prediction. Section 5 summarizes the major findings of this study.

Data
Hydrological data from the Leaf River basin were used in this study. The Leaf River basin is a humid catchment located north of Collins, Mississippi, with a drainage area of 1950 km 2 [17]. The dataset contains forty consecutive years of mean daily precipitation, potential evapotranspiration and streamflow from 1 October 1948 to 30 September 1988, which can represent a wide variety of hydrological conditions.

Model Configuration
HYMOD is a six-parameter hydrological model running on MATLAB platform (R2015a, MathWorks, Inc., Natick, MA, USA) that simulates the rainfall-runoff process of watersheds [30]. The six parameters include the maximum storage of the soil moisture tank C par , the shape parameter b, the ratio of quick flow to total flow α, the number of linear tanks with quick flow N q , the quick flow coefficient K q and the slow flow coefficient K s . HYMOD has a relatively simple model structure, thereby making it suitable for testing the impact of numerical schemes on model output. The configuration of HYMOD is shown in Figure 1. The model contains a soil moisture accounting (SMA) module and a surface flow routing (SUFR) module. A brief description of the mathematical formulations in HYMOD is presented in the following section.

Mathematical Description
SMA converts rainfall excess into the inflow of SUFR. The mass conservation equation for SMA can be expressed as: where Q SMT represents outflow from the soil moisture tank, P represents precipitation, E represents evapotranspiration and C represent the filled storage of the soil moisture tank, which is assumed to be a function of the water depth of the soil moisture tank, H and takes the following form: where H par represents the maximum water depth of the soil moisture tank and b is a shape parameter.
H is used as an indicator of the saturation degree of the river basin, which resembles the groundwater level, soil moisture and so forth. A warm-up period can be used in the simulation when no data can be obtained to prescribe an initial condition for H. SUFR has two components: quick flow routing (QFR) and slow flowing routing (SLFR). QFR can be conceptualized as a Nash cascade process that consists of multiple linear reservoirs with the same outflow coefficient (K q ). The storage equation of the linear reservoir with storage X q , inflow I q and outflow O q can be expressed by: For a Nash cascade process with n linear reservoirs, the mass conservation can be expressed in a vector form as [31]: SLFR is similar to QFR except that only one linear reservoir is utilized; therefore, the governing equation of SLFR is similar to Equation (3) and can be expressed as: where X s , I s , O s and K s are the storage, inflow, outflow and outflow coefficient for the linear reservoir in SLFR. An additional parameter, α, is introduced to determine the fractions of Q SMT that flow to QFR and SLFR using the following equation: Equations (1), (2) and (4)-(6) form the governing equation set for HYMOD.

Analytical Solution to SMA
No numerical solutions are necessary in solving the governing equation for SMA (Equation (1)) because the outflow depends only on the initial and final state of the basin storage (P and E are user-prescribed parameters). In deriving an analytical solution to SMA, an additional equation is required to close the governing equation. Moore [30] assumed the depth of soil moisture tank was increased by rainfall, was depleted by evapotranspiration and generation of direct runoff occurred when rainfall exceeded the storage capacity. Therefore, the following relationship is assumed: Equation (7) indicates that all rainfall excess infiltrates into the soil moisture tank and water depth in the tank rises correspondingly. By taking the derivatives with respect to H on both sides of Equation (2), we obtain: By substituting Equations (7) and (8) into Equation (1), we obtain the analytical form of Q SMT : Equation (9) is the analytical form for Q SMT in this study. It should be noted that Equation (9) is valid only when H < H par . Once H exceeds H par , infiltration no longer occurs and all rainfall excess is converted into surface runoff.

Numerical Schemes
In this section, we formulate the numerical schemes for discretizing the storage equations of the linear reservoirs in both QFR and SLFR. The storage equation of a linear reservoir can be expressed as: By integrating on both sides of the storage equation from time t to t + ∆t, we obtain: For the ith reservoir in a Nash cascade process, Equation (10) can be written as: The integral of X(t) from t to t + ∆t needs to be evaluated with numerical schemes. The explicit Runge-Kutta schemes and the partially implicit Euler schemes are briefly introduced in the following sections.

The Explicit Runge-Kutta Schemes
The explicit Runge-Kutta schemes are a class of numerical schemes with zero implicit factor and variable orders of accuracy. In this study, the first-to fourth-order explicit Runge-Kutta schemes are introduced to evaluate the integral of X(t), which can be expressed as follows.
The first-order explicit Runge-Kutta scheme (RK1): The second-order explicit Runge-Kutta scheme (RK2): The third-order explicit Runge-Kutta scheme (RK3): The fourth-order explicit Runge-Kutta scheme (RK4): where O(∆t p ) represents the truncation error terms, which is order p, resulted from using the numerical schemes to approximate the integrals. According to Equation (3), the derivative dX i dt can be evaluated by:

The Partially-Implicit Euler Schemes
A generalized form of the partially-implicit Euler schemes can be expressed as: where θ represents the implicit factor ranging from 0 to 1. Equation (18) indicates that the partiallyimplicit Euler schemes evaluate the temporal integration time t to t + ∆t as a weighted average of the integrals at the previous (t) and the current (t + ∆t) time steps. When θ = 0, the scheme is called the fully explicit Euler scheme (EXP), which has exactly the same form as RK1 (Equation (13)); when θ = 1, the scheme is called the fully implicit Euler scheme (IMP); when 0 < θ < 1, the scheme is called the partially-implicit Euler scheme (PIM). EXP and IMP can also be viewed as two special types of partially-implicit Euler schemes; therefore, we use the term partially-implicit Euler schemes to refer to all three types of numerical schemes in the rest of this study. The partially-implicit Euler schemes generally have first-order accuracy, except for θ = 1/2. The original HYMOD uses EXP by default.

Analytical Solutions to QFR and SLFR
Both QFR and SLFR can be viewed as the Nash cascade process. Szöllösi-Nagy [31] derived an analytical solution for the Nash cascade process, in which the outflow rate of the ith linear reservoir at time t + ∆t, O i (t + ∆t), can be expressed in terms of its storage at time t as: with where k represents the outflow coefficient and I i represents the inflow of the ith linear reservoir. In the present study, the analytical solution to HYMOD is obtained from Equations (9) and (19).

The Downhill Simplex Method (DSM) for Parameter Optimization
The DSM is a nonlinear optimization technique, which was employed in this study to find optimal parameter combinations that correspond to functional minimums. The optimization starts with an initial simplex consisting of m vertices, each of which represents a randomly generated parameter set. The functional minimum in the parameter space can be found by transforming the simplex according to the functional values at the vertexes and moving the simplex in the downhill direction until the functional value converges to its minimum value. Readers may refer to Nelder and Mead [32] for more details about the DSM.

Model Setup
In this study, we tested the influence of numerical schemes on the output of HYMOD with a series of numerical tests including the simulations of the Nash cascade process, the flow hydrograph of Leaf River basin with a priori parameter set and 2000 random parameter sets, the one-dimensional and two dimensional sensitivity analyses on model parameters and the DSM optimization test.
The Nash cascade process is a basic component of QFR and SLFR. Therefore, we simulated the Nash cascade process to test the influence of numerical schemes on QFR and SLFR. The simulation was carried out with one to ten linear reservoirs in the Nash cascade process. The storage coefficients K q of the reservoirs were set identical at 0.5. One unit volume of pulsed inflow was given to the first reservoir and the outflow hydrograph of the last reservoir was simulated and compared with the analytical solution (Equations (9) and (19)) to assess the performance of the numerical schemes.
The simulation on Leaf River basin was carried out first with a priori parameter set (Table 1) and then with 2000 randomly generated parameter sets that were generated from a uniform distribution in the parameter space. In the priori parameter set, the flow hydrographs simulated by different numerical schemes were compared with that computed by the analytical solution and the observed data to assess the influence of numerical schemes on model output. In the random parameter test, the mean squared error (MSE) was used as a measure of the overall scheme performance over a wide range of model parameter. The entire 40-year sampling period (1 October 1948 to 30 September 1988) was simulated but only the results from January 1970 to January 1972 were plotted for better visibility. The response sensitivity of the six model parameters (C par , b, α, N q, K s and K q ) was assessed by the one-dimensional sensitivity analysis (1DSA) and two-dimensional sensitivity analysis (2DSA). The sensitivity analysis was carried out by varying the values of the parameters being tested within their feasible ranges (Table 1) while holding the other parameters fixed at their priori values.
In the DSM optimization, five model parameters (C par , b, α, K s and K q ) were optimized to find the best combination of model parameters that minimized the MSEs of the computed flow hydrograph. N q was not utilized in the optimization because it was a discrete variable whose value could not be arbitrarily prescribed. Eight initial parameter sets (Table 2) were used to test the performances of the numerical schemes in different regions of the parameter space. A total of 300 steps were carried out during each DSM optimization.

The HYMOD Simulations for the Leaf River Basin with a Priori Parameter Set
The flow hydrographs of the Leaf River basin simulated with the priori parameter set (Table 1) by the explicit Runge-Kutta schemes (RK1 to RK4) were compared with the analytical solutions (Equations (9) and (19)) and the observed data in Figure 2. The analytical solutions deviated from the observed data likely due to inadequate representation by the priori parameter set. Compared with the analytical solutions, RK1 overestimated the flow in the Leaf River basin during high flow conditions and underestimated the flow during low flow conditions, whereas RK2, RK3 and RK4 all yielded better fits to the analytical solution. This is consistent with the findings in the simulation of the Nash cascade process that the first-order explicit scheme is prone to overestimation whereas second-and-higher-order schemes may effectively reduce the error. Similarly, the solutions by the partially implicit Euler schemes (EXP, PIMs with θ = 1/3 and 2/3 and IMP) were compared with the analytical solutions and observed data in Figure 3. The differences between the numerical and analytical solutions varied with the implicit factor θ. The more explicit schemes with smaller implicit factors tended to generate solutions with greater fluctuations in the hydrograph, characterized by higher peaks and deeper troughs, than the more implicit Euler schemes with larger implicit factors and vice versa. In both Figures 2 and 3, some numerical schemes produced hydrographs that were even more close to the observed data than that simulated by the analytical solution, which indicates an interaction between numerical schemes and model parameters. In summary, results from the simulation of Leaf River basin confirm the strong influence of numerical schemes on the output from HYMOD.

The HYMOD Simulations with Random Parameter Sets
The MSEs of the flow hydrographs simulated by the explicit Runge-Kutta schemes (RK1 to RK4) and the partially implicit Euler schemes (EXP, PIMs with θ = 1/3 and 2/3 and IMP) using the 2000 random parameter sets were compared with those of the analytical solutions in Figure 4. The 1:1 slope in Figure 4 represents the numerical solutions that have equivalent performances as the analytical solutions in terms of MSEs. The numerical solutions computed by RK1 showed a considerable scatter around the 1:1 slope, indicating that RK1 is a poor approximation of the analytical solution over a wide MSE range. In contrast, the numerical solutions computed by RK2 were close to the analytical solution for smaller MSE values but became divergent with increasing MSE values, while RK3 and RK4 provided good approximation to the analytical solution over the whole parameter range, which indicates that high-order explicit Runge-Kutta schemes were more reliable in approximating the analytical solutions than their low-order counterparts. For the partially implicit Euler schemes, the implicit factor seemed to determine whether the MSEs produced by the numerical schemes were larger or smaller than those by the analytical solutions. Results showed that the MSEs showed a tendency to decrease with increasing implicit factor over the entire parameter range: EXP (RK1) and IMP produced larger and smaller MSEs than the analytical solutions, respectively and PIMs with θ = 1/3 and 2/3 produced MSEs that were similar to the analytical solutions. Therefore, the two properties of a numerical scheme, that is, order of accuracy and implicit factor, may influence the simulation results in different ways: the latter may determine the direction in which the numerical solution deviates from the analytical solution, whereas the former may determine how far it deviates. The lowest MSE computed by a numerical scheme using the 2000 random parameter sets represents the best performance that can be achieved by the scheme and is therefore of particular importance in model calibration. We tabulated the lowest MSEs for all the numerical schemes and the corresponding parameter sets in Table 3. Results showed that the lowest MSEs computed by the schemes were fairly similar except for RK1 and RK2 which produced slightly higher MSEs. Among the schemes, only RK3, RK4 and PIM with θ = 1/3 identified the same parameter set as that found by the analytical solution, whereas the other schemes pointed to different parameter sets. Despite the different parameter sets identified, the PIM with θ = 2/3 and IMP produced similar MSE minima as the analytical solutions, which indicates that model parameters may sometimes compensate for the inadequacy in numerical schemes. The computational time for 200,000 repeated Nash Cascade simulations with 10 linear reservoirs and 2000 Leaf River basin simulations with random parameter sets was summarized in Table 4 to assess the computational cost of all the numerical schemes. In both cases, the analytical solutions were remarkably computationally more expensive than all the numerical solutions due to its complex form. The computational time for the explicit Runge-Kutta schemes (RK1 to RK4) increased rapidly with the order of accuracy in the scheme because of the increasing number of points at which the integrated function must be evaluated. The partially implicit Euler schemes (PIMs with θ = 1/3 and 2/3 and IMP) were evidently lower than the high-order explicit Runge-Kutta schemes (RK3 and RK4) owing to an optimized algorithm for solving linear equation systems. Therefore, the computational cost of a numerical scheme also varies with its order of accuracy and implicit factor, which should also be considered in selecting numerical schemes for hydrological models.

Sensitivity Analysis
The results of 1DSA were plotted in Figure 5. Results showed that the MSEs of the simulated hydrographs computed by different numerical schemes were sensitive to the variations in α, N q and K q but less so to C par , b and K s . Because α, N q and K q were all associated with QFR, we may conclude that the parameters in QFR had a stronger impact on the output of HYMOD than those in SMA and SLFR. The MSE curves computed by RK1 (EXP) showed U-shaped distributions for α, N q and K q , whereas those computed by other schemes showed monotonic trends with increasing parameter values. For C par , b and K s , the patterns of the MSE curves were relatively similar. Because α, N q and K q were all associated with QFR while C par , b and K s were employed by SMT and SLFR, we may conclude that the effect of numerical schemes was more pronounced with the QFR parameters. The MSE response surfaces computed in the 2DSA for α and K q were plotted in Figure 6. The spatial patterns of the MSE response surfaces for different numerical schemes were similar, except for EXP (RK1). A high MSE region was located in the upper-right corner of the MSE response surface computed by EXP (RK1) but not in those computed by other numerical schemes. The high MSE region indicates that the numerical solution of EXP (RK1) significantly deviated from the analytical solution at large K q and α values due to excessive numerical error, whereas other numerical schemes could approximate the analytical solution over a wider parameter range. The influence of numerical schemes on simulated hydrographs can be seen more clearly in the spatial pattern of MSE difference between the numerical and analytical solutions, which were also plotted in Figure 6. The MSE difference in the upper-right corner of the response surface was effectively reduced by increasing the order of accuracy in the numerical scheme, although this might cause the MSE difference to increase slightly near the right edge of the response surface. The MSE difference could also be reduced by increasing the implicit factor. However, a second zone with significant MSE difference emerged in the central part of the response surface as the implicit factor increased, especially in the case computed by IMP, although this did not necessarily change the pattern of MSE response surface. In summary, low-order explicit schemes may distort the results of sensitivity analysis and this problem can be alleviated but not completely settled by increasing the order of accuracy or adding an implicit factor into the numerical schemes. Figure 6. Two-dimensional sensitivity analysis (2DSA) with different numerical schemes for α and K q . First row: mean square error (MSE) response surfaces for RK1 to RK4; second row: MSE difference between the analytical and numerical solutions for RK1 to RK4; third row: MSE response surfaces for the Euler schemes; fourth row: MSE difference between the analytical and numerical solutions for the Euler schemes.

The DSM Optimization
The optimal parameter sets identified by the numerical and analytical solutions and the MSE variation during each optimization are plotted in Figure 7. In all cases, the MSE values were stabilized at the end of the optimization, indicating that the optimization was successfully accomplished. The optimal parameters varied greatly in most cases, confirming the substantial impact of the numerical scheme on parameter optimization. To assess the overall performances of numerical schemes in the DSM optimization, we rated each optimization with regard to three different levels: if the difference between the optimized parameter computed by the numerical solution and the analytical solution was less than 10% of the parameter range, the optimization was classified as good; if the difference was larger than 10% but less than 30% of the parameter range, the optimization was classified as acceptable; otherwise it was classified as poor. The results were tabulated in Table 5 (G, A and P stand for "good," "acceptable" and "poor").
As shown in Table 5, the performances of the numerical schemes showed a strong dependence on initial parameter sets. Except for PIM with θ = 2/3 which showed good performances with all the initial parameter sets, all the numerical schemes provided acceptable or poor performances in some of the test cases. Both the order of accuracy and implicit factor were shown to influence the effectiveness of the DSM optimization. For the explicit Runge-Kutta schemes (RK1 to RK4), the scheme performance generally improved with increasing order of accuracy: RK1 identified the optimal parameter set in only one of the test cases, whereas RK2 to RK4 were successfully applied in a wider parameter range. For the partially implicit Euler schemes (EXP, PIMs with θ = 1/3 and 2/3 and IMP), the scheme performance first improved and then deteriorated with increasing implicit factor. PIM with θ = 2/3 provided the best overall performance among all the numerical schemes, which also demonstrates that the scheme performance did not vary monotonically with increasing implicit factor. In the optimization process, the MSEs for different numerical schemes converged at a similar speed, indicating that the schemes did not significantly affect the convergence rate of the optimization process. Except for EXP (RK1), all numerical schemes reached a similar MSE level at the end of the optimization, which illustrates how the model parameters compensated for the numerical error and affected the results of the DSM optimization.  Table 5. The scheme performances in the DSM optimization with different initial parameter sets. G, A and P stand for "good," "acceptable" and "poor," respectively.

Source of Numerical Error in HYMOD Simulations
Results of this study reveal a strong dependence of model output on numerical schemes. In the simulation of the Leaf River basin, both the order of accuracy and the implicit factor significantly affected the simulated flow hydrograph but in different ways: the implicit factor determined whether the flow hydrograph was overestimated or underestimated, whereas the difference between the numerical and analytical solutions decreased with increasing orders of accuracy.
The scheme dependence of the numerical results indicates the existence of numerical error, which is originated from the truncation error when using a numerical scheme to approximate the analytical solution to the storage equation of the linear reservoirs. Since the order of accuracy is a measure of the magnitude of truncation error pertaining to a numerical scheme, it is immediately seen that RK1 to RK4 were characterized by truncation errors of O(∆t 2 ) to O(∆t 5 ). Therefore, the numerical error in the flow hydrograph simulated by high-order Runge-Kutta schemes was smaller than that by low-order Runge-Kutta schemes. For partially implicit Euler schemes (Equation (18)), the truncation error can be expressed as: Equation (23) shows that the leading error term can change sign as the implicit factor, θ, in the partially implicit Euler schemes varies between 0 and 1, which causes the flow hydrograph to be overestimated and underestimated.
The Nash cascade process can be used to demonstrate how numerical error varies with numerical schemes. The outflow hydrographs with one to ten linear reservoirs were simulated by the explicit Runge-Kutta schemes (RK1 to RK4) and the partially implicit Euler schemes (EXP, PIMs with θ = 1/3 and 2/3 and IMP) and compared with the analytical solutions computed by Equations (9) and (19). Results (Figure 8) showed that the difference between the numerical and analytical solutions decreases with increasing order of accuracy and first decreases and then increases with increasing implicit factor. This is consistent with the simulation results of the Leaf River basin presented in Figures 2 and 3, which confirms that the source of numerical error in HYMOD simulations is the truncation error of the numerical schemes that were used in simulating the Nash cascade process. Since both QFR and SLFR are conceptualized as Nash cascade processes which differed only in the number of linear reservoirs being employed, numerical schemes may impose their impacts on both of these components. Euler scheme (IMP, θ = 1). ANA represents the analytical solution, θ represents the implicit factor and n is the number of linear reservoirs in the Nash cascade process. Curves from left to right in each plot correspond to n increasing from one to ten.

Implications for Hydrological Modeling
In previous studies, modelers commonly apply an ad-hoc numerical scheme without considering the adequacy of the numerical scheme. Our results demonstrate the need for implementing suitable numerical schemes in hydrological models. In the simulation of the Leaf River basin, different numerical schemes induced a scheme-dependent bias in the predictions of flow hydrograph (Figure 4). In addition, excessive numerical error induced by unsuitable numerical schemes may also distort the response sensitivity surfaces and bring unnecessary challenges to the gradient-based optimization methods. In 1DSA (Figure 5), the numerical schemes changed the monotonicity of the MSE curves for α and K q . In 2DSA (Figure 6), RK1 (EXP) produced a false high-MSE region in the α-K q response sensitivity surface. Kavetski and Clark [28] presented a case in which the response sensitivity surface was so strongly distorted that the optimal region could not be clearly defined. Therefore, the adequacy of the numerical schemes in hydrological models should be carefully considered.
Model parameters may sometimes compensate for the inadequacy in numerical schemes owing to the interaction between them. The simulation results with random parameter sets (Table 3) showed that all of the tested schemes, except RK1 (EXP) and RK2, generated similar flow hydrographs using different parameter sets. In other words, it is possible to obtain similar predictions using poor numerical schemes as those with "good" numerical schemes by tuning model parameters. However, this characteristic is not favorable in model development because structural identifiability is compromised by inappropriate implementation of numerical schemes. For example, the DSM optimization with IMP only provided acceptable or poor performances in all the test cases, which indicates that the best parameter sets identified by IMP is notably different from those identified by the analytical solutions. This will cause discrepancy in optimized model parameters and consequently the model predictions if we do not explicitly report the numerical schemes being used in the study.
The performance of a numerical scheme is closely associated with its numerical error, which should be controlled in obtaining the numerical solutions. The form of Equation (19) suggested that there are two possible solutions to reduce the truncation error of numerical schemes. The first solution is to reduce the size of the time step in the model simulation. However, this is not always feasible in hydrological studies, because the time step in performing a hydrological simulation needs to be compatible with the resolution of the input data. In this study, a daily time step is a natural choice because daily precipitation and evapotranspiration is provided as model input. To adopt a small time step in the simulation, the input data needs to be downscaled with interpolation techniques, which will bring additional source of error into the model output.
The other solution is to choose an appropriate numerical scheme for the hydrological model that faithfully reflects the property of the analytical solutions. However, it is impossible to test the performances of all the numerical schemes for hydrological model studies. Besides, the analytical solution is usually not available for most advanced hydrological models. To seek a suitable numerical scheme without assessing all possible candidates, one needs to understand the relationship between model performance and the common properties of numerical schemes, such as the order of accuracy and the implicit factor. Our results suggested that the use of lower-order explicit schemes should be avoided in modeling practice. This is consistent with the conclusions drawn by Clark and Kavetski [25], who disapproved the use of a fixed-step first-order explicit scheme. Clark and Kavetski [25] also recommended the use of IMP in parameter optimization, which produced smooth response surfaces in their sensitivity analysis. However, IMP did not identify the optimal parameter sets in two of the eight test cases in our study, suggesting that the low-order implicit schemes should also be used with caution in hydrological models. On the other hand, high order explicit Runge-Kutta schemes, such as RK3 and RK4, require higher computational cost than the partially implicit Euler schemes in computing the integral terms in simulating the Nash cascade process. Therefore, the partially implicit Euler schemes with the implicit factor ranging between 1/3 and 2/3 were recommended based on the HYMOD simulations presented in this study.

Conclusions
The importance of numerical schemes in hydrological models has not been well recognized in previous hydrological studies. In this study, we evaluated the effect of numerical schemes on hydrological simulations. HYMOD was employed with two types of numerical schemes (i.e., the explicit Runge-Kutta schemes and the partially implicit Euler schemes) to simulate the flow hydrograph of the Leaf River basin from 1948 to 1988. We did not limit ourselves to the comparison between specific numerical schemes but rather focused on finding a generalized relationship between model performance and the two common properties of numerical schemes (i.e., the order of accuracy and the implicit factor). The effects of the order of accuracy and the implicit factor on the model simulation, sensitivity analysis and parameter optimization were discussed. The following conclusions were drawn: (1) The results of the HYMOD simulations showed a strong dependence on numerical schemes.
Both the order of accuracy and the implicit factor significantly affected the simulated flow hydrograph but in different ways: the implicit factor determined whether the flow hydrograph was overestimated or underestimated, whereas the difference between the numerical and analytical solutions decreased with increasing orders of accuracy. The numerical error in the simulation is originated from the truncation error when using a numerical scheme to discretize the storage equation of the linear reservoirs. For the explicit Runge-Kutta schemes, the magnitude of numerical error decreases with increasing order of accuracy and thus the high-order Runge-Kutta schemes produced better prediction than the low-order Runge-Kutta schemes did. For the partially implicit Euler schemes, the leading error term may change sign as the implicit factor varies, which causes the flow hydrograph to be overestimated and underestimated. (2) Numerical error induced by unsuitable numerical schemes may distort the response sensitivity surfaces of hydrological models and cause discrepancy in optimized model parameters. Owing to the interaction between numerical schemes and model parameters, model parameters may sometimes compensate for the inadequacy in numerical schemes. However, this is not favorable in model development because the model structural identifiability is compromised by inappropriate implementation of numerical schemes. Therefore, the adequacy of the numerical schemes in hydrological models should be carefully considered. (3) Results of the numerical tests have important implications for selecting suitable numerical schemes for hydrological models. The performance of a numerical scheme is closely associated with its numerical error, which can be controlled either by reducing the time step or by choosing a numerical scheme that well approximates the analytical solutions. However, reducing the time step is not always feasible in hydrological studies due to issues in data availability. To design a new scheme for a hydrological model without assessing all possible candidates, one needs to understand the link between model performance and the common properties of numerical schemes, such as the order of accuracy and the implicit factor. The use of lower-order explicit schemes should be avoided and the low-order implicit schemes should also be used with caution in hydrological modeling studies. On the other hand, high order explicit Runge-Kutta schemes require higher computational cost than the partially implicit Euler schemes. Based on the HYMOD simulations presented in this study, we recommended the use of partially implicit Euler schemes with the implicit factor ranging between 1/3 and 2/3. This study demonstrates that the uncertainty in the output of hydrological models can be effectively reduced by choosing numerical schemes with proper orders of accuracy and implicit factors, which provides an optimized balance between model performance and computational cost.