Uncertainty analysis on long-term runoff projection from the Budyko framework and a conceptual hydrological model

Runoff models applied to scenarios for current and future conditions are important for supporting decision-making in water management. However, hydrological processes are subjected to uncertainties related to model parameters, calibration process, model structure, and data measurements (Moges et al. 2021). Although there are several physically-based models available, the complexity of the physical system does not guarantee that those model types have superior performance (Orth et al. 2015); because the amount of information required to represent physical processes in complex models is ultimately related to an increase in model structure uncertainty (Orth et al. 2016).

The Budyko framework is an empirical model used in hydrology that was first developed for climatic characterization (Budyko 1974). Recently, this model have been used to predict water balances under climate change scenarios (Fernandez & Sayama 2015; Xing et al. 2018; Li et al. 2020; Melo et al. 2022) and other applications (Table 1). According to Berghuijs Gnann & Woods (2020), the uses of the Budyko framework often rely on untested assumptions about the conceptual relation between aridity and the partitioning of precipitation into streamflow and evapotranspiration. However, Budyko has provided a benchmark for the long-term water and energy balance with a range of applicability in hydrology (Table 1). Therefore, as simplifications are also desired in hydrological modeling, it may be further explored regarding conditions extrapolating the observed data range.

Fernandez & Sayama (2015) compared future projection from the Budyko framework with a global hydrological model and the results indicate that Budyko did not represent changes in soil moisture in the partitioning of precipitation into runoff. However, although conceptual models do represent more components of the runoff, these are also subject to uncertainties that should be quantified.

These conceptual hydrological models represent natural processes by simplified mathematical equations, where the physics of the system is represented, although there are empiricisms embedded in some parameters (Orth et al. 2015). These can be subjected to equifinality (Savenije 2010) and over-parameterization (Seibert et al. 2019), leading to predictivity uncertainties. However, according to Kwon et al. (2012), even though there are uncertainties related to parameters in conceptual hydrological models, they usually are much smaller than those of climate models.

An uncertainty analysis benefits modeling by identifying limitations in the model, guiding data collection, and quantifying the uncertainties regarding model predictions (Moges et al. 2021). Thus, uncertainty analysis in different methodologies of distinguished complexity to estimate runoff is essential to different scenario assessments.

Considering the importance of uncertainty analysis and the increasing demand for impact assessment, it is necessary to quantify uncertainties related to conceptual and empirical hydrological models to evaluate their applicability. To analyze uncertainties of the models in conditions that extrapolate calibrated data range, climate change scenarios provide representative conditions of climate characteristics differing from the observation range.

Therefore, the purpose of this research is to provide a detailed uncertainty assessment of long-term runoff projections in humid headwater basins for reference upon the reliability of model results. This analysis is also important for assessing how sensitive runoff estimations are under climate change scenarios due to the hydrological model structure. For this, runoff estimation uncertainties were analyzed in climate change scenarios considering the Budyko framework and a conceptual hydrological model in basins of a mountainous region, in southeastern Brazil, while also quantifying uncertainty related to the conceptual calibrated parameters.

The Upper Grande River Basin (UGRB) has a drainage area of 8,758 km² and flows into to the Paraná River Basin. The UGRB has great importance for energy production from hydroelectric plants and water supply in Brazil. The southern basin is characterized by stronger relief, with elevations ranging from 798 to 2,653 m.

The climate, according to the Köppen-Geider classification, is cwb in the Mantiqueira range (south) and cwa in the Campos das Vertentes region (north). Both climates have dry and cold winters and wet summers, with warmer summer temperatures in Cwa. Due to the orographic effect in the Mantiqueira Range, the average annual precipitation is approximately 1,700 mm in the south of the basin, while in the north it is 1,500 mm.

The main water use in the basin is related to energy generation in two hydropower plants: Camargos and Itutinga. Agriculture and livestock are the main economic activities. The land use and land cover most predominant in 2015 were pasture 36.5, forest (Atlantic Forest biome) 23.8%, and agriculture 23.9% (MapBiomas 2022). The soil type most abundant in the south is Inceptisols (57.6%) and in the north is Oxisols (20.2%). According to Araújo et al. (2018), the Inceptisols in the basin have physical characteristics that are adverse to plant development. Therefore, agriculture is more present in the north region.

