Uncertainty-Aware and Non-Negative Hydrological Forecasting Using Gamma-Likelihood Chained Gaussian Processes for Sustainability-Oriented Water Management

1. Introduction

Water constitutes the most fundamental natural resource as it sustains life, ecosystems, and human societies. As global population growth and development continue, the demand for water has increased significantly, drawing international attention to the sustainable management of hydrological resources [1]. The excessive exploitation of these resources has led to multiple challenges that hinder economic development and increase environmental concerns, such as declining river flows. Consequently, advancing understanding of hydrological processes is now a key focus in natural resource research, as it enhances the ability to anticipate extreme events such as droughts and periods of intense rainfall [2]. Moreover, changes in climate patterns have increased the temporal and spatial variability of water resources, making their behavior increasingly complex [3].

In the Colombian context, a country characterized by high hydrological endowment, the real challenge is no longer the availability of water, but its sustainable and efficient management under conditions of pronounced hydroclimatic variability [4]. Recent studies have highlighted that the increasing frequency of extreme events—such as prolonged droughts and intense precipitation—combined with anthropogenic pressures, significantly affects the reliability of municipal water supply systems and ecosystem services [5,6]. Consequently, there is a growing need to transition from traditional deterministic planning approaches toward data-driven and probabilistic frameworks that explicitly incorporate uncertainty and variability in hydrological processes.

Initial approaches to forecasting water resources rely on physically based deterministic models that describe natural laws governing hydrological time series. Despite their theoretical foundation, these models pose major challenges, as obtaining reliable results demands substantial effort and detailed knowledge of numerous system parameters, often leading to restricted practical applicability [7,8]. An alternative strategy employs statistical and machine learning techniques that are entirely data-driven, allowing the capture of complex patterns directly from observations [9]. Within this scope, autoregressive (AR) models have been extensively applied; however, they are frequently unable to represent the stochastic and nonlinear behavior characteristic of water resources [10].

To address these limitations, Artificial Neural Networks (ANNs) have been widely adopted in hydrological studies due to their ability to model nonlinear relationships. Nevertheless, ANNs present several limitations, including overfitting, sensitivity to local minima, and issues related to vanishing and exploding gradients during training [11]. Long Short-Term Memory (LSTM) networks are among the most widely recognized neural architectures for time series modeling, as they are specifically developed to process sequential data. In particular, LSTMs effectively address the vanishing gradient problem inherent in standard recurrent architectures [12,13,14]. Recent research has demonstrated the effectiveness of LSTM-based recurrent neural network for monthly low-flow prediction at the Basantapur station in India [13] and an alternative methodology was presented in [12], where a wavelet transform was first used to filter noise and remove interfering components from hydrological time series. Despite their flexibility, a key limitation of these neural approaches is their inability to provide information about predictive uncertainty [15]. In practical applications, understanding the level of confidence associated with a forecast is essential for informed decision-making and risk management [16].

Gaussian Processes (GPs) offer a robust stochastic modeling alternative by providing full predictive distributions rather than single-point estimates [17]. By linking observations to multivariate normal distributions, GPs effectively capture the random nature and intrinsic non-linearity of hydrological processes [18,19]. The work reported in [18] assessed the effectiveness of several approaches for daily water resource time series forecasting, including Gaussian Process models. The results underscored the capability of GPs to account for data noise and represent nonlinear temporal dynamics. In watershed studies [20] introduced a new method for quantifying uncertainty, employing Gaussian Processes within a Bayesian sensitivity analysis framework. The main goal of this analysis was to detect and localize the factors that most strongly contribute to model uncertainty. Furthermore, hybrid models integrating LSTMs with GPs have been proposed to preserve nonparametric probabilistic characteristics. Within the Colombian context, the study presented in [21] applied a SARIMA model to forecast monthly precipitation. The findings are particularly significant, as they illustrate how such predictions can support a wide range of hydrological analyses and contribute to effective water resource management.

While Gaussian Processes offer strong probabilistic foundations, many real-world hydrological applications require the simultaneous forecasting of multiple correlated time series, a setting commonly referred to as multi-task or multi-output learning. This problem can be addressed by extending GP frameworks to model multiple outputs jointly, thereby explicitly capturing inter-task dependencies. A widely adopted approach relies on collections of latent GPs combined through coregionalization matrices, commonly known as separable kernel sum (SoS) models [22]. Although these models are highly interpretable, their application in hydrology is often hindered by rapid parameter growth as the number of tasks increases, leading to overparameterization and significant challenges in model estimation and validation [23]. Even with mitigation strategies such as variational inference [24] or low-rank approximations, computational complexity typically scales cubically with the number of tasks, limiting practical applicability in highly expressive settings [25]. Moreover, existing multi-output GP models, including Linear Models of Coregionalization (LMCGP) [25], may assign predictive mass to negative values, violating the inherent non-negativity constraint of streamflow data [26].

This paper presents the Chained Correlated Gaussian Process (ChdGP) as a scalable, probabilistic model for the joint forecasting of multiple hydrological time series in Colombia. Unlike traditional GP frameworks, the proposed approach represents likelihood parameters through latent processes, enabling the use of non-Gaussian likelihoods to enforce non-negativity in predictions. Additionally, we employ a hybrid optimization strategy that combines Adam and Natural Gradient (NG) methods to enhance stability and avoid poor local optima in variational inference. The proposed model is validated using data from 23 strategic reservoirs, demonstrating that ChdGP effectively exploits inter-task correlations and heteroscedastic noise the input space, significantly enhancing forecasting accuracy and uncertainty quantification across multiple horizons.

