Sustainable water allocation, drought mitigation, and operational planning require reliable forecasting models that account for hydroclimatic variability while respecting physical constraints. This study proposes Chd-Gamma, a chained correlated Gaussian Process (GP) framework for multi-output hydrological forecasting. The proposed model extends chained GPs beyond independent or single-output settings by embedding their latent likelihood-parameter functions in a Linear Model of Coregionalization. Chd-Gamma also enhances conventional multi-output GP hydrological forecasting by replacing Gaussian likelihood assumptions with a Gamma likelihood, thereby enforcing non-negativity and representing skewed and heteroscedastic storage distributions. The proposed model was contrasted with the well-known Long Short-Term Memory (LSTM) network, the multi-output Linear Model of Coregionalization GP (LMC), and the chained correlated GP with Gaussian likelihood (Chd-Normal) for forecasting the daily useful storage volumes from 23 Colombian reservoirs recorded from 2010 to 2022 across multiple prediction horizons. The results over a two-year testing period show that Chd-Gamma provides the strongest overall performance across the four metrics considered. Chd-Gamma reduced the mean squared error by 80% with respect to LSTM and 20% relative to Chd-Normal. In terms of probabilistic performance, the average Negative Log Predictive Density (NLPD) improved by up to 21%. Compared to LMC, with narrow prediction intervals but low coverage, and Chd-Normal, also narrow but overcovering, Chd-Gamma achieves near-nominal coverage of 0.992 with a moderate increase in interval width, pointed towards the best calibration–sharpness trade-off. These findings demonstrate that Chd-Gamma improves accuracy and uncertainty representation while maintaining physically consistent forecasts, making it suitable for risk-aware reservoir-operation support.

Suitable forecasting of streamflow plays a fundamental role in the sustainable management of water resources. This is especially important in regions with high hydroclimatic variability, such as Colombia. Reliable predictions of water availability are essential to support decision-making on water allocation, risk mitigation, and the long-term resilience of supply systems, including urban aqueduct services. However, traditional deterministic models often require extensive parameterization. They may lack flexibility. Data-driven approaches such as Long Short-Term Memory networks typically fail to provide meaningful uncertainty quantification. In this work, we introduce the Chained Correlated Gaussian Process. This is a scalable probabilistic framework for jointly forecasting multiple hydrological time series. The proposed model represents likelihood parameters through latent Gaussian Processes. This enables the use of non-Gaussian likelihoods to enforce physical constraints on the outputs. In particular, a Gamma likelihood is integrated into a Linear Model of Coregionalization. This approach captures interdependencies among multiple time series while ensuring non-negativity. To achieve computational scalability, we employ a variational inference framework with shared inducing points. This reduces computational complexity from cubic to linear with respect to the number of observations. Additionally, a hybrid optimization strategy is adopted. It combines Adam and Natural Gradient methods to improve numerical stability and avoid poor local optima. The proposed methodology is validated using twelve years of daily observations from 23 reservoirs in Colombia. These reservoirs serve as proxies for water availability dynamics within interconnected resource systems. Experimental results demonstrate that the ChdGP framework consistently outperforms baseline models across multiple forecasting horizons. The Gamma-likelihood configuration shows superior performance. It provides physically consistent predictions, improves uncertainty quantification, and effectively captures asymmetric behaviors and extreme hydrological events. Performance evaluation based on the Negative Log Predictive Density metric confirms the robustness of the proposed approach under non-Gaussian conditions. The results highlight the potential of the proposed framework as a reliable tool for probabilistic forecasting and risk-aware decision-making in sustainable water resource management. It offers a flexible, scalable methodology for modeling complex hydrological systems under uncertainty.