In water resources, the Budyko framework has been used to represent the basin water balance by describing the mean-annual partitioning of precipitation (P) into actual evapotranspiration (AE) and streamflow (Q). This framework is described by the data-driven relationship between the evaporative index (EI = AE/P), and the ratio of the atmospheric water supply (precipitation) to water demand (potential evaporation, PE) which is defined as the aridity index (AI).

Basin characteristics are represented as a single parameter (w). Variations in parameter w has been previously related to land use changes (Donohue et al. 2010, 2007; Oliveira et al. 2019; Melo et al. 2022); climate variability (Jiang et al. 2015; Oliveira et al. 2019; Melo et al. 2022); terrain slope (Sun et al. 2014); and soil water infiltration and storage capacity (Yang et al. 2007).

Potential evapotranspiration (PE) was estimated by the FAO Penman–Monteith equation (Allen et al. 1998) from daily meteorological data. Actual evapotranspiration (AE) was calculated by water balance from 15-year moving average yearly runoff (Q) and precipitation data (AE = P – Q). From AE, PE, and P, Budyko parameters (w) from 15-year moving averages were estimated from Equation (1). To represent each basin as a single parameter, Budyko curve was interpolated considering 15-year average data from 1990 to 2010 by the non-linear least square method. This parameter was then used to estimate runoff in Equation (2). The estimated 15-year average runoff was analyzed against the observed data for model performance validation from 1997 to 2013.

The MHD-INPE is a conceptual distributed hydrological model used to simulate runoff at daily or sub-daily time steps. It has been applied in climate change studies in basin areas varying from 32 km² (Zákhia et al. 2022) to 1.5 million km² (Siqueira Júnior et al. 2015). This model simulates processes of evapotranspiration, vertical soil water balance, surface, subsurface, and groundwater flows, and propagation in channels. Further information on model structure is found in Rodriguez & Tomasella (2016).

In this study, the UGRB was represented by grid-cell areas of approximately 6 km². Soil and vegetation was represented by combining 30 m spatial resolution land use maps from Souza et al. (2020) and hydrological environments defined by the Height Above the Nearest Drainage (HAND) model (Rennó et al. 2008) using SRTM (EROS 2017) elevations. To take into account land use changes, the land use maps were updated every 10 years and HAND environments (Figure 1(d)) were discriminated as valleys (between 0 and 15 m), intermediate (between 16 and 80 m), and plateau (above 80 m). Soil and vegetation parameters can be found in the Supplementary material.

SPEA2 algorithm (Zitzler et al. 2001) was used for parameter optimization. Calibration was performed by changing the following parameters: soil layers depth (D1, D2, and D3) ranging from 0.1 to 50 m, multiplier of saturated hydraulic conductivity (fKss) ranging from 0.1 to 10 times the initial condition, horizontal transmissivity (TSub) ranging from 0.01 to 10,000 m² dia⁻¹, variation of SHC with soil depth (Mu) ranging from 1 to 3 m dia⁻¹ m⁻¹, soil anisotropy (Alpha) ranging from 0 to 10,000, minimum soil storage capacity (fCsi) ranging from 0 to 2%, delay time for superficial (Cs) and base (Cb) flows, both ranging from 0.01 to 1,000 s. The period from 1990 to 2003 was used for calibration, and from 2003 to 2013 for validation.

The objective functions selected for Pareto front solutions were the Nash–Sutcliffe efficiency (Nash & Sutcliffe 1970) calculated on daily simulated and observed data (NSE) and in their logarithm (lNSE). For SPEA2 optimization was attributed an archive and population size of 50 individuals, each representing a set of calibrated parameters in MHD-INPE. The algorithm ran 50 iterations. The initial values for the parameters were 10 m for soil layer depths (D1, D2, and D3), and 1 for the remaining.

Future climate projections were obtained by the Regional Circulation Model (RCM) Eta-Model with 20 km horizontal resolution. This model has simulated RCPs 4.5 and 8.5 scenarios (Riahi et al. 2011; Thomson et al. 2011) in South America (Chou et al. 2014a, 2014b) forced with four Earth System Models (ESM) as boundary conditions: MIROC5, HadGEM2-ES, CANESM2, and BESM. All these were evaluated in this study because in climate change studies a multi-model analysis is preferable to quantify uncertainty in future projections, as it considers differences in model structure and parametrization.

All climate variables were biased and corrected by the quantile-quantile method (Bárdossy & Pegram 2011) using data from 1980 to 2005 for historical reference. The RCP4.5 (intermediate) and RCP8.5 (pessimistic) scenarios were used for runoff estimations from 2006 to 2099.

