Comparative Analysis of Statistical, Deep Learning, and Hybrid Approaches for Daily Rainfall Forecasting in the Ecuadorian Andean Region

1. Introduction

Precipitation is a critical component of the hydrological cycle given its pervasive influence on ecosystems and human livelihoods, including agriculture, drinking water supply, and hydroelectric energy generation [1,2,3]. Globally, extreme hydrometeorological events such as droughts and floods are occurring with increasing frequency and intensity under ongoing climate change [4,5,6]. This trend is particularly pronounced in the Ecuadorian Andes, where complex topography and high spatiotemporal rainfall variability drive flash floods, landslides, and localized droughts. In this context, accurate daily rainfall forecasts are essential, particularly for developing countries seeking to mitigate extreme climate impacts through early warning and adaptive water management [7,8,9].

Precipitation is inherently variable in both space and time, rendering its reliable forecasting a persistent challenge. The highly nonlinear, intermittent, and multi-scale nature of precipitation series challenges classical linear approaches, driving the progressive adoption of advanced time-series techniques, deep learning architectures, and hybrid models [1,2,10,11]. Historically, rainfall forecasting has relied on univariate statistical models such as the Autoregressive Integrated Moving Average (ARIMA) and its seasonal variants, valued for their mathematical rigor, interpretability, and low computational requirements. However, these approaches assume linear relationships and relatively simple dependence structures, which often fail to capture the complex rainfall dynamics found in mountain regions characterized by strong topographic heterogeneity [4,12,13,14].

The emergence of deep learning architectures, particularly Long Short-Term Memory (LSTM) networks [15] and Gated Recurrent Units (GRU) [16], has substantially transformed time-series forecasting by enabling the capture of long-range temporal dependencies and nonlinear patterns in meteorological data [10,17,18,19]. Applications across diverse climatic zones and spatial scales—from the northwestern Himalayas [20] to urban environments in the United Kingdom [21] to national-scale studies [22]—have repeatedly demonstrated that LSTM and GRU models frequently outperform traditional statistical methods for precipitation forecasting, while maintaining acceptable computational efficiency for operational use [17,23,24].

In parallel, hybrid models that decompose precipitation time series into multiple frequency components using signal-processing techniques—notably the Discrete Wavelet Transform (DWT) and Variational Mode Decomposition (VMD)—followed by component-wise prediction through neural networks have emerged as promising solutions. These frameworks aim to overcome the limitations of standalone LSTM and GRU models when handling highly volatile, intermittent precipitation signals characterized by extreme events [25,26,27,28]. Signal decomposition can be particularly advantageous in hydrological systems with strong non-stationarity, by isolating trend, seasonal, and stochastic sub-series to enable more precise modeling of each component [28,29].

Despite these advances, a significant literature gap persists for the Ecuadorian Andean context. Previous studies have applied LSTM-based Bayesian frameworks for hourly rainfall nowcasting in a tropical Andean city [30] and Random Forest models for 30-day accumulated precipitation forecasting over an Ecuadorian Andean basin [31]. However, no systematic, methodologically homogeneous, and reproducible comparison exists among ARIMA, LSTM, GRU, DWT-LSTM, and VMD-GRU focused specifically on univariate daily precipitation series from a meteorological network in the Ecuadorian Andean region. This gap prevents the establishment of a consistent evidence base for determining which model family offers greater robustness and operational utility in a mountain environment characterized by high pluviometric variability [30,31].

This study therefore systematically evaluates the predictive capacity of statistical, deep learning, and decomposition-based hybrid models—including ARIMA, LSTM, GRU, DWT-LSTM, and VMD-GRU—for one-step-ahead daily rainfall forecasting at 18 meteorological stations in Tungurahua Province, Ecuador (2013–2025). Model performance is assessed using a comprehensive dual framework of continuous accuracy metrics (MAE, RMSE, NSE, and KGE) and categorical event-detection metrics (POD, FAR, and CSI) evaluated at three operationally relevant precipitation thresholds (1 mm/day, P90, and P95). By jointly analyzing mean error behavior and extreme-event prediction skill, this study aims to identify the limitations of average-based metrics under highly intermittent rainfall regimes, to evaluate whether hybrid decomposition approaches add predictive value, and to provide actionable guidance for operational rainfall forecasting systems in hydroclimate-complex Andean environments.

2. Materials and Methods

2.1. Area of Study and Meteorological Network

The area of study is located in Tungurahua Province, within the inter-Andean region of Ecuador, and falls within the hydrological influence zone of the Pastaza River basin, a major tributary of the Amazon basin [32,33]. The province is characterized by irregular topography, strong altitudinal gradients (ranging from approximately 1529 to 3875 m a.s.l. across the network), and pronounced hydroclimatic variability driven by orographic processes, local convective dynamics, and large-scale atmospheric circulation patterns (Figure 1).

The provincial capital, Ambato, is classified as Cfb (temperate oceanic) under the Köppen–Geiger classification [34,35,36], though the broader study domain encompasses distinct precipitation regimes shaped by slope orientation, valley exposure, and exposure to moisture-laden air masses of Amazonian origin.

Eighteen meteorological stations were selected from the Tungurahua Provincial Government’s monitoring network to represent the spatial variability of daily precipitation across the study area. Although the network comprises 20 stations, two were excluded due to insufficient record length to support a robust comparative evaluation. The spatial distribution of stations is shown in Figure 1, and their main characteristics are summarized in Table 1.

2.2. Data Description

The dataset consists of daily precipitation records from the 18 selected meteorological stations, provided by the Tungurahua Provincial Government. Daily temporal resolution was adopted for its direct relevance to operational hydrometeorological applications, including water resource management, irrigation scheduling, and flood risk analysis [37]. Available records span the period 2013–2025, although the effective observation period varies by station; stations AS and QP have records available only through May 2023 and August 2024, respectively, and were therefore analyzed over their own reliable observation windows. Each time series exhibits the strongly intermittent behavior typical of Andean precipitation regimes, dominated by zero or near-zero values and punctuated by occasional high-intensity events, as illustrated in Figure 2.

2.3. Data Preprocessing

2.3.1. Quality Control and Data Cleaning

Prior to model training, daily precipitation series were subjected to a structured preprocessing workflow to ensure data consistency and cross-station comparability. Quality control examined records for values exceeding the maximum measurement range of the tipping-bucket rain gauges and for positive values below the instrument resolution threshold of 0.1 mm [38]. No observations outside these limits were detected across the 18 stations; consequently, no data removal or correction was required at this stage.

