Study área

The high Cauca basin is located in southwestern Colombia. The Cauca river follows a south-north direction in an inter-andean valley (Figure 1). In the first 153 km, the channel descends 2.8 km to reach Salvajina reservoir, whose drainage area is 3,652 km². Downstream, the basin reaches 5,111 km² at La Balsa station and finally the area of interest closes at the Juanchito hydrometric station, covering a total of 8,556 km². Both stations were selected considering: 1) the proximity to the reservoir outlet, 2) Juanchito station is the control point for the volume operation in the Salvajina dam and 3) Juanchito is located on the outskirts of the city of Cali, which is the third most populous city in Colombia with 2.4 million inhabitants.

Figure 1 Location of the high Cauca basin and the hydrological network (left). Monthly rainfall in the area tributary to the reservoir and the tributary flow rate to Salvajina (right and above). Monthly storage in the Salvajina reservoir and average flow rate at Juanchito regulation target station (right and below). 

The Salvajina reservoir has a maximum storage capacity of 848 hm³, its operation objectives are flood and estuarine control, guaranteeing flow rates between 900 m³ / s and 130 m³ / s for 95% of the time at Juanchito station (^{Sandoval, Ramirez, & Santacruz, 2011}).

The hydrological regime in the region is strongly affected by the climatic variability associated with the double passage of the Intertropical Convergence Zone (ITCZ), orography, the processes that occur in the Atlantic ocean, the Caribbean sea, the Pacific and the Amazon, but above all, by the influence of the ENSO phenomenon in its extreme phases (^{García, Botero, Bernal, Ardila, & Piñeros, 2012}; ^{Puertas & Carvajal, 2008}; ^{Rueda & Poveda, 2006}).

Mean annual rainfall in the upper Cauca river basin is 1900 mm, and a bimodal cycle with greater rainfall predominates in the periods: March-April-May (MAM) and September-October-November (SON) (Figure 1, right) (^{Sandoval & Ramirez, 2007}). Strong connections have been reported between mean monthly flow rates of various rivers in western Colombia with the El Niño Oceanic Index - 𝑂𝑁𝐼, and other signs of change in the surface temperature of the Pacific ocean, the Pacific Decadal Oscillation - PDO, with composite climatic signals for atmospheric and oceanic variables such as the ENSO - MEI Multivariate Index and exclusively atmospheric variables such as the Chorro del Chocó - CCC and the SOI Southern Oscillation Index, among others (^{Jiménez-Cisneros et al., 2014}; ^{Poveda, Jaramillo, & Vallejo, 2014}; ^{Poveda et al., 2006}).

Reservoir index

^{López and Francés (2013)} propose an IE reservoir index (Equation 1), which establishes the degree of alteration of the hydrological regime based on the percentage of unregulated area, and the percentage of average runoff that cannot be retained in the reservoir.

Where, ARit is the tributary area to the reservoir in km², in year t, AS is the area tributary to the station in km², CS is the mean volume drained in the basin to station S in the year t and in hm³, CR is the total storage capacity of the reservoir in hm³ in year t; and N the number of reservoirs upstream of the station. For a basin with a single reservoir at the head and without other significant anthropic changes over time (eg transfer, increase / decrease of reservoirs, etc.) that affect the variables ARit y CRit, this Index works as a discrete variable or change signal that is intense / weak depending on the proximity of the capacity station of interest to the dam, and has been used in statistical non-stationary flood modeling (^{Liang, Jing, Wang, Binquan, & Zhao, 2017}; ^{Machado et al., 2015}).

The impacts of a reservoir on the rolling of floods depends upon: the available capacity of the reservoir, the size of the works in relation to basin water supply, the uses of stored water and water levels in the reservoir before the channel, etc. (^{Moreno, Begueria, Garcés, & García, 2003}). Available capacity means the free volume or difference between the maximum storage capacity and the volume of water at any given time. Changes in available capacity may limit / favor the management of maximum annual flow rates at the target station. Therefore, in this paper it is proposed to adjust equation 1, replacing the relationship between reservoir capacity and average runoff, by the annual reserve volume coefficient D(t)=CRt-Vat/CRt) that determines the proportion of volume available versus the capacity of the reservoir, so that the index IE2 is:

Where Vat is the mean annual volume stored in the Salvajina reservoir, estimated with data from bathymetry conducted by the CVC in 1985, 2003 and 2011, and the reservoir end-of-day level records. With this change, in addition to the effect of a decrease in regulatory capacity as the tributary area at the capacity station increases (established by the ratio of areas), the index considers volume management as a condition that influences observed maximum annual flow rates.