Budyko framework and MHD-INPE simulations were used to estimate long-term runoff in three basins in the UGRB region: Grande River (GRA), Aiuruoca River (AIU), and Capivari River (CAP). Estimations from both methods were investigated qualitatively and quantitatively for uncertainty analysis evaluation.

Runoff projections uncertainties from both the MHD-INPE and Budyko framework were analyzed by 15-year averages. Runoff estimates from both methods using observed meteorological data were evaluated graphically against observed runoff. Performance measurements such as PBIAS, RSME, and the ratio of the RMSE and the average observed runoff were used to validate the runoff estimates.

On the other hand, the following parameters have small variations within the Pareto front members: depth of the first and third soil layers (D1 and D3), aquifer horizontal transmissivity (Tsub), soil anisotropy (Alpha), and delay time of the superficial flow (Cs), suggesting that these parameters may be sensible to runoff estimation. However, further analysis is necessary to determine parameter sensitivity.

When calculating the efficiency on the logarithm of the flow, the lNSE coefficient increases the sensitivity for the lowest flows, once the NSE is very sensitive to peaks (Moriasi et al. 2015). Thus, the results found for the simulated flow with better NSE favor the applications in flood analysis studies. As the goal of this study is to estimate long-term runoff, headwater basins were calibrated using the best lNSE simulation as it, overall, presented the smallest bias.

Considering the maximum flows at an excess frequency of less than 1%, higher errors were found in the AIU and CAP basins. In these basins, the observed maximum flow occurs with greater intensity than the simulated. However, it is worth mentioning that the maximum flow values measured for the rating curves were 239 and 178 m³ s⁻¹, while the highest values estimated from water level measurements were, respectively, 656 and 345 m³ s⁻¹. Extrapolation of the measured rating curve also occurs in the GRA basin. In addition, the rain gauge network is sparse, particularly at the top of the hills. Due to the orographic effect, higher precipitation occurs at higher altitudes and this is not captured by the existing precipitation gauge network. Therefore, the interpolated precipitation influences both the calibration and the performance of the model.

For the minimum flows above 95% exceedance frequency, small variations in simulated curves occurred in headwater basins. Simulations do not consider water consumption in the basins that may influence minimum flows, such as pumping of groundwater for use in irrigation, which may explain this difference. Data on the granting of water rights in the region, when well monitored, can indicate these consumptions and improve the representation of the watershed in the model.

Despite monitoring limitations in the headwater basins, the hydrological model showed a good performance region in terms of streamflow for the whole study basin. In general, the Pareto front, considering the objective functions NSE and lNSE, showed few variations between the solutions, with larger differences in parameter estimations in headwater sub-basins.

As the Budyko method uses observed runoff data of the period 1997–2003 to estimate the w parameter, in this period a circular inference occurs due to the mutual dependence between parameter and variable. Therefore, the validation period from 2004 to 2006 was considered for performance evaluation in Table 2, in which estimates of both the Budyko method and MHD-INPE simulations overestimated long-term runoff.

Table 2

Statistical performance of the 15-year average runoff estimated for the validation period (2004 to 2006) for GRA basin

Basin .	PBIAS (%) .	RMSE (mm) .	RSME/(%) .
Budyko	5.30	40.95	5.45
MHD-INPE higher NSE	7.20	54.62	7.27
MHD-INPE higher lNSE	5.30	40.46	5.42

Basin .	PBIAS (%) .	RMSE (mm) .	RSME/(%) .
Budyko	5.30	40.95	5.45
MHD-INPE higher NSE	7.20	54.62	7.27
MHD-INPE higher lNSE	5.30	40.46	5.42

View Large

The simulation in MHD-INPE with the highest lNSE showed lower errors in the GRA and CAP basins compared to the Budyko method while the simulation in MHD-INPE with higher NSE showed better performance in the AIU basin. The estimates with the Budyko methodology showed similar results to the simulation performed with the MHD-INPE with higher lNSE in the GRA and AIU basins. Only in the CAP basin, the Budyko method presented greater errors.

The MHD-INPE overestimation of the long-term runoff compared to runoff observed during the validation period was indicated by the PBIAS coefficient calculated in the daily runoff (Table S2, Supplementary material). This coefficient measures the trends in the simulated data relative to the observed data and measures whether the biases of these data are overestimated or underestimated (Moriasi et al. 2015). Therefore, it may indicate limitations in the estimation of long-term runoff.