The remainder of this paper is organized as follows. Section 2 presents the mathematical framework, including the problem formulation, the multi-output Gaussian Process models, and the variational inference approach. Section 3 describes the dataset and experimental setup; and discusses the results, highlighting the implications of probabilistic modeling for hydrological forecasting and sustainable water resource management. Finally, Section 4 concludes the paper by summarizing the main findings, proposing directions for future research.

2. Mathematical Framework

2.1. Problem Setting

Let ${\mathbf{v}}_n\in {\mathbb{R}}^D$ denote the vector of hydrological observations of the D outputs at time instant n. The proposed model uses the complete sequence of these vectors from time $n-T+1$ to time n as input to forecast the resource vector at a future time $n+H$. In this formulation, T corresponds to the model order, while H indicates the forecasting horizon. Accordingly, the input vector ${\mathbf{x}}_n$ is defined as:

$${\mathbf{x}}_n=\left[\begin{array}{c}{\mathbf{v}}_n^{\top }\\ {\mathbf{v}}_{n-1}^{\top }\\ ⋮\\ {\mathbf{v}}_{n-T+1}^{\top }\end{array}\right]\in {\mathbb{R}}^{TD}$$

and the target output vector ${\mathbf{y}}_n$ as follows:

$${\mathbf{y}}_n={\mathbf{v}}_{n+H}.$$

This formulation allows the model to effectively exploit past hydrological observations to generate reliable forecasts of future values. Based on this approach, we construct a dataset $\mathcal{D}={\{{\mathbf{x}}_n,{\mathbf{y}}_n\}}_{n=1}^N=\{\mathbf{X},\mathbf{y}\}$ consisting of N input–output pairs, where $\mathbf{X}\subset \mathcal{X}\in {\mathbb{R}}^{TD}$ represents the dimensional input space. The target is defined as a vector-valued function, with ${\mathbf{y}}_n\in {\mathbb{R}}^D$ containing the observations of all outputs (tasks) corresponding to the same input ${\mathbf{x}}_n$. Thus, ${\mathbf{y}}_n={\left\{y_{d,n}\right\}}_{d=1}^D$, and the complete set of outputs is expressed as $\mathbf{y}={[{\mathbf{y}}_1^{\top },{\mathbf{y}}_2^{\top },\dots ,{\mathbf{y}}_D^{\top }]}^{\top }\in {\mathbb{R}}^{DN}$, where each task-specific output vector d is given by ${\mathbf{y}}_d={[y_{d,1},y_{d,2},\dots ,y_{d,n}]}^{\top }\in {\mathbb{R}}^N$.

2.2. LMCGP for Chained GPs Regression Framework

When the output variable is subject to constraints, such as being restricted to non-negative values, standard probabilistic models like Gaussian Processes may encounter limitations. Deterministic neural networks can naturally enforce such restrictions by applying activation functions like the Rectified Linear Unit at the output layer to ensure non-negative predictions. However, Gaussian Processes require appropriate transformations of random variables in order to properly define their probability distributions. The authors in [26] proposed a framework for handling output constraints by selecting an appropriate distribution likelihood family, with each parameter modeled by a single-output Gaussian Process. Additionally, to extend this formulation to multiple outputs, the method introduced in [27] provides a way to incorporate inter-task correlations while simultaneously managing constraints for each task.

Accordingly, we consider the outputs ${\mathbf{y}}_d$ to be conditionally independent given a set of parameters ${\mathbf{\theta }}_d={[{\theta }_{d,1},{\theta }_{d,2},\dots {\theta }_{d,J_d}]}^{\top }\in {\mathbb{R}}^{J_d}$, where $J_d$ represents the number of parameters that characterize the likelihood distribution of the d-th output. In the same manner, the complete parameter vector is defined as $\mathbf{\theta }={[{\mathbf{\theta }}_1^{\top },{\mathbf{\theta }}_2^{\top },\dots ,{\mathbf{\theta }}_D^{\top }]}^{\top }\in {\mathbb{R}}^J$, with $J={\sum }_{d=1}^D{J}_d$.

Considering a Gaussian likelihood, we can factorize it over the training set instances. The likelihood of observing all output realizations $\mathbf{y}$ given $\mathbf{\theta }$ can be expressed as a product of individual likelihoods:

$$p(\mathbf{y}\mid \mathbf{\theta })=\prod _{d=1}^{D}p({\mathbf{y}}_d\mid {\mathbf{\theta }}_d)$$

(1)

Each j-th parameter of the d-th likelihood distribution ${\theta }_{d,j}$ is a deterministic transformation of a latent Gaussian Process prior realization $f_{d,j}$, given by ${\theta }_{d,j}=g_{d,j}\left(f_{d,j}\right)$. For further detail, we define:

• ${\mathbf{f}}_{d,j}={[f_{d,j,1},f_{d,j,2},\dots ,f_{d,j,N}]}^{\top }\in {\mathbb{R}}^N$ as the vector of latent function evaluations for the j-th parameter of likelihood associated to the d-th output, where $f_{d,j,n}=f_{d,j}\left({\mathbf{x}}_n\right)$.
• ${\mathbf{f}}_d={[{\mathbf{f}}_{d,1}^{\top },{\mathbf{f}}_{d,2}^{\top },\dots ,{\mathbf{f}}_{d,J_d}^{\top }]}^{\top }\in {\mathbb{R}}^{J_{d}N}$ as the stacked vector of all latent function evaluations of likelihood associated to the d-th output
• $\mathbf{f}={[{\mathbf{f}}_1^{\top },{\mathbf{f}}_2^{\top },\dots ,{\mathbf{f}}_D^{\top }]}^{\top }\in {\mathbb{R}}^{JN}$ as the grand vector of all evaluations across all outputs.

