Estimation of Annual Maximum and Minimum Flow Trends in a Data-Scarce Basin. Case Study of the Allipén River Watershed, Chile

Abstract

Data on historical extreme events provides information not only for water resources planning and management but also for the design of disaster-prevention measures. However, most basins around the globe lack long-term hydro-meteorological information to derive the trend of hydrological extremes. This study aims to investigate a method to estimate maximum and minimum flow trends in basins with limited streamflow records. To carry out this study, data from the Allipén River watershed (Chile), the Hydrologiska Byråns Vattenbalansavdelning (HBV) hydrological model at a daily time step, and an uncertainty analysis were used. Through a calibration using only five years of records, 21-year mean daily flow series were generated and the extreme values derived. To analyze the effect of the length of data availability, 2, 5, and 10 years of flows were eliminated from the analyses. The results show that in the case of 11 years of simulated flows, the annual maximum and minimum flow trends present greater uncertainty than in the cases of 16 and 19 years of simulated flows. Simulating 16 years, however, proved to properly simulate the observed long-term trends. Therefore, in data-scarce areas, the use of a hydrological model to simulate extreme mean daily flows and estimate long-term trends with at least 16 years of meteorological data could be a valid option.

1. Introduction

Water availability changes each year due to local, global, natural, and anthropogenic phenomena. These changes and the constant increase in demand, be it for human consumption, irrigation, hydropower, or industrial use, require increasingly efficient water resources planning and management [1]. Moreover, the frequency and magnitude of extreme meteorological events, such as floods and droughts are currently changing [2,3]. Despite scientific and technological progress, the population remains exposed to these events [4]; therefore, trend analysis of hydroclimatic variables has become the study focus for many researchers around the world [3,5,6].

Several studies have shown that changes in flow time series are due to climatic [7,8,9,10] and/or anthropogenic factors [11,12,13]. To analyze the hydrological behavior in a watershed it is important to determine maximum and minimum flow trends; however, ascertaining the history of a basin can be a difficult task when the available hydrometeorological information is limited [14]. Most drainage basins around the globe are ungauged or data-scarce, either because a required variable has not been sampled at the required resolution or because it has not been observed during a period of interest [15]. To address this limitation, different approaches have been developed to estimate stream flows in ungauged areas [16], for example, multiple linear regression (MLR), variations of autoregressive moving-average (ARMA) models, artificial neural networks (ANNs) [17,18], and the widely used conceptual or physically-based hydrological models.

Currently, hydrological models are important components in planning and management of water resources, since they allow the simulation of streamflow series through a simplified representation of hydrological processes. There are various hydrological models, which vary in complexity and utility [19]. Most hydrological models can be classified as conceptual models, with parameters, or at least some parameters, that cannot be directly physically interpreted or measured, making it necessary to estimate them through a calibration process [20,21]. These models can be lumped, semi-distributed, or distributed. Depending on the aim of the study and information availability, lumped models may be preferred, since they provide acceptable results in terms of accuracy [22] and are simpler to interpret and implement than distributed models, especially in scarce-data cases.

In recent years, several studies have reported a decrease in precipitation and in snow and glacier coverage [23,24,25,26,27,28,29] and increasing temperature trends [30,31] in south-central Chile. These changes might affect the way in which water is managed on a regional scale, as well as how extreme events, such as floods and droughts, are faced.

Several countries around the globe currently need to face the lack of hydrological data for water-related management and planning. In Chile, streamflow records have a relatively short length of 30 years on average [32], and gauging stations are unevenly distributed throughout the country [33], limiting basin-and regional-scale analysis of flow trends and making it difficult to study runoff changes [34,35].

Time series reconstruction using a dendroclimatological approach [36,37,38,39,40,41] has been one of the most used methods to conduct hydrological long-term trends analyses. Tree-ring records provide continuous series of past environmental changes for the last several centuries and in some cases, millennia [36]. However, there is still a lack of specific studies related to the analysis of extreme flow trends and, more importantly, such trends in basins with limited hydrological information. Several studies used extensive databases [42,43,44,45] and performed a data quality control considering criteria such as no gaps and a length of time-series of at least 40–50 years. Thus, the objective of this article is to investigate a method to estimate such hydrological extreme trends, with a controlled and known uncertainty, in a basin with limited data availability. With this aim, a combined approach of hydrological modeling with a calibration using 5-years of data, a short simulation window (up to 21 years), and uncertainty analysis of flow trends is used.