2.3.2. Handling of Missing Values

Precipitation records from 15 stations were complete. At stations PL, BA, and HB, the proportion of missing data was minimal (≤0.07%). Missing values were filled using daily precipitation estimates from the NASA POWER reanalysis product [39], matching the corresponding dates and station locations. This satellite-based gap-filling strategy is an accepted approach when missing values are sparse and do not compromise the overall dataset structure [40,41].

2.3.3. Temporal Partitioning and Dataset Construction

Each station’s precipitation series was sorted chronologically and partitioned into training (70%), validation (15%), and test (15%) subsets using a fixed temporal ratio without randomization, preserving chronological order and preventing information leakage. The prediction task was formulated as a one-step-ahead (t + 1) forecasting problem under a univariate framework, using a sliding-window input of 30 consecutive daily observations to predict the next-day precipitation amount. Supervised input-output datasets for all deep learning and hybrid models were constructed using the same window length, chronological alignment, and test dates to ensure that differences in predictive performance reflect model structure rather than data preparation inconsistencies.

2.3.4. Data Scaling and Shared Preprocessing Artifacts

Min-Max normalization (Equation (1)) was applied exclusively to models incorporating deep learning components, fitted solely on the training subset and subsequently applied—without parameter recalculation—to validation and test subsets. This procedure preserves temporal causality and prevents data leakage [27].

$$x^{*}={\displaystyle \frac{x-x_{min}}{x_{max}-x_{min}}}$$

(1)

where x_max and x_min are the maximum and minimum values of the precipitation series estimated from the training subset, respectively.

ARIMA was fitted directly on precipitation data in original units (mm/day), as scaling is unnecessary for its standard formulation and retaining original units facilitates direct meteorological interpretation of model residuals and confidence intervals [42]. To guarantee strict comparability across all model scripts, a shared set of preprocessing artifacts—including chronological split definitions, event thresholds, fitted scaling objects, and sliding-window datasets—was created and reused for all experiments.

2.4. Experimental Forecasting Framework

2.4.1. Common Experimental Protocol

All forecasting models were developed under a unified experimental protocol based on strict temporal ordering to ensure comparability, reproducibility, and experimental fairness. The forecast horizon was fixed at one day ahead (t + 1), and the event analysis framework was standardized across all models. Rainfall occurrence was evaluated at a threshold of p ≥ 1.0 mm/day, while high- and extreme-rainfall events were assessed using P90 and P95 thresholds computed exclusively from rainy-day observations in the training subset. These thresholds were held constant for validation and test period analyses.

Reproducibility was reinforced by setting a global random seed and enabling deterministic configurations in all main computational libraries where the software and runtime permitted. All deep learning and hybrid models were trained from scratch with no weight transfer between architectures. To maintain a comparable training budget, all models were optimized under the same constraints: a maximum of 100 epochs, early stopping based on validation loss with a patience of 10 epochs, best-weight restoration, and a default batch size of 32.

2.4.2. Hyperparameter Optimization and Model Selection

Hyperparameter optimization was conducted exclusively on training and validation subsets. For deep learning and hybrid models, the search was performed using the Optuna framework [43] over a constrained search space encompassing learning rate, number of layers, hidden units, dropout rate, and batch size, with the objective of minimizing validation RMSE and pruning enabled to terminate unpromising trials efficiently. ARIMA model selection followed a reproducible constrained search over (p, d, q) combinations evaluated first by the Akaike Information Criterion (AIC) on the training subset, with final order selected by minimizing validation RMSE. When two or more configurations produced similar RMSE values, Kling-Gupta Efficiency (KGE) served as a tiebreaker, with preference given to the higher KGE configuration.

Following hyperparameter selection, the best-performing configuration for each model was retrained on the combined training and validation data, and final performance was evaluated once on the unmodified test subset. This train-validate-retrain-test workflow avoids optimistic bias while maintaining a uniform model selection procedure across all evaluated approaches. A fundamental principle of the framework was explicit leakage control: all fitting, parameter estimation, and threshold definition operations were restricted to training data during model development, and validation and test subsets were exposed only to transformations learned from training data.

2.5. Forecasting Models

Five univariate daily rainfall forecasting models from three methodological families were evaluated: ARIMA (statistical), LSTM and GRU (deep learning), and DWT-LSTM and VMD-GRU (decomposition-based hybrid). These models represent a spectrum of temporal complexity and predictive capacity, from conventional linear approaches to nonlinear sequence learning and multi-scale signal decomposition, and were selected to reflect both established benchmarks and current state-of-the-art approaches in hydrological forecasting [10,17,25,26,44,45].

2.5.1. Autoregressive Integrated Moving Average (ARIMA)

ARIMA, introduced by Box and Jenkins [42], is expressed as ARIMA(p,d,q), where p is the order of the autoregressive (AR) component, d is the degree of differencing required to achieve stationarity, and q is the order of the moving-average (MA) component. The mathematical formulation of the model is presented in Equation (2).

$${\mathit{\varphi }}_{\mathit{p}}\left(\mathit{B}\right){\left(1-\mathit{B}\right)}^{\mathit{d}}{\mathit{Z}}_{\mathit{t}}={\mathit{\theta }}_{\mathit{q}}\left(\mathit{B}\right){\mathit{a}}_{\mathit{t}}$$

(2)

Where:

${\mathsf{\varphi }}_p\left(B\right)=1-{\mathsf{\varphi }}_{1}B-{\mathsf{\varphi }}_2{B}^2-\cdots -{\mathsf{\varphi }}_p{B}^p$	Backshift operator (B) AR process
${\mathsf{\theta }}_q\left(B\right)=1-{\mathsf{\theta }}_{1}B-{\mathsf{\theta }}_2{B}^2-\cdots -{\mathsf{\theta }}_q{B}^q$	Backshift operator (B) MA process
$B$	Backshift operator
${\left(1-B\right)}^d$	Differentiating operator
d	Order of differencing

2.5.2. Long Short-Term Memory Neural Network (LSTM)

LSTM networks, introduced by Hochreiter and Schmidhuber [44], address the vanishing-gradient problem inherent in standard recurrent neural networks through three gating mechanisms, input (i^t), forget (f^t), and output (o^t) gates, that regulate information flow and memory retention.

