Quantification of Precipitation and Evapotranspiration Uncertainty in Rainfall-Runoff Modeling

Abstract

Mountainous watersheds have always been a challenge for modelers due to large variability and insufficient ground observations, which cause forcing data, model structure, and parameter uncertainty. This study employed Differential Evolution Adaptive Metropolis (DREAM) algorithm which utilizes Markov Chain Monte Carlo (MCMC) approach to account for forcing data uncertainty. A conceptual degree day snowmelt model, MIKE 11 NAM (Nedbor Afstromnings Model), was used to simulate snowmelt runoff from Ilgaz basin, with an area of 28.4 km² area, located in the northern part of Turkey. The mean elevation is around 1700 m and the basin is covered with broadleaf forest and has mainly brown soil with a high water holding capacity. Precipitation and evapotranspiration (ET) values were optimized in combination with model parameters conditioned on observed discharges and corrected values of input data were utilized for calibration and validation. Results showed that the observed precipitation was over-estimated by almost 10%, while evapotranspiration calculated by Penman–Monteith method was underestimated. The mean values of storm and ET multipliers were obtained as 1.14 and 0.84, respectively. When only parameter uncertainty was considered, calibration did not yield Nash–Sutcliffe Efficiency (NSE) greater than 0.64. However, when forcing data uncertainty was incorporated in the DREAM approach, an improved value of NSE (0.84) was obtained. After calibration and treatment of forcing data errors, the model yielded reasonable prediction uncertainty bounds and well-defined posterior distributions of NAM model parameters. Main objectives of the study are to assess the applicability of MIKE 11 NAM model to the selected catchment. In addition, the importance of errors in the input forcing variables to the model is demonstrated.

1. Introduction

A multitude of hydrological models have been developed to simulate physical processes of catchments. Despite the fact that several models maintain a decent sophistication in temporal and spatial complexity, they all transform the complex basin characteristics into naiver theoretical storages and transfer procedures. Characterization of catchment properties is usually accomplished by assessing the model parameters through calibration by tuning the model parameters such that the difference between the observed and measured responses becomes minimum. The parameters thus optimized efficiently represent the complex and distributed watershed properties. A key flaw of such calibration approach is the premise to ascribe all uncertainty sources to model parameters in an implicit manner. However, in addition to the uncertainty associated with model parameters, several other sources also affect model predictions including initial and boundary conditions, input data errors and structural deficiencies in model [1]. Therefore, it is not justified to implicate all sources of uncertainty to model parameters only considering the variable and complex nature of hydrologic system.

Lately, noteworthy studies have been conducted to develop a structured framework for uncertainty treatment [2]. Although a considerable number of studies laid emphasis on enumerating parameter uncertainty only [3,4,5,6,7,8,9,10], a few studies, however, also dealt with uncertainty related to input data [11,12,13,14,15,16], while other studies combined errors related to model structure as well [2,17]. Computational burden and time consumption are some of the issues related to incorporation of input data errors to be explicitly treated. However, proficient algorithms have been developed for complex distribution sampling [9,18] to approximate uncertainty in various sources and processes in modeling process.

Treatment of input uncertainty usually dealt with only rainfall uncertainty as this was considered as a major forcing factor for hydrological models [9]. However, one cannot disregard the influence of errors in other input sources, such as potential evapo-transpiration (ET). Evapotranspiration is a key factor of water accessibility in combination with rainfall [19,20,21]. Therefore, regional energy balance and water balance are highly linked with the unambiguous estimation of ET. The land surface models (LSMs) and satellite remote sensing along with the Penman-Monteith approach for the quantification of ET have uncertainties due to numerous factors [21,22,23,24]. In this study, the ET uncertainty was treated explicitly along with precipitation uncertainty and it is shown that it can prominently increase the efficiency of model outputs.

This paper uses an advanced Markov Chain Monte Carlo (MCMC) approach within a Bayesian framework for the assessment of posterior probability density functions of parameters. This approach is known as differential evolution adaptive Metropolis (DREAM) originally introduced by [16] and reintroduced as a full MATLAB package [25]. The DREAM algorithm simultaneously utilizes several chains for global investigation, while balance and direction of the proposal distribution are automatically tuned during the evolution process. The algorithm is adapted from a global optimization algorithm named as the Shuffled complex evolution Metropolis [9] which maintains the comprehensive stability and ergodicity. Ref. [26] have shown that DREAM can be successfully used for analyzing input rainfall errors, while dealing with additional latent variables for the rainfall events. They were able to decrease the uncertainty in model outputs with the explicit treatment of rainfall errors along with hydrological model parameters. This study extended their approach to incorporate the ET uncertainty by adding latent variables in the model hypothesis and to show the simultaneous effects on model output uncertainty.