2. Study Area and Data

2.1. Study Area

The Allipén River watershed was chosen as a case study location. The Allipén River watershed until the Río Allipén en Los Laureles stream-gauge station (1652 km²) is located in southern Chile and monitored by the Chilean General Water Directorate (DGA; Figure 1). The Allipén River rises in the Andes Mountains. Its climate varies between warm, temperate and rainy in the middle and lower parts of the watershed to cold, temperate and rainy with Mediterranean influence in its upper part (the Andean area). Annual precipitation can reach 3000 mm and mean monthly temperatures oscillate between −3 and 18 °C. The watershed receives mostly pluvial precipitation, with slight snow influence in the upper part, presenting a mixed hydrological regime (pluvio-nival). In the upper watershed volcano-sedimentary rock formations of the Tertiary and Quaternary periods are prominent, while in the Central Valley there are unconsolidated deposits of glacial origin and highly permeable alluvial material [46].

2.2. Data

To implement the hydrological model of the watershed, mean daily precipitation, temperature, and monthly evapotranspiration series representative of the watershed are required. The Los Laureles, Quecheregua, Cunco, Cherquenco, Tricauco, and Malalcahuello rain gauge stations (see locations in Figure 1), managed by the DGA, were used. To estimate mean precipitation in the watershed, the inverse distance weighting (IDW) method was used. Due to the absence of temperature records, 0.25°—resolution (~25 km) gridded data series published by the Department of Civil and Environmental Engineering of Princeton University [47] were obtained and the Thiessen Polygon method was used to interpolate them. Based on these records, potential evapotranspiration was estimated using the Thornthwaite method [48].

Daily flow records from the Río Allipén en Los Laureles station (from 1990 to 2010) were used for the calibration-validation process and to carry out the analyses. To analyze the trend results, the available records from 1951 to 2010 were used.

As this study is focused on maximum and minimum flows, data gaps during the calibration and validation might affect the results. Therefore, prior to the calibration-validation stage, it was verified that no significant data gaps in the seasons of interest were present.

3. Hydrological Model and Methods

In this study, the Hydrologiska Byråns Vattenbalansavdelning (HBV) model [49] is used. This version of HBV is a free MATLAB-coded model that can be obtained from http://amir.eng.uci.edu/software.php. The HBV model and its various versions have been applied in a range of watersheds and climates and several countries [49], including Belgium [50], Sweden [51,52], China [53], and the United States [54], serving as a suitable alternative for flow estimation in a variety of watershed types. Moreover, it has been used at different time scales (e.g., Parra et al. [55]), where the obtained results support the model performance and its ability to adequately represent the main hydrological processes.

3.1. HBV Model Description

The HBV model [49] is a conceptual snow–rain water balance model that can be used as a lumped or semi-distributed model. In this study, the simplified version of the model presented by Aghakouchak and Habib [56] was used. This version simulates daily discharges based on daily precipitation, temperature, and potential evapotranspiration time series [57].

The model consists of three main modules (see conceptual diagram in Figure 2): a snow module, an effective precipitation and soil moisture module, and a response module [57]. The first module controls snow accumulation and melting. Precipitation accumulates as snow when the temperature is below a threshold value (TT); otherwise, the model considers precipitation rain and begins the snowmelt subroutine. The contribution of snowmelt to runoff is estimated through the simple degree-day method [49], depending on the difference between the actual and threshold temperatures.

The second module determines the contribution of precipitation to infiltration and surface runoff. First, it calculates daily potential evapotranspiration (ETPa) by reducing the monthly value based on a correction factor (C), the mean daily temperature and the long-term temperature, and potential evapotranspiration averages [57]. The permanent wilting point (PWP) is a fraction (LP) of field capacity (FC); when moisture surpasses the PWP value, actual evapotranspiration (ET) is equal to daily potential evapotranspiration. In addition, when soil moisture (SM) is below PWP, a linear reduction is applied to evapotranspiration (see upper left corner of Figure 2). Subsequently, the model calculates runoff (∆Q), which depends on precipitation (∆P), soil moisture (SM), field capacity (FC), and an empirical coefficient (β), which determines the relative contribution of rain or snowmelt to runoff (see upper right corner of Figure 2).