Input gate (i^t): This gate determines what portion of information from the current input is stored in the cell state. It considers both the current input and the previous cell state. The corresponding mathematical expression is:

$${\mathbf{i}}_{\mathbf{t}}=\mathsf{\sigma }\left\{{\mathbf{W}}_{\mathbf{i}}\ast \left({\mathbf{h}}_{\mathbf{t}-1}\ast {\mathbf{x}}_{\mathbf{t}}\right)+{\mathbf{b}}_{\mathbf{i}}\right\}$$

(3)

Forget gate (f^t): This gate determines what information from the previous cell state is retained and what is discarded. It considers the previous cell state and the current input. The forget gate equation is:

$${\mathbf{f}}_{\mathbf{t}}=\mathsf{\sigma }\left\{{\mathbf{W}}_{\mathbf{f}}\ast \left({\mathbf{h}}_{\mathbf{t}-1}\ast {\mathbf{x}}_{\mathbf{t}}\right)+{\mathbf{b}}_{\mathbf{f}}\right\}$$

(4)

Output gate (o^t): This gate determines what information from the current cell state is used to produce the output. It considers the current input, the previous cell state, and the information accumulated therein. The output gate expression is:

$${\mathbf{o}}_{\mathbf{t}}=\mathsf{\sigma }\left\{{\mathbf{W}}_{\mathbf{o}}\ast \left({\mathbf{h}}_{\mathbf{t}-1}\ast {\mathbf{x}}_{\mathbf{t}}\right)+{\mathbf{b}}_{\mathbf{o}}\right\}$$

(5)

Where:

i_t, f_t y o_t: vectors of the input, forget, and output gates, respectively, at time step t. $f_t{o}_t$

σ: sigmoid activation function, mapping inputs to the range [0, 1].

W_i, W_f y W_o: weight matrices of the respective gates.

h_(t−1): the concatenated previous hidden state at t − 1.

x_t: the current input vector at time t.

b_i, b_f, b_o: bias terms.

Cell state: the internal memory component of the model that retains prior information and is updated by the following equation [17]

$${\mathbf{C}}_{\mathbf{t}}={\mathbf{f}}_{\mathbf{t}}{\ast \mathbf{C}}_{\mathbf{t}-1}+{\mathbf{i}}_{\mathbf{t}}\ast \widetilde{{\mathbf{C}}_{\mathbf{t}}}$$

(6)

Where:

C_t is the cell state at time t, f_t is the forget gate, i_t is the input gate, C_t−1 is the previous cell state, and $\widetilde{{\mathbf{C}}_{\mathbf{t}}}$ is a candidate cell state computed as:

$${\mathbf{h}}_{\mathbf{t}}={\mathbf{o}}_{\mathbf{t}}\ast \mathbf{t}\mathbf{a}\mathbf{n}\mathbf{h}\left({\mathbf{C}}_{\mathbf{t}}\right)$$

(7)

where h_t is the hidden state at time step t, o_t is the output gate, and C_t is the cell state.

2.5.3. Gated Recurrent Unit (GRU)

GRU is a simplified type of recurrent neural network (RNN) inspired by LSTM, but with fewer parameters, making it highly efficient for capturing and retaining both short- and long-range dependencies in time-series data [21].

The key difference from LSTM is that GRU merges the forget and input gates into a single update gate, enabling dynamic memory control. This design mitigates the overfitting problem that can arise in LSTM when data are limited, while still addressing the long-range dependency problem. The GRU computation can be summarized through Equations (8) to (11), [17,22].

Reset gate (r^t): controls how much of the previous hidden state, h^t₋₁, is discarded. $r_t{h}_{t-1}$

$${\mathrm{r}}_{\mathrm{t}}=\mathsf{\sigma }\left\{{\mathrm{W}}_{\mathrm{r}}\left[{\mathrm{h}}_{\mathrm{t}-1},{\mathrm{x}}_{\mathrm{t}}\right]+{\mathrm{b}}_{\mathrm{r}}\right\}$$

(8)

Update gate (z^t): controls how much of the previous hidden state is carried forward to the current hidden state (h^t)$z_t)h_t$

$${\mathrm{z}}_{\mathrm{t}}=\mathsf{\sigma }\left\{{\mathrm{W}}_{\mathrm{z}}\left({\mathrm{h}}_{\mathrm{t}-1},{\mathrm{x}}_{\mathrm{t}}\right)+{\mathrm{b}}_{\mathrm{z}}\right\}$$

(9)

Candidate activation ${\tilde{h}}_t$ generates new information for the current hidden state.

$${\tilde{h}}_{\mathrm{t}}=\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left\{{\mathrm{W}}_{\mathrm{h}}\left[{\mathrm{r}}_{\mathrm{t}}\ast {\mathrm{h}}_{\mathrm{t}},{\mathrm{x}}_{\mathrm{t}}\right]+{\mathrm{b}}_{\mathrm{h}}\right\}$$

(10)

Where:

$r_t$, z_t: vectors of the reset and update gates, respectively.

${\tilde{h}}_t$: candidate activation.

σ: sigmoid activation function mapping inputs to [0, 1].

$W_r$, W_z and W_h: weight matrices applied to the input; x_t: input vector at time t.

$h_{t-1}$ h: hidden state from the previous time step.

$x_t$: current input at time t.

$b_r$, b_z, b_h: bias terms. $b_z{b}_h$

Cell state update:

$${\mathrm{h}}_{\mathrm{t}}=\left(1-{\mathrm{z}}_{\mathrm{t}}\right)\times {\mathrm{h}}_{\mathrm{t}-1}+{\mathrm{z}}_{\mathrm{t}}\times {\tilde{\mathrm{h}}}_{\mathrm{t}}$$

(11)

where:

$h_t$: Current hidden state, which is a combination of the previous hidden state, $h_{t-1}$, and the candidate activation, which is controlled by the update gate, $z_t$.

2.5.4. Discrete Wavelet Transform (DWT)

The Discrete Wavelet Transform (DWT) is a powerful signal-processing technique that decomposes a time series into multiple resolution levels, separating components with different frequency characteristics. The continuous wavelet transform can be expressed through Equations (12) to (29), [46,47]

$$wt\left(a,b\right)={\displaystyle \frac{1}{\sqrt{a}}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}x\left(t\right)\hspace{0.17em}\mathsf{\phi }\left(a,b\right)\hspace{0.17em}dt$$

(12)