Indices of climatic variability

Bearing in mind that the ENSO phases are defined by the US National Ocean and Atmosphere Agency (NOAA) based on the El Niño Oceanic Index (ONI), in this work the mean annual value of ONI is adopted as a criterion to classify each hydrological year, as follows: one year is considered La Niña if the ONI mean annual value ≤ -0.30; El Niño years occur when ONI ≥ 0.30; and a year is defined as Normal, provided that -0.30 < ONI < 0.30. The low frequency climatic indices used in the modeling correspond to the mean annual values of the El Niño Oceanic Index - ONI or quarterly moving mean of the anomalies of Pacific surface temperature in the Niño3-4 region and the ENSO - MEI Multivariate Index, which corresponds to the First Principal Component of a set of oceanic and atmospheric signals throughout the Tropical Pacific region. All values are estimated based on information published by ^{NOAA (2017)}.

Generalized Additive Models of Location, Scale and Shape parameters-GAMLSS

GAMLSS models assume that the response variable Y (Q_max) has a cumulative probability distribution function, whose parameters may be a function of one or more explanatory variables (climatic and reservoir indices). That is, for yt independent variables at time t=1,2,…,n; there is a function FY(yt│θt), where the parameters θt=µt, σt,υt,τt can change based on a set of m explanatory variables Xmt=x1t,x2t,…,xmt through a monotonic link function gk θk presented in Equation (3) (^{Stasinopoulos, Rigby, Vlasios, Heler, & Bastiani De, 2015}):

Where θk and ηk are vectors of length n; Xk is an array of covariates of order n x m; βk is a vector of parameters of length m; hjk(xjk) represent smoothing in the distribution parameters; and xjk is a vector of covariates for j = 1,2,…, m.

As a smoothing function, the B-splines, pb () which are polynomial parts used in additive models, are evaluated, with the advantage that they minimize the degrees of freedom (df) of the model, through different methods e.g.: Maximum Likelihood (ML), GAIC, etc. More information can be found in ^{Stasinopoulos et al. (2015)}.

Non-stationary return period

The return period T is an indicator of the rarity of flooding, that is, the average time that elapses until, for the first time, a flood exceeds a given value (yp0). It is commonly used to define the design flow rate of a hydraulic work or to reference the level of flood threat. Considering that the statistics of this event follow a geometric distribution, the probability that the first flood exceeds yq0 at time x is fx=px∏t=1x-11-pt,  x=1,2,…,xmax (^{Salas & Obeyskera, 2014}), where f1=p1 and xmax is the time at which pt becomes unity. If the probability pt is constant over time fx=PX=x=1-px-1p, therefore, the return period or Expected Wait Time (EWT) is:

Where p and q are the probability of exceedance and without exceedance, respectively. When the probability pt is non-stationary, then the return period changes over time:

The values of the probability of exceedance 𝑝 𝑡 are obtained from pt=1-FYyq0,θt using the available data. yp0 is a reference flow rate for a non-stationary Frequency Analysis model of which the parameter θt varies in accordance with the covariates Xmt (^{Obeyskera & Salas, 2016}). Considering Equation (5), in the hydraulic design and management of water resources, questions arise: What is the value of T of a historical event, associated with a non-stationary flood model? What is the design flow rate, if the value of T is variable in time? Because of this, ^{Salas and Obeyskera (2014)} propose to determine the non-stationary T value, using a non-homogeneous geometric distribution:

Where the values pt=1-qt, are obtained from pt=1-FYyq0,θt, for the preset yq0 and the non-stationary statistical model. ^{Salas Obeysekera and Vogel (2018)}, and ^{Serinaldi (2015)} present another simplified way of estimating the non-stationary return period, based on the concept of “Average Annual Risk” economic analysis (AAR) expressed as: AARn=p-=1/np1+p2+…+pn. Where it is assumed that for a period n, the risk can be described by the mean of the sequence of probabilities of exceedance (p1,p2,…,pn), due to this, non-stationary T equals:

^{Serinaldi (2015)} argues that one of the most important indicators in the planning and design of hydraulic works is the Risk of Failure 𝑅 or probability 𝑅 that an extreme event will be observed at least once yq0 during the useful life n of a works. This expression being the most appropriate or direct to define the design flow rate against an acceptable level of risk since it depends on the probability of exceedance and not on the return period. When i.i.d. conditions exist, i.i.d. R=PX≤n=FXn=∑x=1nfx=1-(1-p)n (^{Chow, Maidment, & Mays, 1994}); however, under conditions of non-stationarity ^{Salas et al. (2018)} use the concepts Expected Wait Time and Average Annual Risk to determine the risk of non-stationary failure:

Where pt and p- are the probability of exceedance and the probability of multi-year mean exceedance, respectively. Highlighting that REWT and RAAR are solved numerically.

Identification of changes and trends in the flood regime

The non-parametric statistical tests of ^{Mann-Kendall (Mann, 1945}; ^{Yue & Wang, 2002}) and ^{Pettitt (1979)} are applied to the time series to analyze the stationarity. Additionally, in order to address the effects of reservoir operation, the maximum annual flood series were divided into the periods: 1965-1984 and 1986-2015, thereby representing the unregulated regime and the altered regime, respectively.

Adjustment and selection of non-stationary statistical models

GAMLSS models are used to force gradual changes in the location parameter μ based on the ONI, MEI and IE2. covariates. The probability distribution functions selected are: Gumbel (GU), Lognormal (LN) and Gamma (GA). All have two parameters μ, σ, they are suitable to counteract the effects of positive asymmetry, common in series of hydrological extremes. The results will focus on the best models obtained and not on the set of models evaluated.

Two types of statistical models are analyzed: Stationary models (M0), where both distribution parameters are independent; and the models of external covariates (M1), in which the location parameter μ can depend linearly on one or several climatic indices, on the reservoir index or on the combination of both types of explanatory variables. The selection of the models is based on the Generalized Akaike Information Criteria (GAIC) and the Schwartz Bayesian Criteria (SBC) (Equation (10) and Equation (11)):

Where l^ = ln (ML), the maximum likelihood of the model is ML, then k: required penalty associated with the number of parameters of the distribution and df: degrees of freedom. In this work, a penalty k = 3.0 is adopted, in such a way that the degree of complexity of the model does not degrade the ability to describe the behavior of the series and conserves as much as possible the principle of parsimony.

In the absence of a statistic to determine the goodness of fit of the GAMLSS models, ^{Stasinopoulos et al. (2015)} propose to verify the normality of the residuals, for this the Filliben correlation coefficient and the behavior of the graphs without the residuals trend (worm-plot) are evaluated. This ensures that the selected model represents the systematic part and that the remaining (residual) information is white noise.

The return period in non-stationary models

To compare the results obtained in the M0 and M1 models, the EWT and AAR formulations are used. The calculation starts by defining a flow rate yp0 with the stationary model, then the non-stationary model M1 is used to estimate the probability of exceedance, equations 4 to 9 are applied to determine the hydrological risk in terms of non-stationary T and R, considering between other aspects: the moment at which pt becomes unity, the length 𝑛 of records and the predominant pattern of change in M1. It is important to mention that, for each station, R is estimated for the design flow rate of T = 100 years and different values n of working life.

Altered flood regime

Next, the changes in the flood regime in La Balsa and Juanchito are analyzed.

Table 1 confirms that the flow rate at La Balsa decreases during the registration period and this trend is significant at α=0.05. In addition, there are points of change in the mean in 1975 and 1984 that coincide with phase changes of the PDO and the commissioning of the reservoir. Due to the proximity of the station to the reservoir output it is possible to associate the non-stationarity with anthropic alteration. Figure 3b shows the decrease in the magnitude of Juanchito's floods between 1965 and 1984, but this change has no statistical significance (Table 1). On the contrary, between 1986 and 2015, a gradual and significant increase in annual maximum flow rate is observed, although this increase does not necessarily respond to a linear correlation (Figure 3b). There are also identified points of change in the mean in 1984 and the 1970s and 1990s that coincide in the first moment with the commissioning of the reservoir, and the rest with PDO phase changes, a high rate of occurrence of ENSO events (five reported in the 1970s and 1980s) (^{Wolter & Timlin, 1998}), and historical floods of the Cauca river in 1971, 1974, 1975, 1982, 1984, 1988, 1997 and 1999, 2008, 2010 and 2011 (^{Enciso et al., 2016}), that affect the mean of the time series.

Figure 3 Variation of annual floods. From left to right La Balsa and Juanchito; a), b) time evolution; c), d) Connection with the IE2 reservoir index; e)-h) Relationship with the ONI and MEI indices, highlighted in blue and red the flow rates observed in La Niña and El Niño years, respectively. 

The following physical processes that may be related to the trends of increase in the flood regime of Juanchito station between 1986 and 2015 are:

The results obtained in Table 1 highlight the difficulty in accepting / rejecting the i.i.d. hypothesis for the 1965-2015 period at Juanchito station, since the non-parametric tests applied indicate stationarity. However, changes in the slope of the trend lines in the RN and RA periods mask the continuous increase in annual flooding in the last 30 years of the analysis period and therefore, the RA sample requires the non-stationary Flooding Frequency Analysis.

Figure 3 shows the behavior of maximum annual floods as a function of time and reservoir and climate indices, which suggest that although there is a moderate linear correlation (0.30<r<0.70), the variability explained by the selected indices is low. This is understandable, given the complex connection between flow rates and climatic variables that do not necessarily conform to linear models, therefore, it may be necessary to include smoothing functions in non-stationary statistical modeling.

Non-stationary statistical models of the flood regime in the Cauca

River

The results of the statistical modeling of the annual maximum daily flow rates for La Balsa (1965-2015) and Juanchito (1986-2015) are presented below.

Table 2 shows that the stationary model (M0) in La Balsa follows the Lognormal distribution (LN2). The Log transformation helps reducing the positive asymmetry of the observations. Then, the floods in Juanchito conform to the Gamma distribution, which is a function that has a smooth form and does not require the Log function to counteract the asymmetry; This last result is similar to that obtained by ^{Enciso et al. (2016)}, who finds that Juanchito's annual maximum daily flow rates in the 1986-2015 period are adequately adjusted to a Log Pearson III distribution, which is a generalization of the Gamma distribution. Next, it can be seen that the GAIC statistic suggests that the M1 models have lower loss of information. Regarding the quality of fit, in all cases, the Filliben correlation coefficients are higher than the critical values, therefore, the hypothesis of normality in the residuals is accepted and the models are adequately adjusted to the observations.

^{Debele, Bogdanowicz and Strupczewski (2017)} suggest that the selection of the distribution function (fdp) is one of the most important decisions for the proper analysis of the models. We consider that the hypothesis tests establish an adequate fitting, but do not necessarily lead to the selection of the best fdp. For this, GAMLSS models provide a broad range of comparison options between different families of fdp and also the consideration of non-stationarity. However, it is necessary to use the expert’s criterion to select a pfd that controls the effect of the positive asymmetry characteristic of the hydrological series and a model without over-parameterization.

Among non-stationary models with the best fit we have:

It is important to highlight that based on the evidence of non-stationarity, it is necessary to adopt methodologies that incorporate change patterns and allow a comparative analysis of results. In this work, the non-stationary models based on co-variables (M1) show a better representation of the variability of the time series, considering that the majority of the observations are within the band of quantiles 1% to 99% of the models (Figure 4a and Figure 4b). ^{Villarini, Serinaldi, Smith, & Krajewski (2009b)} obtained similar results using GAMLSS for the analysis of floods in a US basin, mentioning that the models manage to capture the wide dispersion and nonlinearity of data for percentiles between 5 and 95%.

Figure 4 Statistical models. From left to right: annual floods at La Balsa and Juanchito; a) and b) the temporal variation of different percentiles in the range 1% - 99% of the M1 model (solid lines), the blue dots correspond to the observations. Theoretically, between 60% and 90% of observations are expected to be within the area covered by the 1% and 99% percentiles of the model; c)-f) contain the worm-like graphics of the residuals of the M0 and M1 models, respectively. 

The non-stationary M1 models show - for all the percentiles plotted - that at certain times the magnitude of the variables obtained is different from that estimated under stationary conditions, e.g. during La Niña years, and that increases in the magnitude of flow rates that can affect flood risk indicators such as the return period and the risk of failure for a particular project life can be identified. Regarding this result, ^{López and Francés (2014)}, when evaluating floods in the Northwest of Mexico, also find a significant influence of the ENSO phenomenon in the inter-annual variability of the flood regime, highlighting increases in magnitude during La Niña. In addition to the above, it is necessary to recognize that a limitation of the results of the co-variable models is the uncertainty associated with ignorance about the future, e.g. the behavior of the explanatory variables beyond the registration period (there are no long-term projections for the ENSO indices that may be incorporated into the predictive power of the models) and that other physical processes (not considered in the present study) may be more significant to describe the variability of floods.