Finally, the response model estimates the flow based on two reservoirs, one above the other. The upper reservoir represents the flow near the surface, while the lower simulates baseflow (groundwater contribution), and both are connected through a percolation rate (k_p). There are three outlets, two in the upper reservoir (Q₀ and Q₁) and one in the lower (Q₂). When the water level in the upper deposit surpasses a threshold value (L), runoff is generated quickly in its upper part (Q₀). The responses of the other two outputs are relatively slower (Q₁ and Q₂), and recession coefficients k₀, k₁, and k₂ are used to ensure that Q₀ has the quickest response and that the response of Q₂ is slower than that of Q₁. For a more detailed description of the model consultation of Bergström [57], Aghakouchak and Habib [56], and Seibert [58] is recommended.

In the Andes of southern Chile, a precipitation enhancement of up to three times the precipitation recorded in the lower areas may occur. Therefore, the model also includes a precipitation adjustment parameter (A). This parameter allows the model to achieve a long-term water balance [59], and thereby correct the underestimation of precipitation as a result of the lack of records in the highest parts of each watershed.

Table 1. Model parameters and initial ranges used for the analysis.

Parameter	Description	Range
Mass balance
A	Precipitation modification parameter	0.8–2.5
Snow module
TT (°C)	Threshold temperature that indicates the initiation of snowmelt (normally 0 °C)	0
Cmelt $(\mathrm{mm}{}^{°}{\mathrm{C}}^{-1}{\mathrm{day}}^{-1})$	Fraction of snow that melts above the threshold temperature (TT) from the beginning of snowmelt.	0.5–7
${\mathrm{S}}_{\mathrm{f}}$	Snow accumulation adjustment factor	0.5–1.2
Moisture module
FC (mm)	Field capacity (storage in the soil layer)	0–2000
$\mathsf{\beta }$	Empirical coefficient that represents the soil moisture variation in the area	0–7
LP	Fraction of field capacity to calculate the permanent wilting point (PWP = LP*FC)	0.3–1
C (${}^{°}{\mathrm{C}}^{-1}$)	Correction factor for potential evapotranspiration	0.01–0.3
Response module
L (mm)	Threshold for quick runoff response	0–100
${\mathrm{k}}_0$ (${\mathrm{day}}^{-1}$)	Quick response coefficient (upper reservoir)	0.3–0.6
${\mathrm{k}}_1$ (${\mathrm{day}}^{-1}$)	Slow response coefficient (upper reservoir)	0.1–0.2
${\mathrm{k}}_2$ (${\mathrm{day}}^{-1}$)	Lower reservoir response coefficient	0.01–0.1
${\mathrm{k}}_{\mathrm{p}}$ (${\mathrm{day}}^{-1}$)	Maximum flow coefficient for percolation	0.01–0.1

presents a brief description of the parameters and initial ranges used, based on the studies of Aghakouchak and Habib [56] and Kollat et al. [60].

3.2. Model Calibration

To obtain a set of parameters that represent the watershed under study, a calibration process based on a regional sensitivity analysis (RSA) [61] was carried out using the Monte Carlo Analysis Toolbox (MCAT) [54]. In this study, 10,000 simulations were run at a daily time step using parameter sets within the valid range of each model parameter, extracted from a uniform distribution (Table 1). The number of simulations was chosen based on the study presented by Sarrazin et al. [62].

To allow a better understanding of the methodology, Figure 3 details the periods in which each process was carried out. Satisfactory results have been obtained using a 5-year or smaller time window for calibration [55]. Harlin [63] recommends using a calibration period between 2 and 6 years to find optimal parameters for the HBV model, and Brigode et al. [64] indicated that calibration periods longer than 3 years do not necessarily lead to more robust parameter sets. In this study, periods of 1, 5, and 5 years were used for model warm-up, calibration, and validation, respectively (grey, yellow, and red periods in Figure 3).

According to Seibert and Vis [65], one year is sufficient to warm up the model. Therefore, in the calibration and validation processes, the first year of analysis (2000 and 2006, respectively) was discarded to eliminate the influence of the starting values on the results.

3.2.1. Regional Sensitivity Analysis