The conditionally independent likelihood is then formulated as follows:

$$p(\mathbf{y}\mid \mathbf{\theta })=p(\mathbf{y}\mid \mathbf{f})=\prod _{d=1}^{D}p({\mathbf{y}}_d\mid {\mathbf{f}}_d)=\prod _{d=1}^D\prod _{n=1}^{N}p(y_{d,n}\mid f_{d,1,n},f_{d,2,n},\dots f_{d,J_d,n})$$

(2)

This formulation introduces J latent parameter functions $f_{d,j}\left(\mathbf{x}\right)$, each governed by a GP prior, therefore:

$${\mathbf{f}}_d\left(\mathbf{x}\right)\sim \mathcal{GP}\left(\mathbf{0},{\mathbf{K}}_d(\mathbf{x},{\mathbf{x}}^{\prime })\right)$$

(3)

with $\mathbf{f}:\mathbf{x}\to {\mathbb{R}}^{J_d}$ representing the vector-valued function responsible for mapping the input space. If all latent functions are assumed independent, the joint prior over $\mathbf{f}$ evaluated at the training inputs $\mathbf{X}$ follows:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\mathbf{f}}_1\\ {\mathbf{f}}_2\\ ⋮\\ {\mathbf{f}}_D\end{array}\right]\sim \mathcal{N}\left(\begin{array}{c}\mathbf{0}\end{array},\left[\begin{array}{cccc}{\mathbf{K}}_1& \mathbf{0}& \dots & \mathbf{0}\\ \mathbf{0}& {\mathbf{K}}_2& \dots & \mathbf{0}\\ ⋮& ⋮& \ddots & ⋮\\ \mathbf{0}& \mathbf{0}& \dots & {\mathbf{K}}_D\end{array}\right]\right)\end{array}$$

(4)

Here, ${\mathbf{K}}_d\in {\mathbb{R}}^{J_{d}N\times J_{d}N}$ denotes the covariance matrix associated with the parameters that characterize the likelihood distribution of the d-th output. Under the assumption of independence, the grand covariance matrix of the multi-output model, $\mathbf{K}\in {\mathbb{R}}^{JN\times JN}$ , present a diagonal by block structure. This formulation is referred to as the Independent Gaussian Process (IGP) model and is equivalent to employing J separate single-output Gaussian Processes. Although this approach can be effective for multi-output prediction, the main motivation for jointly learning multiple outputs is to exploit the dependencies among them, allowing each output to inform and enhance the predictions of the others. Such dependencies can be incorporated by introducing cross-covariance terms to replace the zero entries in the covariance matrix, as suggested in earlier works [25,28], where proposed the Linear Model of Coregionalization Gaussian Process (LMCGP).

Our goal is to represent $\mathbf{f}\left(\mathbf{x}\right)$ with a multi-output Gaussian process that explicitly captures correlations among the latent parameter functions $f_{d,j}\left(\mathbf{x}\right)$. To this end, we adopt the Linear Model of Coregionalization Gaussian Process (LMCGP), which encodes the dependence between any pair of latent variables $f_{d,j}\left(\mathbf{x}\right)$ and $f_{d^{\prime },j^{\prime }}\left({\mathbf{x}}^{\prime }\right)$ through an instantaneous, time-invariant linear combination of a set of independent Gaussian processes ${\left\{u_q\left(\mathbf{x}\right)\right\}}_{q=1}^Q$. Each of these processes is defined by its own covariance function $k_q(\mathbf{x},{\mathbf{x}}^{\prime })$ and shares the same block diagonal structure specified in Equation (4).

Under this formulation, every latent process $f_{d,j}$ can be expressed as follows:

$$f_{d,j}\left(\mathbf{x}\right)=\sum _{q=1}^Q{a}_{d,j,q}{u}_q\left(\mathbf{x}\right)u_q\left(\mathbf{x}\right)\sim \mathcal{GP}\left(0,k_q(\mathbf{x},{\mathbf{x}}^{\prime })\right)$$

(5)

being $k_q$ the kernel function of the independent process q, government by the parameters set ${\Phi }_q$. Due the independency assumption $\mathrm{cov}\{u_q\left(\mathbf{x}\right),u_{q^{\prime }}\left({\mathbf{x}}^{\prime }\right)\}=k_q(\mathbf{x},{\mathbf{x}}^{\prime }){\delta }_{q{q}^{\prime }}$. The cross-covariance function of the latent parameter GP $\mathbf{f}$ is as follows:

$$\begin{array}{cc}\hfill k_{(d,j),(d^{\prime },j^{\prime })}(\mathbf{x},{\mathbf{x}}^{\prime })& =\mathrm{cov}\{f_{d,j}\left(\mathbf{x}\right),f_{d^{\prime },j^{\prime }}\left({\mathbf{x}}^{\prime }\right)\}\hfill \\ \hfill & =\sum _{q=1}^Q\sum _{q^{\prime }=1}^Q{a}_{d,j,q}{a}_{d^{\prime },j^{\prime },q^{\prime }}\mathrm{cov}\{u_q\left(\mathbf{x}\right),u_{q^{\prime }}\left({\mathbf{x}}^{\prime }\right)\}\hfill \\ \hfill & =\sum _{q=1}^Q{a}_{d,j,q}{a}_{d^{\prime },j^{\prime },q}{k}_q(\mathbf{x},{\mathbf{x}}^{\prime })\hfill \\ \hfill & =\sum _{q=1}^Q{b}_{(d,j),(d^{\prime },j^{\prime })}^q{k}_q(\mathbf{x},{\mathbf{x}}^{\prime })\hfill \end{array}$$