The proposed scheme is based on the premise that observed discharge contains all the necessary information which can be extracted to do the backward hydrology [27]. Bayesian total error analysis (BATEA) approach of [14,28] forms the basis of this scheme. However, the quantification of rainfall, model parameters, and ET multipliers is accomplished by utilizing an unlike inference approach to delineate forcing data errors. The simultaneous treatment of rainfall and evapotranspiration uncertainties along with model parameter uncertainties would somehow make the calibration as resource extensive and high dimensional procedure. Therefore, an efficient iteration sampling scheme (MCMC) was used which will provide the constrained posterior distribution despite the multitude of parameters to be optimized. A conceptual rainfall runoff NAM model was used to quantify uncertainty in model simulations. NAM is a lumped, deterministic model which makes use of ten parameters for the simulation of runoff and has been exploited by many researchers [29,30,31,32].

In the following section, first, a short depiction of the conceptual NAM model and DREAM approach is outlined which is essential for the successful implementation of the methodology, followed by the approach used to incorporate rainfall and evapotranspiration multipliers. The approach extends the work of [12,26] and the application of DREAM to mitigate the input uncertainty for the model output of conceptual NAM model is demonstrated through various case studies. A summary is presented as a conclusion.

2. Materials and Methods

2.1. Hydrological System and Data Used

Ilgaz catchment (28.4 Km²) of Turkey is used for this study with a mean elevation 1700 m above the sea level. Snowmelt is the main cause of the peak flow in the catchment, which usually occurs in April. The annual average temperature for the study duration was 5.70 °C with a yearly snow depth accumulation of 36 cm. Likewise, yearly mean flow from the basin flow for the duration of study was 14 hm³. The rainfall in the study area is maximum in the months of June-July with average monthly rainfall reaching to more than 70 mm per day while November is considered as the driest month with average monthly rainfall of as low as 22 mm per day. Overall, the area receives a fair amount of rainfall throughout the year. Figure 1 depicts the site of the study area together with the flow measuring station, meteorological station, and elevation zones. The meteorological station used in the catchment lies at 33°43′05, 4104′15′ E at an elevation of 1775 m. Whereas Kastamonu station is located at 1000 m elevation above mean sea level and at 41°23′51″ N, 33°43′5″ E. The area comprises mainly of coniferous forests with a small share experiencing farming practices.

Daily discharge and rainfall data from 1 October 2012 to 30 June 2017 was utilized in the study (Figure 2). Climatological data were acquired from the weather station situated inside the catchment. However, some of the data were missing, which were consequently acquired from a nearby Kastamonu station, which has an elevation of 1000 m above sea level. Owing to the elevation difference between both the stations and to make us believe in using data from a lower elevation meteorological station, the Spearman rho statistical correlation and Mann–Kendall tests were employed, which exhibited good correlation to approximate the missing data. Daily evapotranspiration values were estimated by using Penman–Monteith method using R environment.

2.2. NAM Rainfall Runoff Model

A lumped, conceptual rainfall runoff code NAM (Nedbor Afstromnings Model) is utilized in this study. The water content from the overland flow, interflow and base flow is usually accounted by this model [33]. NAM simulates the lateral flow in the MIKE 11 module of Danish Hydraulic Institute (DHI), Denmark [1]. Figure 3 delineates the physical processes of the model that govern the runoff simulation and reflecting each sub-catchment as a single unit, hence the parameters and variables typically characterize mean values for the whole sub-catchments. A continuous runoff time series is achieved as an outcome of the model during the modeling period. Therefore, both base and peak flow states are yielded by the NAM model accounting for antecedent soil moisture conditions.

2.3. Markov Chain Monte Carlo Scheme

Markov chain Monte Carlo (MCMC) approach was originally presented by [34] to approximate the expectancy of any function f of a distribution π i.e., E_π f(x). Markov chain forms the basis of such techniques by exploring the search space through a random walk and successively visiting stable solutions [16]. MCMC algorithm executes experimental moves starting from the existing initial state of the chain to the next state z. Random Walk Metropolis (RWM) is known as the simplest kind of MCMC approach. First, a nominee point (z) is tried out based on the current position, x_t−1, and it is also symmetric i.e., q(x_t−1,z) = q(z,x_t−1). Additionally, the nominee point can be rejected or accepted on the basis of Metropolis acceptance probability Equation (1):