The RSA method consists of generating a sample of N points (simulations) in the feasible space of each parameter, obtained from a previously defined distribution (generally uniform). The parameter sets are ranked from best to worst in terms of the chosen objective function. The ranked population is then divided into ten bins of equal size according to their objective values. The objective values associated with each parameter set are then converted to likelihood measures. The cumulative distribution of each group is then plotted, and the highest likelihood parameter distributions are indicated by bold black lines and the lowest likelihood distributions indicated by light gray lines (see Figure 4 as an example). In order to obtain the best representation of all the parameter sets (models), the calibrated parameters values were obtained according to the 50th percentile in the highest likelihood parameter distributions (bold black lines). Parameter sensitivity can be evaluated by assessing the separation of the ten curves. The greater the separation between curves, the greater the sensitivity of the model to the parameter under analysis. By contrast, the less separated the curves, the less sensitivity of the model to the parameter [54].

In order to independently represent the minimum and maximum mean daily flow trends, two calibrations were carried out based on the root mean square error (RMSE) and transformed root Mean square error (TRMSE) objective functions, which are focused on high and low flows, respectively [66,67]. Additionally, NSE values were calculated for both calibration and validation with the purpose of complementing the overall model performance evaluation (not only focusing on either maximum or minimum flows). The objective functions are described in the following sections (Section 3.2.2, Section 3.2.3 and Section 3.2.4). To validate the models, values within the limit of acceptability for the validation period were sought.

3.2.2. Root Mean Square Error (RMSE)

RMSE can be considered a multi-purpose criterion centered on the simulated hydrograph. This function is focused on high streamflow errors [66,67] and is calculated with Equation (1):

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{T}}{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\mathrm{Q}}_{\mathrm{s},\mathrm{t}}-{\mathrm{Q}}_{\mathrm{o},\mathrm{t}}\right)}^2},$$

(1)

where ${\mathrm{Q}}_{\mathrm{s},\mathrm{t}}{\mathrm{and}\mathrm{Q}}_{\mathrm{o},\mathrm{t}}$ are the simulated and observed stream flows in time step $\mathrm{t}$, respectively, and T is the length of the time series. Because RMSE measures differences between simulated and observed stream flows, it has the unit of measurement of the analyzed data (m³/s). An RMSE equal to zero indicates a perfect fit between the simulated and observed series, while the greater the RMSE value, the worse the fit between simulated and observed data. According to Singh et al. [68], RMSE values less than half of the standard deviation of the observed data can be considered low and indicate good model prediction.

3.2.3. Transformed Root Mean Square Error (TRMSE)

This function emphasizes low flow errors using the Box–Cox transformed [60] root mean square error, as shown in Equations (2) and (3):

$${\hat{\mathrm{Q}}}_{\mathrm{t}}=\frac{{\left(1+{\mathrm{Q}}_{\mathrm{t}}\right)}^{\mathsf{\lambda }}-1}{\mathsf{\lambda }},$$

(2)

$$\mathrm{TRMSE}=\sqrt{\frac{1}{\mathrm{T}}{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\hat{\mathrm{Q}}}_{\mathrm{s},\mathrm{t}}-{\hat{\mathrm{Q}}}_{\mathrm{o},\mathrm{t}}\right)}^2},$$

(3)

where ${\hat{\mathrm{Q}}}_{\mathrm{s},\mathrm{t}}$ and ${\hat{\mathrm{Q}}}_{\mathrm{o},\mathrm{t}}$ are the Box–Cox transformed simulated and observed streamflow at time step t, respectively, T is the length of the time series, and $\mathsf{\lambda }=3$ has a similar effect as a log transformation [60]. As the RMSE, we deemed less than half of the standard deviation of the transformed observed data to be the value that indicates good model prediction.

3.2.4. Nash–Sutcliffe Efficiency (NSE)

This index determines the relative magnitude of residual variance in comparison to the variance of the measured data [69]. The NSE is obtained with Equation (4):

$$\mathrm{NSE}=1-\frac{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\mathrm{Q}}_{\mathrm{s},\mathrm{t}}-{\mathrm{Q}}_{\mathrm{o},\mathrm{t}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\mathrm{Q}}_{\mathrm{o},\mathrm{t}}-{\overline{\mathrm{Q}}}_{\mathrm{o}}\right)}^2},$$

(4)