(6)

where $b_{(d,j),(d^{\prime },j^{\prime })}^q=a_{d,j,q}{a}_{d^{\prime },j^{\prime },q}$ encodes the interactions among the outputs models by q-th indepent process, meanwhile, $k_q(\mathbf{x},{\mathbf{x}}^{\prime })$ mange the interaction among input space view from q-th independent process. This kernel formulation, termed as the sum of separable kernels (SoS kernels), enables the cross-covariance function to operate in an extended input space $\{1,2,\dots ,D\}\times \mathbf{X}$. Instead of a scalar kernel function, we now have a matrix kernel function:

$$\mathbf{k}(\mathbf{x},{\mathbf{x}}^{\prime })=\sum _{q=1}^Q{\mathbf{B}}_q{k}_q(\mathbf{x},{\mathbf{x}}^{\prime })$$

(7)

with ${\mathbf{B}}_q\in {\mathbb{R}}^{D\times D}$ as the coregionalization matrix, with elements ${\left({\mathbf{B}}_q\right)}_{(d,j),(d^{\prime },j^{\prime })}=b_{(d,j),(d^{\prime },j^{\prime })}^q$.

The generative process can be summarized as follows. First, the latent vector $\mathbf{f}$ is sampled from a multivariate normal distribution, $\mathbf{f}\sim \mathcal{N}(\mathbf{0},\mathbf{K})$, where $\mathbf{K}\in {\mathbb{R}}^{JN\times JN}$ is the covariance matrix with elements $k_{(d,j),(d^{\prime },j^{\prime })}(\mathbf{x},{\mathbf{x}}^{\prime })$, computed in some set of inputs. Once we obtain $\mathbf{f}$, we evaluate the vector of parameters $\mathbf{\theta }={[{\mathbf{\theta }}_1^{\top },{\mathbf{\theta }}_2^{\top },\dots ,{\mathbf{\theta }}_D^{\top }]}^{\top }$. Subsequently, the output vector $\mathbf{y}={[{\mathbf{y}}_1^{\top },{\mathbf{y}}_2^{\top },\dots ,{\mathbf{y}}_D^{\top }]}^{\top }$ is generated, where each component ${\mathbf{y}}_d$ is obtained by sampling the conditional distribution $p({\mathbf{y}}_d\mid {\mathbf{\theta }}_d)$.

2.3. Variational Inference

The computational complexity of Gaussian Process (GP) models is dominated by the inversion of the covariance matrix $\mathbf{K}$. A common strategy to alleviate this limitation is to introduce a reduced set of $M\ll N$ inducing locations, denoted by $\mathbf{Z}$. Instead of directly working with the latent processes ${\mathbf{f}}_{d,j}$, the model is reformulated in terms of inducing variables associated with a set of Q independent latent processes ${\left\{u_q\right\}}_{q=1}^Q$. For each process $u_q$, the corresponding inducing variables are defined as ${\mathbf{u}}_q\in {\mathbb{R}}^M$, containing the evaluations of $u_q$ at the inducing locations $\mathbf{Z}\in {\mathbb{R}}^M$, which are obtained from the input space $\mathbf{X}$. That is, ${\mathbf{u}}_q={[u_q\left({\mathbf{z}}_1\right),u_q\left({\mathbf{z}}_2\right),\dots ,u_q\left({\mathbf{z}}_M\right)]}^{\top }\in {\mathbb{R}}^M$ becomes the vector of inducing variables for the inducing locations set and $\mathbf{u}={[{\mathbf{u}}_1^{\top },{\mathbf{u}}_2^{\top },\dots ,{\mathbf{u}}_D^{\top }]}^{\top }\in {\mathbb{R}}^{QM}$ the grand vector of all inducing variables across all inducing locations.

Since the joint Gaussian prior distribution can be written as $p(\mathbf{f},\mathbf{u})=p(\mathbf{f}\mid \mathbf{u})p\left(\mathbf{u}\right)$, the latent functions ${\mathbf{f}}_{d,j}$ are conditionally independent given $\mathbf{u}$. Hence, the conditional distribution is expressed as:

$$p(\mathbf{f}\mid \mathbf{u})=\prod _{d=1}^D\prod _{j=1}^{J_d}p({\mathbf{f}}_{d,j}\mid \mathbf{u})=\prod _{d=1}^D\prod _{j=1}^{J_d}\mathcal{N}\left({\mathbf{f}}_{d,j}\mid {\mathbf{K}}_{\mathbf{u},(d,j)}^{\top }{\mathbf{K}}_{\mathbf{u},\mathbf{u}}^{-1}\mathbf{u},K_{(d,j)}-{\mathbf{K}}_{\mathbf{u},(d,j)}^{\top }{\mathbf{K}}_{\mathbf{u},\mathbf{u}}^{-1}{\mathbf{K}}_{\mathbf{u},(d,j)}\right)$$

(8)

Let $K_{(d,j)}\in {\mathbb{R}}^{N\times N}$ denote the covariance matrix corresponding to the j-th latent parameter function for the d-th output, ${\mathbf{f}}_{d,j}$. Additionally, ${\mathbf{K}}_{\mathbf{u},\mathbf{u}}\in {\mathbb{R}}^{QM\times QM}$ represents the covariance matrix for the inducing variables at inducing locations $\mathbf{u}$. Finally, ${\mathbf{K}}_{\mathbf{u},(d,j)}\in {\mathbb{R}}^{QM\times N}$ signifies the covariance matrix between the j-th latent parameter of the d-th output, ${\mathbf{f}}_{d,j}$, and the inducing variables $\mathbf{u}$.