where $\phi$(a, b) is the complex conjugate of the base signal, $\phi$

$$\mathsf{\phi }\left(a,b\right)={\displaystyle \frac{1}{\sqrt{a}}}\mathsf{\phi }\left({\displaystyle \frac{t-b}{a}}\right)w$$

(13)

However, CWT’s scale parameter’s dyadic discretization, which is the DWT, has a linear connection to the shift parameter’s step size.

$$\{\begin{array}{c}a=2^j\hfill \\ b=k{2}^j\hfill \end{array}$$

(14)

This transforms Equations (12) and (13) into the following expressions:

$$wt\left(j,k\right)={\displaystyle \frac{1}{\sqrt{2^j}}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}x\left(t\right)\hspace{0.17em}\mathsf{\phi }\left(j,k\right)\hspace{0.17em}dt$$

(15)

$$\mathsf{\phi }\left(j,k\right)={\displaystyle \frac{1}{\sqrt{2^j}}}\mathsf{\phi }\left({\displaystyle \frac{t-k{2}^j}{2^j}}\right)$$

(16)

Equations (15) and (16) are simplified, using the Mallat algorithm, into:

$${\mathsf{\phi }}_{\mathrm{j},\mathrm{k}}\left[\mathrm{t}\right]=2^{{\displaystyle \frac{\mathrm{j}}{2}}}\sum _{\mathrm{k}}{d}_{\mathrm{j},\mathrm{k}}\hspace{0.17em}\mathsf{\phi }\left(2^{\mathrm{j}}\mathrm{t}-\mathrm{k}\right)$$

(17)

$${\mathrm{\varnothing }}_{\mathrm{j},\mathrm{k}}\left[\mathrm{t}\right]=2^{{\displaystyle \frac{\mathrm{j}}{2}}}\sum _{\mathrm{k}}{\mathrm{a}}_{\mathrm{j},\mathrm{k}}\hspace{0.17em}\mathrm{\varnothing }\left(2^{\mathrm{j}}\mathrm{t}-\mathrm{k}\right)$$

(18)

Where:

$${\mathrm{d}}_{\mathrm{j},\mathrm{k}}=\sum _{\mathrm{k}}\mathrm{g}\left(\mathrm{n}\right)\hspace{0.17em}{\mathrm{a}}_{\mathrm{j}-1,\mathrm{k}}$$

(19)

$${\mathrm{a}}_{\mathrm{j},\mathrm{k}}=\sum _{\mathrm{k}}\mathrm{h}\left(\mathrm{n}\right)\hspace{0.17em}{\mathrm{a}}_{\mathrm{j}-1,\mathrm{k}}$$

(20)

Where ${\mathrm{d}}_{\mathrm{j},\mathrm{k}}$ are the high-frequency coefficients, ${\mathrm{a}}_{\mathrm{j},\mathrm{k}}$ are the low-frequency coefficients, $g\left(n\right)$ and $h\left(n\right)$ are the high- and low-pass filters, respectively, expressed as:

$$\{\begin{array}{c}h\left(n\right)=\langle \phi ,{\phi }_{1,n}\rangle \\ g\left(n\right)= \langle \phi ,{\phi }_{1,n}\rangle \end{array}$$

(21)

Where,

$$\sum _{n}h\left(n\right)=\sqrt{2}$$

(22)

$$\sum _{n}g\left(n\right)=0$$

(23)

For DWT with Daubechies wavelet four as the scaling function, the following matrix applies [46]:

$$\left[\begin{array}{ccc}{a}_0{a}_1{a}_2{a}_{3}0& \dots & 000\\ ⋮& \ddots & ⋮\\ a_1-a_{0}000& \dots & 0a_3-a_2\end{array}\right]\left[\begin{array}{c}{y}_1\\ ⋮\\ y_n\end{array}\right]=\left[\begin{array}{c}{s}_{j-1}\\ ⋮\\ d_{j-1}\end{array}\right]$$

(24)

Where $“a”$ are the low-frequency coefficients, $“y”$ is the raw signal, $“s”$ and $“d”$ are the transformed signal, which can also be calculated as:

$$s_{{j-1}^{\left(1\right)}}\left[n\right]=s_j\left[2n\right]+\sqrt{3}\hspace{0.17em}{s}_j\left[2n+1\right]$$

(25)

$$d_{j-1}{\left[n\right]}^{\left(1\right)}=s_j\left[2n+1\right]+{\displaystyle \frac{1}{4}}\sqrt{3}\hspace{0.17em}{s}_{j-1}\left[n\right]-{\displaystyle \frac{1}{4}}\left(\sqrt{3}-2\right)\hspace{0.17em}{s}_{j-1}^{ \left(1\right)}\left[n-1\right]$$

(26)

$$s_{j-1}{\left[n\right]}^{\left(2\right)}=s_{j-1}^{ \left(1\right)}\left[n\right]+d_{j-1}^{ \left(1\right)}\left[n+1\right]$$

(27)

$$s_{j-1}\left[n\right]={\displaystyle \frac{\sqrt{3}-1}{\sqrt{2}}}\hspace{0.17em}{s}_{j-1}^{ \left(2\right)}\left[n\right]$$

(28)

$$d_{j-1}\left[n\right]={\displaystyle \frac{\sqrt{3}+1}{\sqrt{2}}}\hspace{0.17em}{d}_{j-1}^{ \left(1\right)}\left[n\right]$$

(29)

Further details on other wavelet types can be found in [48].

2.5.5. Variational Mode Decomposition (VMD)

Variational Mode Decomposition (VMD) is a technique that decomposes a complex time series into a discrete number of intrinsic mode functions (IMFs), each characterized by a specific center frequency and limited bandwidth. VMD formulates the decomposition problem as a constrained variational optimization task, resulting in greater noise robustness, reduced mode mixing, and a stronger mathematical foundation [27]

The constrained variational problem can be expressed through Equations (30) to (35) [49],

$$\left\{{\displaystyle \sum _{k=1}^K}{\Vert {\partial }_t\left[\left(\delta \left(t\right)+{\displaystyle \frac{j}{\pi t}}\right)·u_k(t)\right]e^{j{\omega }_{k}t}\Vert }_2^2\right\}$$

(30)

$$s.t.{\displaystyle \sum _{k=1}^K}{u}_k=f\left(t\right)$$

(31)

Where $\left\{u_k\left(t\right)\right\}=\left\{u_1\left(t\right), u_2\left(t\right),\dots , u_K\left(t\right) \right\}$ and $\left\{{\omega }_k\right\}=\left\{{\omega }_1, {\omega }_2, \dots {\omega }_K\right\}$ are shorthand notations for decomposed sub-sequences and their center frequencies, respectively; $\delta \left(t\right)$ is the Dirac distribution, the symbol $\ast$ denotes convolution, and $e^{-j{\omega }_{k}t}$ is a phasor describing the rotation of the complex signal in time, with $j^2=-1$. The variational problem is addressed efficiently using the alternate direction method of multipliers (ADMMs). The modes $u_K\left(t\right)$ are updated by Wiener filtering in the Fourier domain with a filter tuned to the current center frequency; see Equation (32). Then the center frequencies ${\omega }_k$ are updated as the center of gravity of the corresponding mode’s power spectrum, expressed as Equation (33), and finally the Lagrangian multiplier $\mathsf{\lambda }$ enforcing exact constraints is updated as the dual ascent by Equation (34). The updating procedure is repeated until the convergence condition is satisfied, as in Equation (35).

$${\widehat{u}}_k^{n+1}\left(\omega \right)\leftarrow {\displaystyle \frac{\widehat{f}\left(\omega \right)-{\sum }_{i<k}{\widehat{u}}_i^{n+1}\left(\omega \right)-{\sum }_{i>k}{\widehat{u}}_i^n\left(\omega \right)+\frac{{\widehat{\lambda }}^n(\omega )}{2}}{1+2\theta {(\omega -{\omega }_k^n)}^2}}$$

(32)

$${\omega }_k^{n+1}\leftarrow {\displaystyle \frac{{\int }_0^{\infty }\omega {\left|{\widehat{u}}_i^{n+1}\left(\omega \right)\right|}^{2}d\omega }{{\int }_0^{\infty }{\left|{\widehat{u}}_i^{n+1}\left(\omega \right)\right|}^{2}d\omega }}$$

(33)

$${\widehat{\lambda }}^{n+1}\left(\omega \right)\leftarrow {\widehat{\lambda }}^n\left(\omega \right)+\tau \left(\widehat{f}\left(\omega \right)-\sum _k{\widehat{u}}_k^{n+1}\left(\omega \right)\right)$$

(34)

$$\sum _k{\displaystyle \frac{{\Vert {\widehat{u}}_k^{n+1}\left(\omega \right)\Vert }_2^2}{{\Vert {\widehat{u}}_k^n\Vert }_2^2}}<\epsilon$$

(35)

Here, ${\widehat{u}}_k^{n+1}\left(\omega \right)$, $\widehat{f}\left(\omega \right)$, and ${\widehat{\lambda }}^{n+1}\left(\omega \right)$ represent the Fourier transforms of $u_k^{n+1}\left(t\right)$, $f\left(t\right)$ and respectively; $n$ is the iteration, $\theta$ is a quadratic penalty term, $\tau$ is the iterative factor that indicates VMD’s noise tolerance, and $\epsilon$ denotes the convergence tolerance.

2.6. Model Training and Leakage Control

All experiments were implemented in Python within a Google Colab environment using TensorFlow/Keras for deep learning components and the pmdarima library for ARIMA. A global random seed was set for Python, NumPy, and TensorFlow to minimize stochastic variability and improve reproducibility. The full computational environment—including Python version, library versions, and package list—was documented for each experiment. Leakage prevention safeguards included: (i) strictly temporal, non-randomized data splits; (ii) all fitting operations (scaling, ARIMA order selection, decomposition parameterization) restricted to the training set; (iii) percentile-based event thresholds defined exclusively from training-period rainy days; and (iv) sliding-window alignment ensuring each prediction at t+1 depended solely on observations available at time t. Consistency checks verified correct date alignment, absence of missing values in model inputs, and proper separation among subsets.

2.7. Performance Evaluation

A multi-metric evaluation framework was adopted to capture different dimensions of forecast quality, recognizing that a single metric cannot adequately characterize model behavior under intermittent, heavy-tailed precipitation regimes [17,23].

2.7.1. Continuous Metrics

Model performance for precipitation amounts was quantified using four complementary metrics: Mean Absolute Error (MAE, Equation (36)), Root Mean Square Error (RMSE, Equation (37)), Nash-Sutcliffe Efficiency (NSE, Equation (38)), and Kling-Gupta Efficiency (KGE, Equations (39)–(42)).

$$\mathrm{MAE}={\displaystyle \frac{{\displaystyle \sum _{i=1}^n}|yi-{\widehat{y}}_i|}{n}}$$

(36)

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{i=1}^n}{\left(yi-{\widehat{y}}_i\right)}^2}{n}}}$$

