Figure 3 shows box-plots of storage volume (V), monthly rainfall (R), total phosphorus (TP), total nitrogen (TN), chlorophyll-a (C), cyanobacteria (B) and Secchi transparency (T) measured from 2013 to 2021 for each reservoir of the cascade system: R1 - Castanhão, R2 - Curral Velho, R3 - Pacoti, R4 - Riachão, and R5 - Gavião. The storage volumes and their standard deviations in percent of total capacities (see Table 1) for R1, R2, R3, R4, and R5 were 9±14%, 69±12%, 36±19%, 39±15%, and 86±5%, respectively. This indicates a relatively large variability of volume for the first four reservoirs (R1, R2, R3 and R4), which is a remarkable characteristic of the reservoirs in the region (Rabelo et al., 2021, 2022), while the last reservoir in the cascade (R5) was strategically kept with a high and approximately fixed volume to guarantee water supply for the capital city of Fortaleza and the Pecém industry and port complex. It is also important to highlight that the low storage volume of R1 (9%), the largest reservoir in Ceará, is a response of low inflows, high evaporation rates and high withdrawals to the reservoir cascade. The monthly rainfall presented an increase from about 20 to 70 mm along the cascade, since the system starts in a semiarid area (R1-R4) and then enters a tropical one (R5), as it approaches the coastal zone, where the annual rainfall overcomes the semiarid limit of 800 mm, similarly to previous studies in the region (Freire et al., 2021; Wiegand et al., 2021; Guimarães and Lima Neto, 2023). On the other hand, only a slight variation of up to about ±30% in the mean values of water quality parameters were observed among the reservoirs, especially for the last three (R3-R5), contrasting with previous studies in non-semiarid reservoirs, in which a significant improvement in water quality was noted along the cascade (Pompêo et al., 2017; Soares and Calijuri, 2022). This was attributed to the low water residence time (< 5 days throughout the cascade) and, consequently, to the low capacity of pollutant removal along the man-made Integration channel (IC), which is lower than that observed in cascade systems of reservoirs interconnected by natural rivers (flow velocities about 10-fold lower). Moreover, the internal phosphorus loads both from fish-farming and release from bed sediments may algo play a relevant role in preventing an improvement in water quality from up-to-downstream reservoirs (Moura et al., 2020; Lima Neto et al., 2022). Note that the internal loads in the studied reservoirs are higher than those reported for non-semiarid reservoirs due to warmer water temperature, higher trophic state, more enriched sediments and longer anoxic periods (Rocha and Lima Neto, 2022).
Although the variability of the mean values of the water quality parameters along the cascade was not significant (see Fig. 3), a large temporal variation was observed. Figure 4 shows time series of total phosphorus (TP) concentration for each reservoir of the cascade. A significant increase in TP concentration was observed from 2013 to 2021 in R1-R4 reservoirs (p < 0.05). This increase was caused by diffuse and point-source contributions from soil, agriculture, sewer, livestock and fish-farming to R1, as pointed out by Rocha and Lima Neto (2021, 2022), and this effect propagated to the downstream reservoirs through the IC. However, as reported byLima Neto et al. (2022) and Carneiro et al. (2023), the high seasonal and interannual variation of reservoir water level in tropical semiarid regions strongly impacts TP dynamics, as a result of abrupt changes in the stratification, oxygenation and internal load patterns. This suggests that three processes affect TP concentration in the individual reservoirs: external load from the IC, external load from the local catchment, and internal load from fish-farming and bed sediments. Contrastingly, Gavião reservoir (R5) showed statistically equal values for TP concentration over time (p > 0.05). This occurred because R5 is the oldest reservoir studied herein (see Table 1) and had already reached a hypereutrophic level (TP > 0.05 mg.L− 1) at the beginning of the time series shown in Fig. 4.
Figure 4 also indicates that the temporal variation of TP concentration among the reservoirs presented a different behavior, with the peak values occurring at different times. This implies that the presence of a single variable (V, R, TP, TN, C, B or T) in the empirical equation could not explain the TP concentration in the downstream reservoir satisfactorily (R2 < 0.40, see Moriasi et al., 2015). This non-linear behavior requires indeed the use of multiple regression models to predict TP along the cascade. Table 2 shows the selected predictive equations (highest R²) for the determination of TP concentration from up-to-downstream reservoirs according to the combinations C-1, C-5, C-6 and C-7. A clear trend of R2 decline with the distance between the upstream and downstream reservoirs (RD) was observed, as the R2-values for the combinations C-1, C-5, C-6 and C-7 were 0.66, 0.32, 0.22 and 0.12, respectively (see Fig. 5). Similar R2 trends were also obtained for the other combinations: C-2 (R2 = 0.34), C-8 (R2 = 0.19) and C-9 (R2 = 0.06); C-3 (R2 = 0.49) and C-10 (R2 = 0.37); and C-4 (R2 = 0.23). On the other hand, the combination C-1 resulted in a higher R² than C-4, even though its RD is larger (48 vs. 14 km). This confirms two effects: (1) R1 was the most influential reservoir in the cascade; and (2) the closer the reservoir pairs (lower RD), the higher the predictive ability of the model. These results are consistent with those obtained from Soares and Calijuri (2022) by using a coupled hydrodynamic-biogeochemical model, in which the impact of water quality changes in the first (and largest) reservoir decreased progressively along the subsequent cascade reservoirs.
Table 2
Selected multiple regression equations for the response variable TP for different combinations between reservoir pairs, as defined in Fig. 2.
Combination (reservoir distance - RD) | Equation | R2 |
C-1 (48 km) | TPR2 = 0.231(TPR1) − 3.8E-05(VR1) − 0.008(TN R1) + 0.001(CR1) − 3.4E-08(BR1) + 0.017(TR1) − 1.2E-04(RR1) + 0.106 | 0.66 |
C-5 (172 km) | TPR3 = -0.150(TPR1) − 4.0E-05(VR1) + 9.5E-05(CR1) + 1.8E-08(B R1) + 0.005(T R1) − 1.9E-04(R R1) + 0.131 | 0.32 |
C-6 (178 km) | TPR4 = 0.151(TPR1) − 2.3E-05(VR1) + 0.047(TNR1) − 1.5E-05(CR1) − 6.2E-08(BR1) − 0.009(TR1) − 3.6E-04(RR1) + 0.034 | 0.22 |
C-7 (192 km) | TPR5 = 0.094(TPR1) − 9.2E-06(VR1) + 0.009(TNR1) + 5.5E-04(CR1) − 1.9E-08(BR1) + 0.023(TR1) − 1.4E-04(RR1) + 0.032 | 0.12 |
The adoption of different time delays (TD) between upstream and downstream reservoirs promoted an improvement in the R2-values. Table 3 shows selected multiple regression equations with different TD for the combination C-7. Considering the R2 values, as well as p-values and upper and lower limits for each estimated coefficient of the explanatory variables, the best result (R2 = 0.69) was obtained for TD = 21 months (∼ 630 days) (see Fig. 6). In this case, R2 is within the reference values (0.65 ≤ R2 ≤ 0.80) for the model to be considered “good”, according to Moriasi et al. (2015).
The optimum time delays (highest R2) for the combinations C-1, C-2, C-3, C-4, C-5, C-6, C-7, C-8, C-9 and C-10 were TD = 30, 180, 90, 90, 180, 270, 630, 180, 270 and 180 days, respectively. It is interesting to note that these time delays corresponded respectively to the following ratios of the cumulative water residence time for each of the above combinations: TD/WRT = 1.57, 0.51, 2.01, 1.37, 0.54, 0.72, 1.43, 0.51, 0.64 and 1.63. This results in an average ratio of TD/WRT = 1.09±0.57, suggesting that time delays of the order of the cumulative water residence times of the reservoirs promote an increase in the R2-values. This also suggests that previous data from upstream reservoirs can be used to predict current and future total phosphorus concentration in downstream reservoirs. Observe that the methodology proposed in the present study is simpler than the coupled hydrodynamic-biogeochemical model used by Soares and Calijuri (2022) to predict water quality dynamics in cascade reservoirs interconnected by rivers.
Table 3
Selected multiple regression equations for the response variable TP for the combination C-7 and different time delays, highlighting the best fit (R2 = 0.691) for TD = 21 months.
Time delay (TD) (months) | Equation | R2 |
0 | TPR5 = 0.094(TPR1) − 9.2E-06(VR1) + 0.009(TNR1) + 5.5E-04(CR1) − 1.9E-08(BR1) + 0.023(TR1) − 1.4E-04(RR1) + 0.032 | 0.12 |
3 | TPR5 = 0.018(TNR1) + 9.0E-04(CR1) + 3.05E-08(BR1) + 0.022(TR1) – 3.0E-04(RR1) + 0.012 | 0.25 |
6 | TPR5 = -0.290(TPR1) − 2.57E-05(VR1) + 6.56E-08(BR1) + 0.012(TR1) + 1.1E-04(RR1) + 0.113 | 0.11 |
9 | TPR5 = -0.444(TPR1) − 4.8E-05(VR1) − 5.7E-04(CR1) + 0.021(TR1) + 2.2E-04(RR1) + 0.158 | 0.23 |
12 | TPR5 = − 0.068(TPR1) + 0.014(TNR1) − 7.0E-04(CR1) − 4.89E-08(BR1) − 0.005(TR1) + 0.117 | 0.12 |
15 | TPR5 = − 0.242(TPR1) − 2.86E-05(VR1) − 0.046(TNR1) + 0.022(TR1) + 1.3E-04(RR1) + 0.170 | 0.29 |
18 | TPR5 = − 0.207(TPR1) − 0.039(TN R1) + 1.202E-07(B R1) − 2.20E-04(R R1) + 0.163 | 0.44 |
21 | TPR5 = − 0.268(TPR1) − 4.6E-05(VR1) − 0.024(TNR1) + 2.03E-07(BR1) + 0.052(TR1) + 0.083 | 0.69 |
24 | TPR5 = − 3.0E-05(VR1) + 0.010(TNR1) + 1.2E-07(BR1) + 0.069(TR1) + 9.8E-05(RR1) − 0.032 | 0.67 |
27 | TPR5 = 0.248(TPR1) + 0.016(TNR1) + 5.0E-08(BR1) + 0.055(TR1) − 0.060 | 0.59 |
30 | TPR5 = − 7.0E-06(VR1) − 0.022(TNR1) + 7.9E-04(CR1) + 4.3E-08(BR1) + 0.051(TR1) + 0.007 | 0.37 |
Figure 7 shows a comparison of measured data and values predicted by the equation selected for the combination C-7 (first and last cascade reservoirs) with the optimum time delay of TD = 21 months (see Table 3). Based on Moriasi et al. (2015), a good performance (0.65 ≤ R2 ≤ 0.80) can be seen considering the period of 2013–2021 (R2 = 0.69), while a satisfactory performance (0.40 ≤ R2 ≤ 0.65) was obtained by considering the complete time series (also including data of 2022) (R2 = 0.62). Even though a small drop in R² occurred, the residues found were reasonably small (±36%) and randomly scattered around zero. This gives credence to the adjusted equation, which can be used to predict future TP concentrations by considering TD = 21 months. Moreover, different scenarios of changes in the flow rate of the cascade system can also be assessed by assuming TD = WRT.