The posterior joint distribution over the latent variable $\mathbf{f}$ and inducing variable $\mathbf{u}$ takes the following form:

$$p(\mathbf{f},\mathbf{u}\mid \mathbf{y})=p(\mathbf{f}\mid \mathbf{u})p(\mathbf{u}\mid \mathbf{y})$$

(9)

One of the main issues of the Chained GP is that exact posterior inference is generally intractable due to the presence of non-Gaussian likelihoods. We use variational inference to compute a lower bound $\mathcal{L}$ for the marginal log-likelihood $\log p\left(\mathbf{y}\right)$, and for approximating the joint posterior distribution $p(\mathbf{f},\mathbf{u}\mid \mathbf{S})$ as follows:

$$p(\mathbf{f},\mathbf{u}\mid \mathbf{S})\approx q(\mathbf{f},\mathbf{u})=p(\mathbf{f}\mid \mathbf{u})q\left(\mathbf{u}\right)=\prod _{d=1}^D\prod _{j=1}^{J_d}p({\mathbf{f}}_{d,j}\mid \mathbf{u})\prod _{q=1}^{Q}q\left({\mathbf{u}}_q\right)$$

(10)

where $q\left({\mathbf{u}}_q\right)=\mathcal{N}({\mathbf{u}}_q\mid {\mathbf{\mu }}_q,S_q)$ are variational distributions assumed Gaussian whose parameters ${\{{\mathbf{\mu }}_q,S_q\}}_{q=1}^Q$ must be optimized. The lower bound $\mathcal{L}$ for $\log p\left(\mathbf{y}\right)$ is obtained as follows:

$$\begin{array}{cc}\hfill logp\left(\mathbf{y}\right)& =log\int \int p(\mathbf{y}\mid \mathbf{f})p(\mathbf{f}\mid \mathbf{u})p\left(\mathbf{u}\right)d\mathbf{f}d\mathbf{u}\hfill \\ \hfill & =log\int \int q(\mathbf{f},\mathbf{u}){\displaystyle \frac{p(\mathbf{y}\mid \mathbf{f})p(\mathbf{f}\mid \mathbf{u})p\left(\mathbf{u}\right)}{q(\mathbf{f},\mathbf{u})}}d\mathbf{f}d\mathbf{u}\hfill \\ \hfill & \ge \int \int q(\mathbf{f},\mathbf{u})log{\displaystyle \frac{p(\mathbf{y}\mid \mathbf{f})p(\mathbf{f}\mid \mathbf{u})p\left(\mathbf{u}\right)}{q(\mathbf{f},\mathbf{u})}}d\mathbf{f}d\mathbf{u}=\mathcal{L}\hfill \end{array}$$

(11)

We can further simplify $\mathcal{L}$ to obtain:

$$\begin{array}{cc}\hfill \mathcal{L}& =\int \int q(\mathbf{f},\mathbf{u})log{\displaystyle \frac{p(\mathbf{y}\mid \mathbf{f})p(\mathbf{f}\mid \mathbf{u})p\left(\mathbf{u}\right)}{q(\mathbf{f},\mathbf{u})}}d\mathbf{f}d\mathbf{u}\hfill \\ \hfill & =\int \int q(\mathbf{f},\mathbf{u})logp(\mathbf{y}\mid \mathbf{f})d\mathbf{f}d\mathbf{u}-\int q\left(\mathbf{u}\right)log{\displaystyle \frac{q\left(\mathbf{u}\right)}{p\left(\mathbf{u}\right)}}d\mathbf{u}\hfill \\ \hfill & =\int q\left(\mathbf{f}\right)logp(\mathbf{y}\mid \mathbf{f})d\mathbf{f}-\sum _{q=1}^Q\mathrm{KL}\{q\left({\mathbf{u}}_q\right)\Vert p\left({\mathbf{u}}_q\right)\}\hfill \\ \hfill & ={\mathbb{E}}_{q\left(\mathbf{f}\right)}\{logp(\mathbf{y}\mid \mathbf{f})\}-\sum _{q=1}^Q\mathrm{KL}\{q\left({\mathbf{u}}_q\right)\Vert p\left({\mathbf{u}}_q\right)\}\hfill \\ \hfill & =\sum _{d=1}^D{\mathbb{E}}_{q\left({\mathbf{f}}_{d,1}\right),\dots ,q\left({\mathbf{f}}_{d,J_d}\right)}\{logp({\mathbf{y}}_d\mid {\mathbf{f}}_{d,1},\dots ,{\mathbf{f}}_{d,J_d})\}-\sum _{q=1}^Q\mathrm{KL}\{q\left({\mathbf{u}}_q\right)\Vert p\left({\mathbf{u}}_q\right)\}\hfill \\ \hfill & =\sum _{d=1}^D\sum _{n=1}^N{\mathbb{E}}_{q\left(f_{d,1,n}\right),\dots ,q\left(f_{d,J_d,n}\right)}\{logp(y_{d,n}\mid f_{d,1,n},\dots ,f_{d,J_d,n})\}-\sum _{q=1}^Q\mathrm{KL}\{q\left({\mathbf{u}}_q\right)\Vert p\left({\mathbf{u}}_q\right)\}\hfill \end{array}$$

(12)