(37)

Nash-Sutcliffe Efficiency (NSE) quantifies how well simulated and observed data fit a 1:1 line, and is evaluated using Equation (38), [31].

$$\mathrm{N}\mathrm{S}\mathrm{E}={\displaystyle \frac{{\sum }_{i=1}^n{\left(yi-{\widehat{y}}_i\right)}^2}{{\sum }_{i=1}^n{\left(yi-{\mu }_{\widehat{y}}\right)}^2}}$$

(38)

Where:

n: total number of observations in the test set.

yᵢ: observed value.

${\widehat{y}}_i$: model-predicted value.

Kling-Gupta Efficiency (KGE) is calculated using Equations (39) to (42):

$$KGE=1-\sqrt{{\left(r-1\right)}^2+{\left(\beta -1\right)}^2+{\left(\gamma -1\right)}^2}$$

(39)

Donde:

Correlation (r):

$$r={\displaystyle \frac{cov(y,\widehat{y})}{{\sigma }_y{\sigma }_{\widehat{y}}}}$$

(40)

Bias ratio (β):

$$\beta ={\displaystyle \frac{{\mu }_{\widehat{y}}}{{\mu }_y}}$$

(41)

Variability ratio (γ):

$$\gamma ={\displaystyle \frac{{CV}_{\widehat{y}}}{{CV}_y}}={\displaystyle \frac{{\sigma }_{\widehat{y}}/{\mu }_{\widehat{y}}}{{\sigma }_y/{\mu }_y}}$$

(42)

Where:

y: real observed value.

${\widehat{y}}_i$: model-predicted value.

${\mu }_y$: mean of the real variable.

${\mu }_{\widehat{y}}$: mean of the simulated variable.