$$\mathsf{\alpha }\left({\mathrm{x}}_{\mathrm{t}-1},\mathrm{z}\right)=\{\begin{array}{l}\min \left[\frac{\mathsf{\pi }\left(\mathrm{z}\right)}{\mathsf{\pi }\left({\mathrm{x}}_{\mathrm{t}-1}\right)},1\right] \mathrm{if} \mathsf{\pi }\left({\mathrm{x}}_{\mathrm{t}-1}\right) > 0 \\ 1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \mathrm{if} \mathsf{\pi }\left({\mathrm{x}}_{\mathrm{t}-1}\right)=0\end{array}$$

(1)

where π (.) symbolizes the probability density function (pdf) of the aimed distribution. Accordingly, the chain travels and advances to z if the proposal is recognized, or it rests at its recent location x_t−1.

The proposal of RWM can jump both ways in equal probability and hence defined as the symmetric property of RWM. Equation (1) was further exploited by [35] to integrate asymmetrical proposal distribution as q(x_t−1,z) ≠ q(z,x_t−1), which prevents the converse leaping of the proposal. This addition is termed as Metropolis Hastings (MH) algorithm, and many prevailing MCMC sampling structures have their basis on MH algorithm. MH algorithm competence generally relies on the orientation and scale of the proposal distribution. For instance, too many candidate points could be discarded owing to the larger width of the proposal distribution, and as a result, slow convergence and ineffective mixing of the chain toward target distribution may occur. Conversely, a narrow distribution could accept all candidate points. Therefore, proposal distribution selection is acute and governs the efficacy of MCMC simulations [36].

Most of the problems relating to suitable jumping distribution can be handled by tuning the proposal distribution automatically. This strategy exploits the sampling history information to constantly adjusting the size and shape of the proposal distribution and advancing the sampler competently toward a restraining distribution. Examples of such methodologies are Adaptive Proposal (AP), Delayed Rejection Adaptive Metropolis (DRAM) algorithms, and Adaptive Metropolis (AM) [37,38,39]. The covariance adaptation strategy used in AM, AP, and DRAM is suitable for comparatively simple inference problems and does not work well when heavy tailed and multifaceted posteriors are included. Therefore, a single chain is incapable to deal with these complications. Owing to the above mentioned limitations of MCMC schemes, Differential Evolution-Markov Chain (DE-MC) method was introduced by [40] which runs N different Markov Chains simultaneously in parallel. DE-MC potentially resolves two main problems in MCMC sampling. Appropriate balance and direction of the proposal distribution are automatically selected and heavy-tailed and multimodal target distributions are efficiently accommodated.

2.4. Differential Evolution Adaptive Metropolis (DREAM)

Vrugt et al. 2009 demonstrated that an accustomed and self-adaptive evolution scheme with subspace randomized sampling could enhance efficiency of DE-MC scheme. This suggested scheme preserves a comprehensive stability and is known as Differential Evolution Adaptive Metropolis (DREAM). This method is an extension of the Shuffled Complex Evolution Metropolis algorithm [9] and preserves comprehensive stability and ergodicity while demonstrating outstanding competence on multifaceted, extremely non-linear, and high dimension distributions [18]. The detailed DREAM algorithm can be found in [16] and therefore is not included in this manuscript.

The DREAM methodology initiates with a preliminary population points to extract possible answers. Multiple regions of highest attraction are approached by using a number of various chains with unlike starting points. The search progress of individual chains is globally shared within the population of points. Therefore, the survivability of individual chains is improved and the generated MCMC sampler conducts a vigorous and competent pursuit of the parameter space. MCMC samplers are prone to the outlier chains, which are required to be excluded for improved convergence. Interquartile range statistic, IQR = Q₂ − Q₁, is calculated for this purpose and Q₁ and Q₂ symbolize the upper and lower quantiles of the N different chains, respectively. The removal of outliers affects the detailed balance of the process and therefore, it is usually accomplished during the burn-in period.

2.5. Inference Methodology

Imitation of the system under consideration by the assessed parameters is vital in catchment modeling. Most of the model parameters are approximated through calibration of catchment models. Figure 4 reveals numerous sources of uncertainty and depicts the common calibration methodology. The key notion in calibration is to relate simulated and observed outcomes and improve the parameter accuracy by curtailing the error between the outcomes. Preferably, the observed and simulated values dissimilarity should be as close to zero, which becomes difficult due to the existence of errors at various phases of the modeling scheme. Likely, the initial presumed parameter ranges might contain an intrinsic uncertainty and the uncertainty in the forcing data measurement as well. Likewise, the errors due to poor structure of model can cause uncertainty in the modeled results and eventually, the observed values or calibration data might be subjected to measurement discrepancies. The incongruity between simulated and observed outcomes is typically ascribed to the parameter ambiguity and the influence of forcing data errors is generally ignored which can potentially increase overall uncertainty of the model outputs.

2.6. Precipitation Forcing Data Error

A hydrological model, whether conceptual or physical, can be approximated as Equation (2)

Y = f (θ, ŋ, D,)

(2)

where Y = {y₁, …, y_n} is the simulated response of the model to forcing D (rainfall, ET, etc.,) and initial conditions ŋ. Function $f:D, \mathsf{ŋ}\to Y$. represents the model hypothesis which can be either deterministic or stochastic bounded by the parameter vector θ. The effect of initial conditions ŋ usually plummets as the time elapses after the start of simulation and therefore a warm up period is defined at the start of simulation to lessen the effect of state value initialization [16]. The suitability of function f is usually guaranteed by matching its output with observed measurements Ŷ = {ŷ₁, …, ŷ_n} and the residual vector is achieved by taking the difference between Y and Ŷ:

$${\epsilon }_i(\theta | \hat{\mathrm{Y}}, D, \mathsf{ŋ})=y_i (\theta | D, \mathsf{ŋ})-{\hat{y}}_i i=1, \dots , n$$

(3)

The error in Equation (3) should be principally close to zero but it cannot be materialized due to various error sources in the observed forcing data and initial conditions and due to the structural deficiencies in the model. The inappropriate selection of parameter vector $\theta$ also contributes to the total error. The sum of squares of residuals (SSR) is a commonly utilized measure to be minimized with simultaneous tuning of parameters without explicitly treating other sources of errors:

$$SSR(\theta |\hat{\mathrm{Y}},\mathbf{D},\mathsf{ŋ})={\Sigma }_{i=1}^n{\epsilon }_i(\theta |\hat{\mathrm{Y}},\mathbf{D},\mathsf{ŋ}{)}^2$$

(4)

Equation (4) can be optimized through numerical optimization methods, but these methods usually provide only the finest optimized parameter set without furnishing the requirement of the posterior probability density function of $\theta$, p$(\theta | \hat{\mathrm{Y}}, D, \mathsf{ŋ})$. Bayesian paradigm is a fruitful way to associate various probability distributions through Bayes theorem and once coupled with the Monte Carlo sampling, it can handle various error sources systematically. The suitability of this technique in quantifying the errors in environmental modeling has been shown by many studies [9,16,18,41,42,43]. Equation (3) can be written as a log-likelihood function after making some assumption [16] and the final form becomes as follows:

$$l(\theta | \hat{\mathrm{Y}}, D, \mathsf{ŋ})=-\frac{n}{2}\ln \left(2\pi \right)-\frac{n}{2}\ln \left({\sigma }_e^2\right)-\frac{1}{2}{\sigma }_e^{-2} X {{\displaystyle \sum }}_{i=1}^n{\left(y_i\left(\theta |D, \mathsf{ŋ}\right)-{\hat{y}}_i\right)}^2$$

(5)

Equation (5) does not consider the forcing data uncertainty explicitly. However, it can be modified to incorporate the precipitation uncertainty by a similar approach used by [26] who showed that the forcing error could be treated in hydrological modeling by allocating precipitation multipliers to each rainfall event rather than making every single precipitation value as a latent variable. So, first, the distinct rainfall events are recognized from the careful observation of hydrograph and hyetograph. In every event where k = 1, …, ϕ is allotted a dissimilar rain multiplier δ_k, and then integration is done in the vector of model parameters to be optimized. So, Equation (5) can be written in a modified form as shown in Equation (6):

$$l(\theta ,{\delta }_k| \hat{\mathrm{Y}}, D, \mathsf{ŋ})=-\frac{n}{2}\ln \left(2\pi \right)-\frac{n}{2}\ln \left({\sigma }_e^2\right)-\frac{1}{2}{\sigma }_e^{-2} X {{\displaystyle \sum }}_{i=1}^n{\left(y_i\left(\theta |D, \mathsf{ŋ}\right)-{\hat{y}}_i\right)}^2$$

(6)

2.7. Evapotranspiration Forcing Data Error

The above approach was extended in this study to treat the uncertainty in evapotranspiration data explicitly. To avoid over parametrization and complexity in the calibration scheme, ET values were considered constant for each month, although daily values were used in the hydrological model. This approach would confine the ET parameters to 12 per year which can be optimized using the prevailing algorithm. Each month j = 1, …, ζ was assigned a different multiplier ξ_j and the corresponding values were integrated in Equation (6) in a similar manner as rainfall multipliers were added. So, Equation (7) would now take the following form:

$$l(\theta ,{\delta }_k,{\mathsf{\xi }}_j| \hat{\mathrm{Y}}, D, \mathsf{ŋ})=-\frac{n}{2}\ln \left(2\pi \right)-\frac{n}{2}\ln \left({\sigma }_e^2\right)-\frac{1}{2}{\sigma }_e^{-2} X {{\displaystyle \sum }}_{i=1}^n{\left(y_i\left(\theta |D, \mathsf{ŋ}\right)-{\hat{y}}_i\right)}^2$$

(7)

2.8. Case Studies

To exemplify the robustness of the method, the NAM conceptual model was utilized for streamflow forecasting. NAM is usually developed with ten parameters which are described with their prior uncertainty ranges in

Table 1. Description and initial ranges of NAM rainfall runoff model.

Parameter	Range	Definition
umax (mm)	1–50	Surface storage water content (maximum)
lmax (mm)	1–1000	Root zone storage water content (maximum)
CQOF	0–1	runoff coefficient for overland flow
ckIF (h)	0.01–1000	Interflow time constant
ck1,2 (h)	3–100	Overland flow routing time constant
Tof	0–0.99	Overland flow threshold value for root zone
Tif	0–0.99	İnter flow threshold value for root zone
ckbf (h)	0.01–5000	Time constant for routing base flow
Tg	0–0.99	Groundwater recharge threshold value for root zone
Csnow (mm/d/°C)	0.5–10	Degree day coefficient

.

Keeping in view the addition of so many latent variables in terms of precipitation and evapotranspiration multipliers, we excluded some of the non-sensitive parameters of NAM model from the optimization algorithm. Based on the expert judgment and extensive sensitivity analysis results, which were carried out by the authors in their previous work, three parameters were not included, namely, root zone threshold value for overland flow (TOF), root zone threshold value for interflow (TIF), and root zone threshold value for groundwater recharge (TG). These parameters were set to a fixed value from their prior ranges.

We selected first 4 years of data for calibration and the remaining data were utilized for validation. The time series data in hand, albeit not in good quantity, provided quite reasonable results for a small catchment under consideration. In addition, using longer time series data would cause computational issues as with every added water year, and a greater number of multipliers would be incorporated into the calibration scheme. For the calibration period of 4 years, a total of ϕ = 170 storm events and ζ = 48 evapotranspiration months were distinguished. The prior uncertainty ranges or lower and upper bounds for the storm multipliers were adopted from [26] as [0.25 2.50] while the range for the evapotranspiration multipliers was set as [0.25 1.25]. The justification for the range for evapotranspiration multipliers will be discussed in the coming sections.

Table 2. Prior uncertainty ranges of the parameters for optimization.

Parameter	Range	Definition
Umax (mm)	1–50	Max W.C in the surface storage
Lmax (mm)	1–1000	Max W.C in root zone storage
CQOF	0–1	Overland flow runoff coefficient
CKIF (h)	0.01–1000	Time Constant for Interflow
CK1,2 (h)	3–100	Time constant for routing overland flow
δk k = 1, …, ϕ	0.25–2.5	Storm multipliers
CKBF (h)	0.01–5000	Time constant for routing base flow
ξj j = 1, …, ζ	0.25–1.25	Evapotranspiration multipliers
Csnow (mm/d/°C)	0.5–10	Degree day coefficient

presents the ranges of NAM model parameters and storm and ET multipliers to be optimized.

The prior distribution of parameters in Table 2 was set as uniform, as is the usual practice during iterative calibration schemes. However, ref. [28] raised objections for using the uniform prior for storm multipliers, as they believed that the resulting parameter estimation could not be well posed which is pertinent in inverse problems. But [26] did not seem to second this surmise while treating input rainfall uncertainty using the DREAM methodology and showed that DREAM could reasonably handle the less informative prior for the rainfall multipliers. Therefore, the prior uniform distribution was used both for storm and evapotranspiration multipliers in this study.

In the first stage, NAM parameters were estimated without incorporating forcing data uncertainty to set a benchmark for our further scenarios (case 1). Then, the estimation of NAM parameters along with precipitation multipliers is done (case 2) and finally evapotranspiration multipliers were estimated (case 3) to see the impact on model efficiency and corresponding uncertainty ranges.

3. Results and Discussion

3.1. Classical DREAM Simulation (Estimation of NAM Parameters Only)

In this case, the NAM parameters (modified) were assessed explicitly. Marginal posterior probability distributions of the estimated NAM parameters are shown in Figure 5. Almost 500,000 simulations were generated and the final posterior distributions were obtained utilizing the last 25 percent of the samples after reaching a convergence threshold.

Prevalent asymmetry can be seen in the rightly skewed marginal probability distribution of most of the parameters. However, CKBF and C_snow parameters posterior distribution is relatively normal depicting the mean value close to the middle. The presence of other sources of inherent uncertainties is evident from the deviance of Lmax, CQOF, Umax, CK12, and CKIF from the normal distribution and therefore, instigates to explicitly judge those sources, which might be because of model structure or forcing data. Yet, a prominent aspect of these posterior distributions is the occupancy of shorter regions, irrespective of the prior ranges.

95% prediction uncertainty band are shown in Figure 6 portraying the rendering of model parameter uncertainty into prediction uncertainty. Explicit representation of validation and calibration periods are shown in Figure 6a,b. Even though the simulated uncertainty bands shadowed a matching pattern with the observed stream discharges (calibration period), some large incompatible portions were still present and it looked like the uncertainty bands could not hold maximum observation points. Unexpectedly, the patterns were highly incompatible during the validation period with a broader shaded area, demonstrating the inefficiency of assessed parameters in the validation period. The shaded portion was comparatively thinner during the calibrating stage, but it failed to capture the maximum flows most of the time and the model results were underestimating the observed flow. These results triggered the intuition of having further error causes that could reinforce the surmise of handling the input forcing data errors in an explicit manner.

3.2. Estimation of Storm Multipliers along with NAM Parameters

The simultaneous estimation of storm multipliers and NAM parameters was done to explicitly inculcate precipitation uncertainty. For the calibration period, a total of δ = 170 precipitation multipliers were optimized from the selected range [0.25 2.5] using Latin hypercube sampling. In this case, DREAM was prepared to optimize d = 177 parameters as δ = 170 precipitation multipliers were fused with seven NAM parameters. A population size of N = 2d was set, based on the literature [25]. For this case, 500,000 simulation runs were performed and the last 20% simulations were used after DREAM converged to the threshold limit. For this case, the evapotranspiration values were kept constant for the whole calibration period just to observe the sensitivity of evapotranspiration for Ilgaz catchment. Figure 7 reveals the marginal posterior density graphs with explicit treatment of precipitation uncertainty.

The marginal pdfs of the parameters were well defined and most of the parameters sustained quite normal distribution, except Umax, which was still in the truncated form with most of the density toward the left side. In addition, after explicitly treating the precipitation uncertainty, the mode of parameters had also shifted to different values and it was evident that the estimated value for the overland flow runoff coefficient CQOF was now under physically viable range, which was previously very low. The values for the degree-day coefficient Csnow were now properly distributed around the mean and the mode was in a reasonable range. Interestingly, the mode of Umax remained almost the same toward the lower value and it seemed that it was still hungry for the extension of the range toward the lower side of the hypercube that cannot be further relaxed without harming the rather realistic behavior of the model.

Similarly, the 95% prediction uncertainties are shown in Figure 8a,b. It is observed that the exclusive treatment of precipitation errors not only compensated for the mismatches between the peaks, but it also provided quite reasonable and narrower 95% prediction bounds (Figure 8a). The peaks that were dominantly underestimated in Figure 6a, were somehow in good agreement, although the peak for the water year 2015 still could not be reached. Moreover, there were also some small portions where the model was showing underestimation. The validation period did not conform to the observed stream flows and the 95% prediction intervals were wider with less coverage of the observed data. These discrepancies implied that there was still some margin that other sources of errors might still be contributing toward the model prediction uncertainty, which could be treated such as evapotranspiration and model structure error.

Figure 9 presents the histograms of all the precipitation multipliers combined. It is evident that the observed precipitation was largely underestimated, as the mean value of precipitation multipliers was approximately 1.37. This means that the observed precipitation was significantly underestimating the actual precipitation by almost 37%. As we have set the evapotranspiration values as constant for the whole calibration period just to check the sensitivity of model to evapotranspiration, the mean value of 1.37 might be for the compensation for that fact. Figure 6a also supports the evidence of the existence of some considerable errors in the input forcing, investigation of which deemed unavoidable.

For stream flow prediction, the similar concept can be deployed to treat the precipitation uncertainty for the validation period. This can be materialized by creating the ensemble of precipitation values for the validation period. From the precipitation multiplier histogram, a single multiplier for individual storms would be sampled, followed by combining this vector of multiplier with observed rainfall data. In this way, an ensemble of precipitation values would be achieved which can be grouped with the already estimated model parameter values and stream flow hydrographs outside the calibration period can be realized. The approximated stream flow hydrographs include an explicit treatment of precipitation data error.

3.3. Explicit Treatment of Evapotranspiration Uncertainty

In this section, the results of simulations with an explicit treatment of forcing data (precipitation + evapotranspiration) error are discussed. Equation (9) was used to explicitly incorporate evapotranspiration uncertainty along with precipitation and model parameter uncertainty. The evapotranspiration data error mainly affects the dry months’ period and it will be rather easier if evapotranspiration multipliers for those months are estimated. However, it was decided to include the whole water year instead of using the dry season only. The justification for this approach could be the weakness of the model simulation for the low flow seasons and since evapotranspiration is usually governed by many other factors than sunshine itself.

Figure 10 presents the posterior parameter distributions of NAM parameters. The simultaneous estimation of precipitation and evapotranspiration multipliers did not affect the posterior pdfs of NAM parameters, except that the mode of overland flow runoff coefficient CQOF increased and shifted closer to 0.5. Moreover, the maximum water content in the surface zone storage shifted to a higher value with overall high probability density region (HPD) in a very small range.

The explicit treatment of evapotranspiration uncertainty resulted in over-estimation in the measured evapotranspiration values, which surpassed the actual evapotranspiration by almost 14%. The mean value for evapotranspiration turned out to be 0.86 with some of the multipliers still in the close vicinity of 01. This can be seen in the combined histogram figures of evapotranspiration multipliers (Figure 11). Although the incorporation of evapotranspiration multipliers in the monthly evapotranspiration values did not cause a significant effect on the NAM parameters distribution and streamflow prediction, however, the uncertainty in the precipitation values was considerably lessened with a new mean value of 1.14 (Figure 12). It means that the explicit treatment of evapotranspiration data largely compensated for the uncertainty in the precipitation data errors for the Ilgaz catchment that supported our argument that both input sources should be given an explicit treatment.

The translation of parameter uncertainty toward stream flow prediction can be depicted from Figure 13 which was simply 95% prediction uncertainty range for calibration and validation periods. It seemed that the explicit treatment of ET just compensated for the over estimation of precipitation values, thus improving the output hydrograph. However, it still improved certain underestimation for simulated discharges. For instance, in April 2013, the peaks were a bit matching, while before, there was an abrupt drop. The gap in the month of January 2015 reduced considerably. The same was the case with August 2015 where the peaks tended to match. However, the 2015 peak (April) was not reached by the model. The mismatch for the month of May 2016 was also compensated for.

Finally,

Table 3. Summary statistics during calibration and validation periods for different scenarios *.

	Calibration Period			Validation Period
	NSE	RMSE	Correlation	NSE	RMSE	Correlation
Case Study 1	0.61	0.32	0.7	0.45	0.4	0.5
Case Study 2	0.778	0.22	0.88	0.70	0.24	0.8
Case Study 3	0.84	0.18	0.93	0.77	0.20	0.83

* Case study 1: classical calibration problem without explicit treatment of forcing data uncertainty. Case study 2: simultaneous estimation of parameter and precipitation data uncertainty. Case study 3: simultaneous estimation of parameter, precipitation, and evapotranspiration data uncertainty.

shows the summary statistics of all the case studies considered separately for calibration and validation periods. The tabulated values correspond to the mean discharge values from the ensemble of outputs obtained from the posterior pdfs derived through DREAM using 4 years of calibration period. For the evaluation periods, stream flow was simulated by generating a precipitation ensemble for each storm by utilizing various storm multipliers drawn from the marginal posterior distribution graphs. The acquired precipitation records were combined with the NAM parameters posterior pdfs (obtained through calibration) to be subsequently used for prediction. The values in Table 3 revealed that case study 3 could achieve the best performance when precipitation and evapotranspiration uncertainty were treated explicitly. Case study 1 showed the least performance, even though an efficient iterative method (DREAM) was utilized for the calibration purpose. This led to the need for extra measures for the explicit treatment of other sources of errors, such as input or model structure uncertainty. In this study, the focus was given to input data uncertainty and precipitation and evapotranspiration errors were treated simultaneously with the NAM parameters. The results were appealing, as for the calibration (case study 2), the NSE value escalated to 0.778 from 0.6 falling in the acceptable region. The RMSE value also dropped significantly from 0.32 to 0.22. The results for the validation period for case study 1 were not at par with quite a disagreement between observed and simulated discharges and a wider range of 95% prediction uncertainty intervals. This was compensated to some extent in the validation results for case study 2, where the NSE approached 0.7 from 0.45 while a significant decrease in RMSE value was also observed. Lastly, for case study 3, the efficiency measures were improved slightly, especially for the validation period as the integration of evapotranspiration error largely compensated for precipitation uncertainty.

Table 4. Posterior maximum probable (MAP) and mean values of modified NAM parameters, precipitation multipliers, and evapotranspiration multipliers for different scenarios *.

	Case Study 1		Case Study 2		Case Study 3
	MAP	MEAN	MAP	MEAN	MAP	MEAN
Umax	5	5.2	5.05	5.1	10.2	10.3
Lmax	1.9	2.1	1.76	1.81	1.76	0.81
CQOF	0.01	0.01	0.14	0.14	0.48	0.45
CKIF	145	156	70	75	70	75
CK12	12	12.5	20	23	21	24
CKBF	420	425	685	692	685	691
CSNOW	1.6	1.62	3.1	3.16	3.10	3.16
δ	-	-	-	1.37	-	1.14
ξ	-	-	-	-	-	0.86

* Case study 1: classical calibration problem without explicit treatment of forcing data uncertainty. Case study 2: simultaneous estimation of parameter and precipitation data uncertainty. Case study 3: simultaneous estimation of parameter, precipitation, and evapotranspiration data uncertainty.

contains the optimized and mean values of the modified NAM parameters, precipitation multipliers, and evapotranspiration multipliers. The table shows that the inclusion of precipitation and evapotranspiration multipliers changed the mode for most of the parameters while uncertainty range of the some of the NAM parameters also increased.

4. Conclusions

The estimation of model parameters and structural or input forcing errors require a proficient and robust algorithm, which can tune the posterior probability density function of the parameters. The MCMC algorithm used in this paper can competently estimate the posterior pdf from a high-dimensional and complex parameter space. The MATLAB package DREAM used in this study utilizes MCMC approach, runs simultaneous chains, and automatically tunes the scale and orientation of the proposal distribution during the sampling procedure. The DREAM has the capability to maintain the detailed balance and ergodicity and better depicts the efficiency in non-linear, multimodal, and complex distributions. The applicability of DREAM was projected by application to a small catchment, which has the contribution from snow melt runoff too. Therefore, the usefulness of the method for the catchments having inherent complexity and un-evenness is also demonstrated.

A conceptual rainfall runoff model, containing ten parameters, is used for streamflow prediction. To tackle over-parametrization and to reduce the computational burden, some of the parameters were set to the fixed values through extensive sensitivity analysis. As a major objective of this study, rainfall and evapotranspiration uncertainty was treated explicitly which significantly tuned the posterior distribution of model parameters. It not only helped in finding the optimized parameters of the model, but also highlighted the percentage of error in the forcing data.

Three case studies were considered in this study. First, the conceptual model parameters were calibrated only followed by combined rainfall and model parameters and combined rainfall, evapotranspiration, and model parameters. These case studies were set to set a benchmark for prediction results and the posterior distribution of parameters. Seven parameters were used in the case study 1 while in the second case, 170 rainfall multipliers were detected and incorporated in the model. Next, the evapotranspiration multipliers (12/year) were also added in the modeling hypothesis to materialize case study 3. For case study 2, rainfall showed an underestimation by almost 37%, which was a significant number. However, the treatment of rainfall uncertainty tuned the posterior distribution and lessened the width of uncertainty bounds for streamflow prediction. The inclusion of evapotranspiration multipliers did not have a significant effect on model parameters. However, the rainfall uncertainty was minimized and the mean value for the rainfall multipliers was reduced to 1.14. The evapotranspiration error treatment likewise showed that the measured values of ET were over estimated by 14%.

Although the inference methodology could give reasonable results and reduce the errors due to forcing data uncertainty, the improvement in the methodology still was a question. Keeping in mind the computational efficiency, the evapotranspiration multipliers were set to 12/year which could be further refined. The number of simulation runs may also affect the posterior distributions of parameters. The present study provided satisfactory results for the catchment under consideration, but the regionalization of the study approach is highly required. The main limitation of the study include the curtailed amount of time series data and a small study region. Nevertheless, the model still performed well in both calibration and validation showing the efficacy of MIKE 11 NAM model over the study area. Moreover, the applied methodology in this research could also be performed for other geographical regions in the world in order to make the work more generalized. It can offer a useful insight to the water resources managers about the errors present in the forcing data and allow them to make substantial measures to effectively collect and acquire the climatic inputs for hydrological models while minimizing the inherent errors.