where the approximate marginal posterior for ${\mathbf{f}}_{d,j}$ is $q\left({\mathbf{f}}_{d,j}\right)=\int p({\mathbf{f}}_{d,j}\mid \mathbf{u})q\left(\mathbf{u}\right)d\mathbf{u}$ which have a close form giving by

$$q\left({\mathbf{f}}_{d,j}\right)=\mathcal{N}\left({\mathbf{f}}_{d,j}\mid {\mathbf{K}}_{\mathbf{u},(d,j)}^{\top }{\mathbf{K}}_{\mathbf{u},\mathbf{u}}^{-1}\mu ,K_{(d,j)}-{\mathbf{K}}_{\mathbf{u},(d,j)}^{\top }{\mathbf{K}}_{\mathbf{u},\mathbf{u}}^{-1}(\mathbf{S}-{\mathbf{K}}_{\mathbf{u},\mathbf{u}}){\mathbf{K}}_{\mathbf{u},\mathbf{u}}^{-1}{\mathbf{K}}_{\mathbf{u},(d,j)}\right)$$

(13)

with $\mu ={[{\mu }_1^{\top },{\mu }_2^{\top },\dots ,{\mu }_Q^{\top }]}^{\top }$ and $\mathbf{S}$ denoting a block-diagonal matrix with blocks given by $S_q$, and $\mathrm{KL}\{q\left({\mathbf{u}}_q\right)\Vert p\left({\mathbf{u}}_q\right)\}$ representing the Kullback-Leibler divergence between two Gaussian distributions, with a closed form solution. The learnable parameters in our model encompass the variational parameters ${\{{\mathbf{Z}}_q,{\mu }_q,S_q\}}_{q=1}^Q$, and kernel parameters ${\{{\mathbf{B}}_q,{\Phi }_q\}}_{q=1}^Q$. The computational complexity of this model is now $\mathcal{O}\left(JNQ{M}^2\right)$.

Now, consider a test point set ${\mathbf{X}}_{*}\in \mathcal{X}$. Assuming a good approximation of the variational posterior $p(\mathbf{u}\mid \mathbf{y})\approx q\left(\mathbf{u}\right)$, the posterior of latent parameter function vector at test points ${\mathbf{f}}_{*}$ is

$$q\left({\mathbf{f}}_{*}\right)=\int p({\mathbf{f}}_{*}\mid \mathbf{u})q\left(\mathbf{u}\right)d\mathbf{u}.$$

(14)

The predictive distribution for a new observation ${\mathbf{y}}_{*}$, given the observed data $\mathbf{y}$, can thus be approximated as:

$$p({\mathbf{y}}_{*}\mid \mathbf{y})\approx \int p({\mathbf{y}}_{*}\mid {\mathbf{f}}_{*})q\left({\mathbf{f}}_{*}\right)d{\mathbf{f}}_{*}.$$

(15)

which integrates over the latent function vector ${\mathbf{f}}_{*}$ to account for its uncertainty in the prediction of ${\mathbf{y}}_{*}$.

The optimization procedure for ChdGP relies on the estimation of several groups of parameters. These include the kernel hyperparameters, such as characteristic lengthscales and the output scale associated with the exponential quadratic kernel, the likelihood parameters like data noise ${\left\{{\sigma }_{Nd}^2\right\}}_{d=1}^D$, and the variational distribution parameters ${\{{\mathbf{Z}}_q,{\mu }_q,S_q\}}_{q=1}^Q$. According to [29], there is a strong dependency between the variational parameters and the other parameters, making the model extremely sensitive to any changes in these variables. Additionally, since the ELBO loss function is generally non-convex, stochastic gradient optimization methods often converge to poor local minima.

Recent studies have shown that optimization with Natural Gradient (NG) can effectively address the challenges described above. In particular, the work presented in [30] highlights the benefits of applying NG for variational GPs optimization, demonstrating that it can be implemented with little additional complexity and how to combine it with the Adam optimizer , providing a hybrid optimization framework. Following this approach, we adopt an optimization scheme in which the variational parameters ${\{{\mathbf{Z}}_q,{\mu }_q,S_q\}}_{q=1}^Q$ are updated using Natural Gradient (NG), while the remaining model parameters are optimized with Adam; we denote this strategy as Adam + NG. This hybrid method enables optimization performance comparable to that obtained when training an exact model using Adam alone.

2.4. Model Setup

The proposed methodology constructs the LMCGP covariance function, as defined in Equation (7),using the widely applied squared exponential kernel, shown in Equation (16). This kernel facilitates smooth data mapping:

$$k_q\left(\mathbf{x},{\mathbf{x}}^{\prime }\mid {\Phi }_q\right)=\exp \left(-{\displaystyle \frac{1}{2}}{(\mathbf{x}-{\mathbf{x}}^{\prime })}^{\top }{\Phi }_q^{-2}(\mathbf{x}-{\mathbf{x}}^{\prime })\right)$$

(16)

where the diagonal matrix ${\Phi }_q=\mathrm{diag}{\left\{{\Delta }_{ql}\right\}}_l^{DT}\in {\mathbb{R}}^{DT\times DT}$ contains the length scale factors ${\Delta }_{ql}\in {\mathbb{R}}^{++}$ from each input dimension. Upon examining Equation (16), one might question the absence of an output scale parameter. However, a detailed look at Equation (7) shows that the elements in the matrix ${\mathbf{B}}_q$ perform the function of rescaling the exponential term in the $k_q$ kernel.

The performance of the models is evaluated using the Negative Log Predictive Density (NLPD), due to its flexibility in accommodating various likelihood functions, which is particularly useful when Gaussian assumptions are not appropriate. The NLPD metric is defined as:

$$\begin{array}{cc}\hfill \mathrm{NLPD}& =-{\displaystyle \frac{1}{N_{*}}}logp({\mathbf{y}}_{*}\mid \mathbf{y})\hfill \\ \hfill & =-{\displaystyle \frac{1}{N_{*}}}\sum _{n=1}^{N_{*}}logp({\mathbf{y}}_{n_{*}}\mid \mathbf{y})\hfill \end{array}$$

(17)

3. Results and Discussions

3.1. Dataset Collection

The hydrological forecasting task used to validate the proposed model is based on time-series data of streamflow contributions from 23 reservoirs in Colombia. This dataset is particularly relevant within the context of sustainable water resource management, as streamflow dynamics reflect the availability, variability, and regulation of water resources that ultimately support multiple uses, including potable water supply, agriculture, and energy production. From a systemic perspective, reservoir inflows constitute a key indicator of water availability within the water–energy nexus, enabling the analysis of how hydroclimatic variability propagates through interconnected resource systems.

The streamflow contributions were recorded on a daily basis from 1 January 2010 to 28 February 2022, yielding a total of 4442 observations. Although these contributions originate from volumetric hydrological processes, they are reported in kilowatt-hours (kWh) by hydroelectric operators, as this unit directly represents the effective energy potential associated with water inflows under operational conditions. This representation does not alter the underlying hydrological information but instead provides a functional transformation that integrates hydraulic head, turbine efficiency, and system constraints. Consequently, modeling streamflow in kWh allows capturing the operational relevance of water availability, facilitating the analysis of how variability in hydrological inputs impacts resource allocation decisions. In the context of this research, this perspective is consistent with a broader sustainability framework, where efficient water use is evaluated not only in volumetric terms but also in terms of its contribution to essential services within interconnected systems. The spatial distribution of the reservoirs is shown in Figure 1, which provides geographical context for the forecasting problem. It is worth noting that the reservoirs are primarily located along mountainous regions to take advantage of streamflow generated by precipitation events.

Figure 2 illustrates the distribution of streamflow contributions for each reservoir. First, a substantial variation in streamflow levels across reservoirs can be observed, which is mainly attributed to differences in their power generation capacities. Second, the majority of reservoirs exhibit non-Gaussian behavior, including platykurtic distributions (e.g., reservoir C), negative skewness (reservoir A), positive skewness (reservoir B), heavy-tailed behavior (reservoir F), and bimodal patterns (reservoir O). These heterogeneous distributions arise from factors such as irregular rainfall events, extended dry periods characteristic of Colombian climate conditions, and operational strategies aimed at either conserving water or maximizing energy production. In addition, reservoir Q, which consistently reports zero streamflow, provides a useful benchmark for evaluating forecasting performance in trivial scenarios. For model validation, the final two years of observations are used as the testing set, while the preceding ten years constitute the training set, resulting in $N=3554$ training samples.

3.2. Model Setup and Hyperparameter Tuning

The model setup takes into account that the time series cannot take negative values, meaning that a zero probability must be assigned to impossible scenarios. Fortunately, the Chained GP framework offers the possibility of using non-Gaussian likelihoods to model the data under the latent process, termed ChdGamma. Therefore, we propose using a Gamma likelihood function, which has positive support, as given by:

$$p(\mathbf{y}\mid \mathbf{f})=\prod _{d=1}^D\prod _{n=1}^N\mathcal{G}amma(y_{d,n}\mid g_{d,1}\left(f_{d,1,n}\right),g_{d,2}\left(f_{d,2,n}\right))$$

(18)

where the deterministic functions $g_{d,1}(·)=g_{d,2}(·)=ln(exp(·)+1)$ represent the softplus function, ensuring positivity in the modeling of the shape and scale parameters for the Gamma likelihood. This approach allows the model to respect the inherent non-negativity constraints of the data more effectively, thereby providing a more accurate and realistic representation of the underlying processes.

As baseline models, the LMC and Chained GP with Gaussian likelihood (ChdNormal) are adopted. These models consider a Gaussian likelihood function given by:

$$p(\mathbf{y}\mid \mathbf{f})=\prod _{d=1}^D\prod _{n=1}^N\mathcal{N}(y_{d,n}\mid g_{d,1}\left(f_{d,1,n}\right),g_{d,2}\left(f_{d,2,n}\right))$$

(19)

where the identity function $g_{d,1}(·)$ is used to represent the mean parameter, while the softplus function for $g_{d,2}(·)=ln(exp(·)+1)$ ensures the positivity of the variance parameter. The choice of the softplus function is based on empirical findings that it better prevents overflow compared to other functions, such as the exponential function. All models hereafter presented are optimized using the Adam + NG strategy. Prior to training, the time series are normalized to the interval $[0,1]$ to handle their non-negative nature and improve numerical stability.

We tuned the hyperparameter Q for the chained GP models over a range from 1 to 46, evaluating the performance using the NLPD metric. It is important to note that other metrics mentioned earlier are no longer suitable due to the non-Gaussian assumption. The results of this search are shown in Figure 3 for both the Gaussian and Gamma likelihood models. Unlike the baseline Gaussian model (in Figure 3a), the Gamma likelihood one (in Figure 3b) turns more stable as Q increases, indicating better data learning without overfitting. Based on such observation, we selected $Q=15$ and $Q=26$ as the optimal number of independent GPs for the Gaussian and Gamma models, respectively.

3.3. Performance Analysis