Figures 4c-4f evaluate the normality of the residuals of the different models based on their configuration. In a Q-Q graph with no tendencies, which is usually shaped like a worm, one seeks for the continuous red line (trend) to be similar to a straight line, parallel and close to the horizontal axis. On the contrary, when the residuals show accentuated configurations either in S or U, they indicate high asymmetry and / or kurtosis (^{Buuren & Fredriks, 2001}). In both stations, it is evident that the models meet the normality condition. However, of all of them, the M1 model of La Balsa has fewer deviations from that assumption.

Analysis of changes in the return period

Figure 5 and Table 3 show the variations obtained when estimating the period of return T and risk of failure R, for the selected stationary and non-stationary models. Figures 5a and 5b contain the variation of non-stationary T as a function of stationary T, following the methods of Expected Wait Time (EWT) and Average Annual Risk (AAR) (^{Salas et al., 2018}). This paper establishes that the maximum annual flow rates at La Balsa have a long-term tendency to decrease, so that the probabilities of exceedance also suggest a decreasing pattern. Figure 4a establishes the variation of non-stationary 𝑇, highlighting two things: first, for 1≤T≤5 years, the non-stationary model indicates that TEWT<T. This corresponds to an increase in the probability of observing minor floods. However, when using the AAR method, no differences can be seen between non-stationary T and T≤20 years, (TAAR≅T), this situation can be associated with the fact that in the first 20 years of records (1965-1984) the hydrological series was not affected by the reservoir, so the mean of the probability of exceedance remains constant and similar to that of the stationary model. Second, from the aforementioned inflection points, it is appreciated that TEWT>T and TAAR>T, thus, the decrease in the magnitude of the floods during the altered regime (1986-2015) leads to the decrease in hydrological risk expressed as the probability of observing a flood equal to that of the reference or increase in non-stationary T. Said behavior of increase and subsequent decrease of the hydrological risk can be seen directly in Figure 4c where the risk of failure R associated with a T=100 years flood flow rate or the probability of observing, at least once, a flood equal to that of the reference in a period of 𝑛 years of working life. In the case of La Balsa, a REWT,AAR larger than the stationary is observed for works with n<20 years and a lower risk, from 30 years of working life. Another aspect to note is that the non-stationary 𝑅 values in each station are identical regardless of the method of determination, so their use could be a more direct way of establishing the design flow rate of a hydraulic works for a predefined acceptable level of risk.

Figure 5 Variation of the Return Period and Risk of Failure of the statistical models of La Balsa (left) and Juanchito (right). Above, evolution of non-stationary T as a function of stationary T. Below, changes in risk as function of working life. 

In the case of annual floods at Juanchito between 1986 and 2015, an increasing trend has been described (Figure 3b). When comparing the non-stationary return period, Figure 5b shows: 1) When 1<T<7 years TEWT>T and TAAR>T are obtained. Therefore, the co-variable model indicates that the most frequent floods have a lower probability of occurrence in comparison with the stationary model. After seven years, the opposite occurs: TEWT<T and TAAR<T. Although graphically this difference is modest, the decrease of TEWT against T ranges between 14% and 50%, whereas between TAAR and T varies between 15% and 39%, with a slight increase when T is 1,000 years. Among the possible explanations for the evolution of TAAR are the gentle slope of change of the annual maximum flow rates, its significance (α=0.10) and the short length of the time series (30 years). In any case, the information obtained highlights a slight increase in the probability of occurrence of extreme maximum floods. Which leads us to consider that although the management of the reservoir has positive effects on the most common floods, there are external factors (climatic variability, increase in precipitation, changes in the unregulated basin, etc.) that negatively affect the objective of regulation of the rarest and most extreme floods. Similar information is presented in Figure 5d, in which case REWT, AAR<T for a working life of n≤ 25 years, but from said n, the risk failure of the non-stationary model exceeds that of the stationary R.

All the changes in the hydrological risk indicators (Qmax, T y R), are of interest for flood management at Juanchito, not only because this is the target station for volume control in the Salvajina reservoir; but because the unplanned urban growth in the city of Cali, as in many other cities in developing countries, has led to intensive building developments in the flood plains of both banks of the Cauca River. About 900 thousand inhabitants of the city are located in a flood risk area on the right bank of the Cauca river. Flood control is managed through the operation of the reservoir and lateral levees built for a T=30 year flood. Despite these structures, the changes in the hydrological risk identified are corroborated by historical flow rates observed in 2008, 2010, 2011 and 2017; and more than ten floods in the Juanchito area recorded in the last decade (^{El País, 2011}; ^{Enciso et al., 2016}), highlighting the need to provide for a natural disaster risk management program that incorporates a systemic view of the problems and includes change patterns associated with climate variability.