where ${\mathrm{Q}}_{\mathrm{s},\mathrm{t}},{\mathrm{Q}}_{\mathrm{o},\mathrm{t}}, \mathrm{and} {\overline{\mathrm{Q}}}_{\mathrm{o}}$ are the simulated stream flows in time step $\mathrm{t}$, observed stream flows in time step $\mathrm{t},$ and the average of the observed stream flows, respectively. Meanwhile, T is the length of the time series. NSE values vary between 1 and −∞, with 1 being the optimum value. Normally the value for a model to be considered suitable is NSE = 0.6 [55,70].

3.3. Trend Calculation and Analysis

After the calibration and validation processes, all daily stream flows were simulated, and the annual minimum and maximum values in the 1990–2010 period were extracted for simulations. To assess the effect of the length of the data series on the estimation of trends in the annual maximum and minimum stream flows, three situations were studied by eliminating 2, 5, and 10 years of simulated stream flows. All the possible trends according to the possible combinations (green period in Figure 3) were estimated. This methodology was adopted because information at basin scale in Chile generally does not consist of more than 30 years of records and may even consist of less than 10 years of data. In addition, since performing the entire procedure of randomly eliminating from 1 up to 19 years would significantly increase the computational resources needed, a simpler (reduced) approach was used, removing 2, 5, and 10 years from the analyses. Thus, totals of 210 (2 years removed), 20,349 (5 years removed), and 352,716 (10 years removed) possible combinations were analyzed. Then the trends of all the resulting time series (combinations) were calculated by the least squares method (green period in Figure 3). This is the most popular method used to determine the position of the trend line of a given time series. The trend line is technically called the best fit. In this method, a mathematical relationship is established between the time factor and the variable given [71]. The equation of the trend line can be expressed as

$$y=ax+b,$$

(5)

where $a$ represents the trend of the variable when $x=0$ and $b$ represents the slope of the trend line. If $b$ is positive, the trend line will be upward, and if $b$ is negative the trend line will be downward.

To analyze the influence of the length of each data set, the generalized likelihood uncertainty estimation (GLUE) method [72] was applied to the resultant estimated trends. Unlike the GLUE method, in which uncertainty bands are derived from the outputs that only produce adequate results based on a performance measure, in this study the median and the uncertainty bands between 5% and 95% of all the calculated trends were calculated. Finally, the results of each situation were compared with the long-term trend derived from observed streamflow data of the maximum and minimum flows over the entire period of records from the Río Allipén en Los Laureles station (1951–2010; blue period in Figure 3).

4. Results

4.1. Model Calibration and Validation

Figure 5 presents the results of the RSA in the calibration process focused on high and low flows (TRMSE and RMSE, respectively). In both cases, high sensitivity of the results to variations in parameters β and L is observed. β represents the contribution of precipitation and snowmelt to soil storage and runoff, and L is the limit for the quick response. Most stream flows in the watershed are of a pluvial origin; therefore, the model does not present sensitivity to parameters associated with the snow accumulation and snowmelt processes (module) (C_melt and S_f).

Along with the model sensitivity to β and L, in the calibration based on RMSE, it is observed that the parameter that most influences the model performance is FC, which represents field capacity. In addition, the results based on TRMSE indicate that the model is sensitive to K₂ (slow-response recession coefficient). The highest stream flows are a result of high precipitation in winter periods, when soil storage can approach its capacity; therefore, the process that takes on the greatest importance during high stream flows is soil moisture (FC). Meanwhile, in the recession and stable-streamflow periods (spring-summer), groundwater runoff and the low snowmelt input are represented in the model by the increase in sensitivity to the slow response (K₂).

Table 2. Objective function values and the reference performance rating in parenthesis.

Process	Period	Peak Flows		Low Flows
		RMSE	NSE	TRMSE	NSE
Calibration	2001–2005	45.1 (<54.8)	0.81 (>0.6)	0.89 (<2.0)	0.77 (>0.6)
Validation	2007–2010	39.4 (<40.1)	0.73 (>0.6)	0,91 (<1.1)	0.84 (>0.6)

shows the values of the objective functions calculated for the calibration and validation, while

Table 3. Parameter values obtained after calibration and validation focused on high and low flows (separately).