We compare the performance of the LMC, ChdNormal and ChdGamma models after retrained with the newly rescaled data. Figure 4 compares the Negative Log Predictive Density (NLPD) across various prediction horizons for the models. Firstly, the ChdNormal model consistently outperforms the LMC model across all scenarios, which can be attributed to its enhanced ability to model data noise using a Gaussian Process (GP). A particularly notable difference is observed at horizon H = 7, likely due to a week periodicity that the ChdNormal model can effectively exploit. In turn, the ChdGamma model consistently exhibits a lower NLPD across all horizons, indicating superior robustness in handling short-term and long-term predictions over LMC and ChdNormal. The Gamma likelihood, employed in the ChdGamma model, is particularly effective for non-negative and asymmetric data, typical in hydrological datasets, and avoids assigning predictive mass to impossible negative values, thereby improving accuracy and confidence.

For a qualitative analysis, Figure 5 presents the forecasting results of both chained models on four selected reservoirs using a one-day-ahead horizon ($H=1$). The confidence interval corresponds to the probability mass enclosed between the 2.5th and 97.5th percentiles, and the median of the predictive distribution is now also plotted. On the first two selected reservoirs (first two rows) and compared to the traditional Gaussian likelihood, ChdGamma yields narrower confidence intervals, indicating greater confidence in the model’s predictions. Additionally, in these plots, the mean and median predictions are very close, indicating a highly symmetric distribution concerning the mean due to the absence of extreme peaks. On the last two reservoirs, the time series exhibit abrupt peaks. Similar to the ChdNormal forecast, the Gamma likelihood model attempts to account for these peaks with a wider confidence interval. However, unlike the baseline, the ChdGamma setting does not allocate predictive distribution mass to negative values, which is more appropriate for non-negative time series. Consequently, the mean tends to be higher than the median in the presence of peaks, as the predictive distribution favors positive values while recovering symmetry and a narrower confidence interval during periods without peaks.

4. Concluding Remarks and Future Work

This work introduces the Chained Correlated Gaussian Process (ChdGP) model, which models likelihood parameters via GPs. This approach offers greater flexibility and enables the use of non-Gaussian likelihood functions to account for the natural constraints on the outputs, in this particular case, streamflow levels. The results of this study demonstrate that the proposed ChdGP with a Gamma likelihood constitutes an advancement in the probabilistic modeling of hydrological time series, particularly in contexts characterized by high variability, nonlinearity, and structural constraints, such as non-negativity constraints. Beyond the improvements in predictive accuracy, the main contribution of this work lies in its ability to provide a coherent stochastic representation of water availability dynamics, which is essential for decision-making under uncertainty in water resource systems.

From a methodological perspective, the transition from deterministic, point-based forecasting to fully probabilistic frameworks represents a critical step toward addressing the inherent uncertainty of hydrological processes. The superior performance of the ChdGP model, especially when equipped with the Gamma likelihood, confirms that explicitly incorporating distributional constraints and heteroscedasticity yields more realistic and physically consistent predictions. This is particularly relevant in hydrological systems, where extreme events such as peak inflows and drought periods play a disproportionate role in system reliability and risk management. As shown in the experimental results, the proposed model dynamically adapts its predictive variance, widening confidence intervals during high-variability periods while maintaining sharp estimates under stable conditions, thereby providing meaningful uncertainty quantification. A key insight of this work is the relevance of modeling interdependencies across multiple hydrological series. Thanks to the multi-output structure in the Linear Model of Coregionalization, the proposed framework captures shared and task-specific patterns across reservoirs, improving forecasting performance across all horizons. This finding is particularly important in the Colombian context, where hydrological systems are spatially interconnected and influenced by common climatic drivers. Consequently, joint modeling not only enhances predictive accuracy but also provides a more systemic understanding of water dynamics, which is essential for integrated resource management.

From an application standpoint, the ChdGamma model’s ability to generate full predictive distributions opens new opportunities for risk-aware decision-making. For instance, probabilistic forecasts can be directly integrated into optimization models for reservoir operation, early warning systems, and demand management strategies in urban water supply systems. This is particularly relevant in the context of increasing hydroclimatic variability, where traditional deterministic forecasts fail to capture the range of possible scenarios and may lead to suboptimal or vulnerable decisions. Although provided in energy units (kWh), the dataset fundamentally represents transformations of water inflows within the water–energy nexus. Therefore, the proposed modeling approach transcends its immediate application in hydropower forecasting and can be interpreted as a proxy for water availability dynamics. In this sense, the results are directly aligned with the broader objective of promoting sustainable water resource management, as established in the underlying research project. Reliable forecasting of hydrological inputs enables more efficient allocation of water across competing uses, including potable water supply systems, where variability in source availability directly affects service continuity, infrastructure operation, and long-term planning. As highlighted in the manuscript, reservoir inflows encapsulate the combined effects of climatic variability, watershed processes, and operational decisions, making them a valuable indicator of systemic pressures on water resources.

Future research should focus on extending the proposed framework toward integrated water management applications. In particular, coupling probabilistic hydrological forecasts with exogenous variables such as precipitation, water demand, and large-scale climate indices could further enhance predictive performance and interpretability, enabling the identification of consumption patterns and the design of adaptive allocation strategies under uncertainty. Additionally, integrating socio-environmental variables and alternative distributions or hierarchical likelihood structures for categorical data could improve the representation of extreme events and multimodal behaviors observed in some reservoirs, while allowing identification and classification of sustainable consumption profiles. Finally, deploying the proposed methodology within decision-support systems could facilitate its practical adoption by water utilities and policymakers, contributing to the development of resilient and sustainable water management strategies in Colombia and similar contexts.