Gupta et al. (2009) used the KGE coefficient as an objective function in the multiobjective calibration in a conceptual hydrological model and observed improvements in the bias and variability of the simulations. Therefore, in studies aiming to analyze long-term runoff, it is suggested to use functions such as KGE or PBIAS as an objective function in the SPEA2 algorithm in an attempt to identify solutions that better represent long-term runoff.

The Budyko method considers the variation of potential evapotranspiration in the estimation of runoff, which assumes there are no limitations to vegetation water loss, such as stomatal control or relative soil moisture, which are essential factors to estimate the actual evapotranspiration (Allen et al. 1998). Thus, the Budyko estimates generated by a fixed w parameter do not consider the variation of these factors or changes in land use and land cover.

The influences of climate and land use and land cover were related to the variation of the parameter w in several watersheds in previous studies (Jiang et al. 2015; Melo et al. 2022; Oliveira et al. 2022, 2019). For this study, when setting the parameter w, it is assumed that the watershed has stationary characteristics based only on the observed series. Melo et al. (2022) used the decomposition method of the Budyko curve in the same study area and observed that the parameter w is related to variations in the agricultural area and, in some basins, with the precipitation and maximum temperature. Thus, variations in parameter w due to changes in land use explained the variation in runoff of up to 70 mm in the CAP basin between 1985 and 2015. Therefore, the variation observed in Figure 6 during the validation period can be explained by this non-stationary characteristic of w and its effect on the estimated runoff.

On the other hand, the MHD-INPE represents the variations in land use with maps updated every 10 years, as well as represents water storage in the upper layer of the soil as a limiting factor for the evapotranspiration of plants. However, the performance of the hydrological model in the validation period also overestimates the observed runoff. There is, therefore, the possibility that the 10-year interval for updating the land use land cover map does not adequately represent the changes in the basin. Shorted intervals of land use maps should be evaluated to test this hypothesis.

Thus, between the years of 2011 and 2015 this region was affected by several droughts (Cunha et al. 2019; Tomasella et al. 2023), this could have influenced the capacity of the models to estimate the long-term runoff. In addition, giving the exceptional characteristic of the drought, observed discharge values had to be extrapolated since water stage values were beyond the measurements (ratting curve). Consequently, discharge values are also affected by uncertainties.

Besides uncertainties in input data, hydrological models with a large number of parameters are subjected to equifinality and overparameterization (Savenije 2010; Seibert et al. 2019), leading to uncertainties in the prediction of runoff. The Pareto front, on the other hand, supports the analysis of equifinality uncertainties because it represents a set of parameters that generate dominant performance considering multiple objective functions. As seen by the results presented in this study, the uncertainties regarding the calibrated parameters on the Pareto front for estimating long-term runoff do not cause considerable variations during the validation period.

However, in an uncertainty analysis study, Tibangayuka et al. (2022) observed that the sensitivity of model parameters may vary depending on the watershed. Therefore, further sensitivity or Monte Carlo analyses would be needed to identify and classify the sources of uncertainties presented in MHD-INPE simulations.

According to Seibert et al. (2019), in a model in which there is a wide variety of free parameters, there is a greater possibility of variations in the model behavior that do not represent the real characteristics of the watershed, especially in predictions outside the observation range. In this regard, less parameterized models can be beneficial for prediction studies and, therefore, it will be evaluated in the next item of this study.

Qualitatively, both methods (Budyko and MHD-INPE) presented similar trends and behaviors for estimating long-term runoff. In general, the Budyko estimates were higher than the MHD-INPE estimates under conditions with higher precipitation, which can be visualized by the maximum runoff values estimated by the multi-model ensemble (colored area in Figure 8). In the minimum values of the ensemble, it is observed that, at certain moments, the Budyko method estimates lower runoff than MHD-INPE. However, the estimated medians of the set showed few differences between the methods evaluated. These differences are quantified in Table 3, differentiating the Budyko estimate from the MHD-INPE.

Differences between scenarios occur mainly due to the RCP8.5 scenario estimating a greater reduction in average precipitation. Overall, the Budyko method estimates higher runoff for the maximum values of the multi-model ensemble in the RCP4.5 scenario. A lower runoff estimate was also observed for the minimum values in the multimodel ensemble in the RCP8.5 scenario.