${\sigma }_y:$ standard deviation of observed values.

${\sigma }_{\widehat{y}}$: standard deviation of simulated or predicted values.

2.7.2. Event-Based Categorical Metrics and Rainfall Thresholds

To assess the models’ ability to correctly detect rainfall occurrence and high-intensity events, continuous predictions were transformed into binary categorical forecasts using three precipitation thresholds: p ≥ 1.0 mm/day (rainfall occurrence), P90 (high-rainfall events), and P95 (extreme-rainfall events). The P90 and P95 thresholds were computed from rainy-day observations in the training subset and held constant for validation and test periods (Equations (43)–(45)). Three categorical metrics were computed from the resulting contingency table (Hits H, False Alarms F, Misses M):

Probability (POD):

$$POD={\displaystyle \frac{H}{H+M}}$$

(43)

False Alarm Ratio (FAR):

$$FAR={\displaystyle \frac{F}{H+F}}$$

(44)

Critical Success Index (CSI):

$$CSI={\displaystyle \frac{H}{H+M+F}}$$

(45)

POD (Probability of Detection) ranges from 0 to 1, with 1 indicating perfect event detection. FAR (False Alarm Ratio) ranges from 0 to 1, with 0 indicating no false alarms. CSI (Critical Success Index) ranges from 0 to 1 and provides an integrated measure of detection skill that penalizes both misses and false alarms [51].

3. Results and Discussion

3.1. Hydroclimatic Characterization of the Station Network

To establish a hydroclimatic baseline for the subsequent model evaluation, the rainfall regime was first characterized at each station. Table 2 presents key descriptors of the daily precipitation distributions across the 18 stations.

Precipitation regimes across the Tungurahua station network exhibited high spatiotemporal variability. Coefficients of variation (CV) exceeded 147% at all stations, confirming that all series are highly variable relative to their means—a threshold classically associated with high variability in precipitation [51,52]. CV values exceeded 200% at stations BA, CB, AP, FC, HB, JA, and PL (located between 2270 and 3314 m a.s.l.), reflecting a predominance of dry days interrupted by intense rainfall episodes.

Skewness values ranged from 2.79 (MC) to 11.61 (BA), confirming strongly positive (right-skewed) distributions at all stations. Fisher kurtosis values exceeded 10 at all stations, reaching extreme values above 100 at PL (141.08), HB (96.71), BA (250.28), and RV (80.61), indicating strongly leptokurtic distributions with heavy right tails—a well-documented feature of Andean precipitation regimes [53,54]. These distributional characteristics have critical implications for model performance, as the dominance of low and zero precipitation values creates severe class imbalance that challenges all evaluated forecasting approaches, particularly at high-intensity thresholds (Figure 3).

Daily P90 thresholds ranged from 3.20 mm (CB) to 19.83 mm (RV), while P95 thresholds ranged from 5.90 mm (CB) to 28.12 mm (RV). The highest P90/P95 values were observed at RV (82.12% rainy days), PF, and PY (83.68% rainy days), indicating that these stations experienced more frequent and severe high-intensity rainfall episodes. In contrast, CB, AP, FC, and BA showed rainy-day frequencies at or below 55%, reflecting more intermittent regimes with longer dry spells punctuated by moderate-to-intense events.

The Walsh-Lawler Seasonality Index (SI) ranged from 0.16 (QP) to 0.33 (PF), classifying most stations as ‘equable with a short dry season’ and MC, QP, and HB as ‘very equable.’ The Precipitation Concentration Index (PCI) ranged narrowly from 8.64 to 9.54, with all stations classified as having uniform monthly distribution, indicating that precipitation is not strongly concentrated in a few months of the year. Nonetheless, seasonal patterns exist: the wettest month was most frequently June, while September was most commonly the driest month. Stations exhibited unimodal (8 stations: CH, PY, ET, AS, BA, GL, PF, and—marginally—others), bimodal (9 stations), or relatively uniform (1 station: QP) annual precipitation regimes (Figure 4), consistent with the diversity of hydroclimatic typologies reported for the Ecuadorian and Colombian Andes [51,55].

The diversity of seasonal regimes has direct implications for model design and calibration. Stations with clearly defined unimodal seasonality (e.g., CH, PY, AS) may benefit from explicit seasonal components in statistical models or recurrent architectures capable of capturing annual cycles. In contrast, bimodal stations (MC, PS, CM, PL, FC, CB, HB, JA, AP, RV) require models capable of capturing multiple periodicities, potentially justifying multi-scale decomposition approaches such as DWT or VMD [26,27]. The pronounced hydroclimatic heterogeneity across the station network—with stations separated by as little as 10–30 km exhibiting fundamentally different precipitation dynamics due to orographic effects and slope orientation [55], argues strongly against spatially uniform forecasting strategies and in favor of adaptive, station-specific model selection.

3.2. Overall Continuous Predictive Performance

LSTM and GRU achieved the best overall continuous predictive performance across the station network (Figure 5), producing the lowest cross-station mean MAE and RMSE values. Among all evaluated models, LSTM emerged as the most consistently robust, combining low error magnitudes with the highest mean NSE, a mean bias close to zero, and comparatively lower cross-station performance dispersion (Figure 5e). These results indicate that LSTM most reliably reproduced the temporal behavior of the observed daily precipitation series.

This pattern is consistent with the broader hydrological deep learning literature: LSTM and GRU gating mechanisms enable the retention of relevant temporal information across varying time lags, facilitating the representation of dependencies that are structurally difficult to capture with linear models [17,45]. Several studies have highlighted this advantage, particularly when the target process is sequential, non-stationary, and exhibits complex autocorrelation structures [21,22]. For routine operational water resource management in Tungurahua Province, LSTM is therefore the most robust option for continuous daily rainfall prediction, while GRU emerges as a well-balanced alternative offering comparable performance at lower computational cost—an important consideration for resource-constrained operational contexts.

ARIMA showed the weakest overall performance among the five models, with the highest mean RMSE and MAE and the lowest NSE (Figure 5), reflecting the inability of its linear structure to capture the complex, nonlinear daily precipitation dynamics under Andean mountain regimes, particularly at stations with bimodal seasonality [10,17]. This weaker continuous prediction skill is consistent with prior comparative studies in complex rainfall environments [20,24].

