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.
Bastidas-Pantoja et al. (2026) studied this question.