Climate scenarios exhibit larger amplitudes of average precipitation than those conditions observed during the calibration period, which means that runoff estimates fall into a scenario not previously observed in the basin. The aridity index in the observational period was lower than 1.0 and the average rainfall was around 1,500 mm, while the climate scenarios, the AI reached values up to 3.0 and precipitation ranged from 600 to 2,000 mm.

The sensibility of watersheds to variations in w parameter should be considered when analyzing the results. As demonstrated in Melo et al. (2022), during the analyzed period of data the w parameter had varied up to 0.11 around the average value. The Budyko parameter w can be related to the variation of agriculture land coverage, average annual precipitation, and the long-term mean maximum temperature. They noticed that the variation in runoff due to these variables is, in general, inversely proportional to the variation explained only by the aridity index. Therefore, by not considering the response of vegetation to climatic conditions, the Budyko method might underestimate the actual evapotranspiration in the wettest years and overestimate it in the driest years.

Melo et al. (2022) also showed future long-term runoff and estimated variations in the parameter w with precipitation and maximum temperature and a greater amplitude in the estimated runoff was observed in the multimodel ensemble. This indicates that by extrapolating the climatic conditions in which variation of the w parameter was calibrated, greater uncertainties are added to the estimates. Thus, by representing the watershed by a fixed parameter w in this study, the uncertainties of this nature are minimized.

According to Mianabadi et al. (2020), Budyko assumes stationary conditions in the long-term water balance, but this does not occur in many basins, especially in arid and semi-arid regions. In the study of Wu et al. (2018) the errors in the estimation of the actual monthly and annual evapotranspiration by Budyko were explored, and the variation in soil water storage was observed as a dominant contribution to errors in more arid climatic conditions. Therefore, although this study considers annual averages of 15 years, the influence of soil water storage in scenarios with a higher aridity index should be better explored.

Uncertainties related to the estimates made by the conceptual hydrological model should also be considered. Despite the small variation estimated by different parameterizations with dominant performance during the calibration period, there are still uncertainties related to the structure of the model, the fixed parameters, and the data collected, among others. Especially considering a model where a wide variety of parameters are used, besides the calibrated ones.

Despite the differences raised, both methods estimated runoff with considerable similarities. Still, the median of the multi-model ensemble in both methods presented very similar estimates quantitatively. Therefore, despite the simplifications of the Budyko method, it seems to be a viable alternative for first-order analysis studies in places where there is limited information about the basin characteristics. This can be especially important for impact assessment in many watersheds all over Brazil, which lacks additional information regarding soil characteristics and precise land use land cover information.

The applicability of empirical models has also been recently explored in machine learning techniques. An advance of using empirical machine learning techniques is related to short-term predictions, which is not possible under the Budyko framework due to the water balance simplifications. A recent study indicated a better performance of machine learning methods in comparison to conceptual hydrological models for predicting extreme floods (Nearing et al. 2024). It is important to emphasize, however, that machine learning models make predictions based on the training acquired using past data, which can introduce uncertainties given the fact that global warming is driving the climate system to unprecedented conditions. Since the Budyko framework is based on simple physical considerations (i.e. the continuity equation) which are strictly valid under any climate conditions, the usefulness of the approach to constraint future predictions is undeniable, particularly in areas with scarce historical data and/or knowledge.

Among the parameters selected in the multi-objective optimization of the hydrological model, differences were observed in daily runoff extremes, showing uncertainties in the prediction of future extreme events or applications for disaster forecasting alerting and in the integrated management of water resources. However, for long-term estimations few uncertainties were observed between the simulations of the Pareto front using the multi-objective optimization in MHD-INPE during the validation period and in future scenarios.

In general, the future projections of long-term runoff by the Budyko method were quantitatively higher than those of the MHD-INPE under conditions with higher precipitation and lower aridity index. The opposite was observed in lower precipitation and higher aridity index. Thus, despite the quantitative differences in the estimation of long-term runoff, both methods qualitatively represented the same patterns in runoff variation. Overall, uncertainties related to climate models are greater than those of hydrological model complexity or parameters.

It should be noted that the Budyko curve offers wide possibilities of application in studies of hydrological impacts in hydrological basins. Its applicability to estimate future impacts in climate change scenarios offers an alternative to complex hydrological models in first-order analyses and/or large-scale applications. Therefore, this methodology can be better explored since it requires minimal information and is easily implemented.

This work was supported by Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG) (Grant number APQ-00709-21); Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001; Federal University of Lavras (UFLA) (Grant number 907957/2020); Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) (Grant Numbers 304695/2020-3 and 305295/2021-7).

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.