The hybrid models DWT-LSTM and VMD-GRU produced intermediate error levels but negative NSE and KGE values, indicating weaker overall agreement with observations than the standalone recurrent architectures. The failure of these hybrid models to outperform simpler deep learning architectures suggests that, for daily precipitation series of approximately 12 years with high intermittency, the additional complexity of multi-scale decomposition does not translate into predictive improvement. This may be attributed to error accumulation during component reconstruction, fragmentation of temporal information into sub-series with insufficient training data for individual neural networks, and the tendency of VMD to generate spurious modes under limited data conditions [31,45,56].

3.3. Station-Wise Variability in Model Performance

The network-scale performance hierarchy described above was not spatially homogeneous, as revealed by the station-wise ranking analysis (Figure 6). LSTM most frequently ranked first or second in MAE and RMSE across the 18 stations, consistent with its superior network-level performance. However, GRU and DWT-LSTM outperformed LSTM at several individual stations, confirming that no single universally optimal model exists and that the best-performing architecture depends on the local hydroclimatic context [10,17].

Stations where GRU outperformed LSTM likely feature precipitation patterns with shorter effective temporal dependencies, a context in which GRU’s simplified architecture with fewer parameters provides better generalization [21,56]. The occasional superiority of DWT-LSTM in NSE and KGE at certain stations suggests that wavelet decomposition can improve representation of multi-scale precipitation variability, particularly at stations with complex patterns combining long-term seasonal trends and high-frequency fluctuations [26,28].

Bias behavior also varied across stations and models (Figure 6e). LSTM showed nearly balanced bias, GRU showed slight underestimation (attributable to its single update-gate structure, which can dampen extreme event magnitudes [21]), DWT-LSTM showed overestimation (likely arising from reconstruction errors when high-frequency sub-series predictions are recombined [26,28]), and ARIMA showed moderate overestimation. The integrated ranking (Figure 6f) confirmed LSTM as the most consistently competitive model, followed by GRU, with ARIMA generally remaining among the weaker models and VMD-GRU showing limited cross-site stability.

These findings have direct implications for operational provincial forecasting networks: spatially homogeneous strategies applying a single model to all stations are suboptimal and likely to underperform at individual stations. An adaptive multi-model system—applying LSTM at bimodal stations with complex long-range dependencies (e.g., MC, PS, RV), GRU at stations with simpler unimodal regimes where parsimony improves generalization (e.g., CH, QP), and DWT-LSTM at stations with pronounced multi-scale variability (e.g., PL, BA)—can improve aggregate network accuracy while maintaining computational efficiency.

3.4. Rainfall Event Detection at the 1 mm/Day Threshold

Accurate estimation of rainfall magnitude does not necessarily translate into reliable detection of rainy-day occurrence. At the 1 mm/day threshold, LSTM, GRU, DWT-LSTM, and VMD-GRU achieved high POD values (Figure 7a), reflecting the capacity of deep learning models to capture the nonlinear atmospheric patterns preceding rainfall onset [17,25]. In contrast, ARIMA showed lower and more variable POD, consistent with its linear structure’s reduced sensitivity to the abrupt transitions characteristic of precipitation intermittency [12,18].

Despite their high sensitivity, neural network models—particularly VMD-GRU—produced elevated FAR values (Figure 7b), while ARIMA exhibited the lowest FAR. Neural networks trained by minimizing mean squared error (MSE) tend to generate small positive predictions under uncertain conditions, which, when compared against a 1 mm binary threshold, results in false alarms when observed precipitation is zero [23,31]. ARIMA’s autoregressive structure relies on the persistence of dry conditions and is therefore less prone to predicting spurious wet-to-dry transitions, though at the cost of missing real events.

When the POD-FAR trade-off was integrated through CSI (Figure 7c), ARIMA showed better median overall detection skill at the 1 mm threshold, while GRU and LSTM remained moderately competitive. This result—apparently contradictory to the superior continuous-metric performance of deep learning models—illustrates an important methodological limitation: optimizing for magnitude accuracy (MSE minimization) does not guarantee skill in binary occurrence classification. The trade-off between sensitivity and specificity has critical implications that depend on the operational context: in flood early warning and reservoir management, high POD is prioritized to avoid missing precipitation events that could trigger hazardous conditions, and a moderate FAR is accepted as an unavoidable operational cost [8,32].

3.5. Detection of High- and Extreme-Rainfall Events

Model skill deteriorated dramatically when the detection threshold was elevated from 1 mm/day to the P90 and P95 percentiles (Figure 8). At the P90 threshold, POD values were close to zero for most models across the 18 stations, while FAR values remained high (medians at or above 0.8–1.0 for GRU, DWT-LSTM, ARIMA, and VMD-GRU), representing the worst possible scenario: models fail to detect real high-intensity events while simultaneously generating abundant false alarms when predicted magnitudes exceed the P90 threshold without corresponding to observed extreme precipitation.

CSI values were very low across the station network for all models at P90. ARIMA showed slightly higher POD and CSI than the deep learning architectures—a result that apparently contradicts the superior continuous-metric performance of LSTM and GRU, and challenges the assumption that greater architectural complexity necessarily yields better extreme-event predictive capacity. ARIMA’s relative advantage at high thresholds can be explained by its autoregressive structure: under conditions of temporal persistence or clustering in the upper tail of the precipitation distribution, the linear autoregressive components can partially capture serial correlation among consecutive intense events [1,13]. Furthermore, by operating directly on the original precipitation series without the intermediate decomposition and reconstruction steps required by hybrid models, ARIMA avoids the cumulative errors associated with those processes [26,27].

The overall collapse in predictive skill was even more severe at the P95 threshold, where POD (Figure 8d) and CSI (Figure 8f) were practically zero for all models across the station network, and FAR values approached 1.0. ARIMA showed the only marginally detectable deviation from zero POD. This catastrophic failure at P95 reflects multiple inherent challenges that univariate approaches cannot adequately resolve. The 95th percentile defines the upper 5% of the daily precipitation distribution, corresponding to approximately 18 extreme events per station per year—yielding roughly 216 extreme event examples over the 12-year training period. This severe class imbalance between normal and extreme days fundamentally limits the ability of any supervised learning algorithm trained with standard loss functions to learn robust representations of extreme-event precursor conditions [29,53]. Neural network models optimized by MSE minimization naturally prioritize accurate prediction of frequent low-to-moderate precipitation events at the expense of rare extreme episodes [17,23].