Parameter	Peak Flows	Low Flows
A	1.2	1.2
TT (°C)	0	0
$\mathrm{Cmelt}(\mathrm{mm}{}^{°}{\mathrm{C}}^{-1}{\mathrm{day}}^{-1})$	4.01	3.65
FC (mm)	1280	1030
$\mathrm{C}({}^{°}{\mathrm{C}}^{-1})$	0.15	0.18
$\mathsf{\beta }$	0.50	0.60
L (mm)	115	180
${\mathrm{k}}_0$$({\mathrm{day}}^{-1})$	0.44	0.44
${\mathrm{k}}_1$$({\mathrm{day}}^{-1})$	0.15	0.14
${\mathrm{k}}_2$$({\mathrm{day}}^{-1})$	0.015	0.008
${\mathrm{k}}_{\mathrm{p}}$$({\mathrm{day}}^{-1})$	0.11	0.12
LP	0.68	0.68
${\mathrm{S}}_{\mathrm{f}}$	0.83	0.85

shows the parameter sets obtained in each ease. The reference values of the objective functions are indicated in parentheses. For the RMSE and TRMSE values, Parra et al. [55], Singh et al. [58], and Moriasi et al. [60] suggest that the results must be less than half of the standard deviation of the observed and transformed data, respectively, and the NSE values must be greater than 0.6. According to the calculated indices, the parameter sets in both cases are suitable for simulating the stream flows of the Allipén River watershed.

In Table 3, it is observed that L is greater for low stream flows than for high stream flows and k₂ is lower for low stream flows than for high stream flows, which shows that the model calibrated with TRMSE emphasizes streamflow stability through the contribution of groundwater (slow process) to the total streamflow, and the model calibrated with RMSE is focused on quick responses (peaks) through the contribution of surface runoff to the total streamflow.

According to the sensitivity analysis, the objective functions, and the obtained parameter sets, the model is deemed suitable for use in the subsequent analyses. Figure 6 shows a comparison of the observed and simulated stream flows in the calibration and validation processes for each case.

4.2. Trend Uncertainty Analysis

Figure 7a–c shows the maximum daily streamflow trends, the median of the trends obtained by eliminating 2, 5, and 10 years of information, and the uncertainty bands for each case. The same is shown in Figure 7d–f for the annual minimum daily stream flows. It is observed that extreme events tend to increase in magnitude, as the maximum streamflow trend is positive, and the minimum streamflow trend is negative.

In Figure 7a,b,d,e it is observed that the medians of the trends are similar to the observed streamflow trend. Trends of the annual minimum daily stream flows present greater uncertainty than those of the annual maximum daily stream flows. As the uncertainty of trends calculated using 16 and 19 years is similar to that of the trend estimated from data, it proves suitable to use either 16 or 19 years of simulated flows to estimate the basin annual maximum and minimum flow trends. Meanwhile, for the trends estimated with 11 years of simulated data, uncertainty increases, and the slope of the median differs from the trend obtained from the streamflow records available in the watershed.

In south-central Chile, a decreasing precipitation and increasing temperature trend has been recorded [24,25,26,30]. Higher temperatures may lead to earlier snowmelt and to a shift of high flows from spring to winter. Earlier recession flows will be therefore expected with lower flows in summer, explaining the negative trend of the low flows observed in Figure 7.

5. Conclusions

The lack of stream flow records is one of the most common drawbacks when conducting hydrological analyses. Results of the estimation of annual minimum and maximum flow trends using hydrological modeling in watersheds with limited hydrological (meteorological and streamflow) information were presented in the present study. To estimate these trends, a 21-year period of observed and simulated stream flows was used. To account for data scarcity, trends were estimated in three situations (by eliminating 2, 5, and 10 years of simulated stream flows) and accounting for the possible combinations of the gaps.

According to the results of the trend analysis, minimum extreme flows tend to decrease, while maximum extremes tend to increase. The Allipén River presents a mixed hydrological regime where higher temperatures increase the snow melt, contributing to the increasing trend of high flows. Additionally, to derive the long-term trends of annual peak and low flows it is necessary to simulate at least 16 years of flows; otherwise, the resulting uncertainty is high, with levels at which even opposing extreme flow trends are obtained.

This study presented an approach to estimate flow trends in basins with limited stream-flow information. A minimum recommendable amount of data (16 year of data) is given in order to obtain results under controlled uncertainty. The presented method is not limited to the use of other hydrological models, objective functions, and/or sensitivity analysis. Despite the results presented here, more robust analyses can be obtained if basins with different characteristics and more extensive gaps of data are analyzed. Future research might complement the current research by including regional analyses and different trend estimation methods (e.g., parametric vs. nonparametric).