Furthermore, extreme precipitation typically results from complex atmospheric configurations—involving interactions among synoptic-scale dynamics, mesoscale convective systems, and local orographic forcing—that cannot be inferred from precipitation history alone. Incorporating additional meteorological predictors (atmospheric humidity profiles, wind fields, sea surface temperature anomalies, large-scale circulation indices) has been shown to significantly improve extreme-event detection in complex terrain [20,57]. These results underscore the fundamental inadequacy of purely univariate approaches for operationally reliable extreme-event rainfall forecasting in the Ecuadorian Andes.

3.6. Integrated Cross-Domain Synthesis of Model Performance

Table 3 presents an integrated synthesis of model performance across all evaluation domains. No single model dominated across all forecasting objectives, and the optimal operational approach depends on the specific application context and station characteristics.

LSTM’s leadership in continuous prediction is attributable to its sophisticated memory cell architecture, with gating mechanisms that capture long-range nonlinear temporal dependencies—a fundamental characteristic for modeling seasonal persistence and autocorrelation in precipitation series [32,45]. GRU ranked second and demonstrated the most balanced behavior across all evaluation domains, making it particularly attractive for provincial-scale operational networks where a single model must perform acceptably across multiple forecasting objectives simultaneously [21].

DWT-LSTM showed an intermediate profile: weaker baseline performance for continuous prediction (3rd, integrated rank 3.06 ± 1.04) and 1 mm detection (4th), but comparatively better rankings for P90 and P95 detection (2nd in both). This differential behavior suggests that wavelet decomposition enables the LSTM component associated with high-frequency detail sub-series to capture abrupt intensity peaks characteristic of extreme events, even when the overall continuous metric performance is penalized by reconstruction errors [26,28,31].

ARIMA’s paradoxical first-place ranking for extreme-event detection (P90 and P95) relative to the deep learning architectures should be interpreted cautiously: as demonstrated by near-zero POD and FAR ≥ 0.8 at both thresholds across the station network (Figure 8), this relative superiority reflects not genuine predictive skill but rather the conservative prediction behavior of the autoregressive structure, which tends to revert toward historical means and thereby reduces the probability of spurious extreme-threshold crossings at the cost of near-zero detection rates. ARIMA’s ranking advantage at extreme thresholds therefore represents relative rather than absolute superiority [13,18]. VMD-GRU showed the lowest overall robustness consistently across all evaluation domains, occupying last or penultimate positions, and failed to show local compensating advantages at any individual station. VMD’s requirement for a priori specification of the number of intrinsic modes K may have led to spurious mode generation under limited and highly intermittent data conditions, fragmenting the precipitation signal into components without hydroclimatic significance [27,49].

3.7. Limitations and Future Research Directions

This study acknowledges several important limitations that constrain the generalizability of the findings and motivate future research. First, the approximately 12-year record length provides an insufficient number of P95 extreme-event examples for robust neural network training under standard loss functions, contributing to the near-zero detection skill observed at extreme thresholds. Second, the strictly univariate framework—employing daily precipitation as the sole predictive variable—excludes complementary meteorological predictors (atmospheric humidity, wind fields, temperature, large-scale circulation indices) that have been demonstrated to substantially improve extreme-event prediction [20,57]. Third, evaluation was confined to the one-step-ahead (t + 1) horizon; performance degradation at operationally relevant extended horizons (t + 3 to t + 14 days) was not assessed.

Future research should prioritize: (i) multivariate integration of atmospheric precursor variables to improve P90/P95 extreme-event detection; (ii) implementation of event-sensitive loss functions—including weighted MSE, quantile regression loss, and focal loss—that explicitly up-weight rare extreme events during training; (iii) probabilistic and ensemble forecasting frameworks that provide calibrated uncertainty estimates; (iv) multi-horizon evaluation (1–14 days) to characterize temporal performance degradation; and (v) transfer learning or data-augmentation strategies to overcome the limited extreme-event sample size. The application of physics-informed neural networks that incorporate hydrological process knowledge as structural constraints represents a particularly promising research direction for complex mountain environments.

4. Conclusions

This study provided a systematic and methodologically rigorous comparison of five univariate daily rainfall forecasting models—ARIMA, LSTM, GRU, DWT-LSTM, and VMD-GRU—across 18 meteorological stations in Tungurahua Province, Ecuador (2013–2025), using a comprehensive evaluation framework integrating continuous accuracy metrics (MAE, RMSE, NSE, KGE) and categorical event-detection metrics (POD, FAR, CSI) at three precipitation thresholds (1 mm/day, P90, P95). The following principal conclusions are drawn:

LSTM achieved the best overall performance for continuous daily rainfall prediction and rainfall-occurrence detection at the 1 mm/day threshold, while GRU delivered comparable skill at substantially lower computational cost. Both deep learning architectures significantly outperformed ARIMA in continuous accuracy metrics, demonstrating the superiority of nonlinear recurrent architectures over classical linear statistical models for daily rainfall prediction in topographically complex Andean environments.

Model performance was not spatially uniform: station-wise rankings varied considerably across the 18-station network, reflecting the strong hydroclimatic heterogeneity of the Ecuadorian Andean region. No single model ranked first at all stations, confirming that spatially uniform single-model forecasting strategies are suboptimal and that adaptive multi-model systems tailored to local hydroclimatic typology are warranted.

All evaluated models failed to provide operationally reliable detection of high-intensity (P90) and extreme (P95) rainfall events, producing near-zero POD and FAR values at or above 0.8 across the station network. Although ARIMA achieved the best relative ranking for extreme-event detection, its absolute predictive capacity was insufficient for operational early warning or reservoir management. This universal failure at high-intensity thresholds underscores the fundamental inadequacy of univariate approaches—regardless of architectural complexity—for reliably forecasting rare extreme precipitation events from precipitation history alone.

Hybrid decomposition models (DWT-LSTM, VMD-GRU) did not outperform the simpler deep learning architectures in continuous prediction; VMD-GRU showed the lowest overall robustness across all evaluation domains. These results suggest that, for daily precipitation series of approximately 12 years with high intermittency, multi-scale decomposition does not provide a net predictive benefit and may exacerbate error accumulation during reconstruction.

Future operational rainfall forecasting systems in the Ecuadorian Andes should prioritize multivariate predictor frameworks incorporating atmospheric precursor variables, event-sensitive loss functions, and probabilistic uncertainty quantification to address the critical limitations identified in this study, particularly for extreme-event prediction at operationally relevant thresholds.