Evaluating Rainfall Timing and Volume of Gridded Precipitation Products

1. Introduction

The accurate estimation of precipitation is a fundamental prerequisite for hydrological modelling, agricultural forecasting, and flood risk management [1,2]. In East Africa, where the spatiotemporal variability of rainfall is heavily influenced by the interplay of the Intertropical Convergence Zone [3], complex topography, and local convective instability, the scarcity of ground-based gauge networks presents a bottleneck [4]. Consequently, hydrologists increasingly rely on GPPs, including satellite-based estimates such as IMERG [5] and CHIRPS [6], and reanalysis products like ERA5 [7], to serve as input for hydrological models [8].

However, the utility of these GPPs is often limited by systematic discrepancies when compared to ground observations. Traditional validation studies employ point-to-grid comparisons using error metrics such as Root Mean Square Error (RMSE), Mean Bias Error (MBE), and the Kling-Gupta Efficiency (KGE) on raw precipitation data [9,10]. While traditional metrics provide a necessary assessment of volume accuracy, they are often insufficient for diagnosing the underlying mechanisms of a model’s failure. For instance, a dataset may achieve a satisfactory daily volume by simulating persistent, low-intensity drizzle rather than the actual observed short-duration, high-intensity rain events [11,12,13]. On the other hand, a high-resolution satellite might capture a physical rainfall event correctly, but standard metrics will score it as a total failure simply because the rain event was recorded slightly in the wrong place or at the wrong time [14].

This rigidity gives rise to the “double penalty” problem—where a GPP dataset is penalised once for missing the exact timing or location of a rain event, and again for predicting the wrong intensity, hiding the dataset’s true skill [15]. Research in forecast validation has suggested moving beyond rigid point-to-point matching toward fuzzy, neighbourhood-based evaluations [16] and pattern-based evaluations [17]. By evaluating predictions within a spatial or temporal window of tolerance, these neighbourhood approaches disentangle timing accuracy from volume accuracy. Distinguishing whether a dataset fails due to a temporal phase misalignment or simply because it suffers from volumetric attenuation is critical for developing effective bias-correction strategies [18].

In this study, we perform an analysis of four widely used precipitation products: ERA5, IMERG, TAMSAT, and CHIRPS, against the Trans-African Hydro-Meteorological Observatory (TAHMO) network, which at the time of writing, has 332 stations across the East African region [19]. While recent studies have successfully leveraged the TAHMO network to expose volumetric biases in satellite products [20], these evaluations have predominantly relied on standard, raw-space linear metrics aggregated at daily or monthly scales. Consequently, the sub-daily physical mechanisms driving these errors often remain obscured. To address this, our methodology is two-fold. First, we address the extreme skewness of daily precipitation distributions by applying variance-stabilising transformations, enabling the assessment of signal consistency within a variance-stabilised “transformed space” of the data [21,22]. Second, we evaluate the timing consistency of these signals. By partitioning the data into 30-day overlapping windows, we apply a time-shifting technique combined with a Gaussian temporal filter. This allows us to identify the best possible timing alignment for the rainfall events within each specific temporal window.

This approach allows us to answer two critical questions regarding the observed differences between datasets: First, are these deviations driven by a fundamental lack of physical signal correlation, or by a systematic temporal misalignment? Second, to what extent do the GPPs and point-scale observations differ regarding the diurnal precipitation cycle in East Africa? By extending our analysis to the diurnal scale, we aim to identify the structural mechanisms that influence the performance characteristics of GPPs. The findings presented here offer a perspective on GPP selection, suggesting that the selection of the “appropriate” GPP depends on the use-case, considering and prioritising either the accurate timing and frequency of precipitation events [23] or the conservation of total rainfall volume required for processes such as downscaling [24].

The rest of the paper is organised as follows: Section 2 describes the data and the methods used in the analysis to obtain the results in Section 3. In Section 4 we discuss the results. Section 5 concludes the paper.

2. Materials and Methods

2.1. Study Area and Data Acquisition

The study domain encompasses the East African region, covering Kenya, Uganda, Tanzania, Rwanda, and Burundi (Figure 1) between 2018 to 2024. This region is positioned approximately between latitudes 5°N and 12°S, and longitudes 29°E and 42°E; it is characterised by highly complex topography that heavily influences local weather patterns. Elevations range from sea level along the eastern Indian Ocean coastline to regions exceeding 5,000 meters above sea level (such as Mount Kenya and Mount Kilimanjaro). Furthermore, features like the Great Rift Valley and Lake Victoria create distinct, localised microclimates [25,26]. This topographical variation drives the complex, localised rainfall patterns—and the highly varied diurnal precipitation cycles—that make precipitation evaluation in this region particularly challenging.

2.1.1. Station Network (TAHMO)

Reference precipitation data was obtained from the Trans-African Hydro-Meteorological Observatory (TAHMO) network as shown in Figure 1. TAHMO data is available at a 5-minute resolution, aggregated to both hourly and daily for our analyses [19].

2.1.2. Gridded Precipitation Products (GPPs)

We evaluated four widely used GPPs representing different estimation methodologies: Reanalysis (ERA5), Satellite-Infrared (TAMSAT), Satellite-Microwave merged (IMERG), and Blended Satellite-Gauge (CHIRPS). Table 1 summarises the product characteristics.

2.2. Data Pre-Processing and Quality Control

2.2.1. Station Selection and Temporal Alignment

To ensure robust statistical comparisons, a filtering protocol was applied to the TAHMO network:

• Daily Completeness: For daily aggregation, a station was considered valid only if at least 21 out of 24 hourly records were present.
• Series Continuity and hourly completeness: Stations missing more than 30% of hourly data over the 7-year study period were excluded.
• Temporal Intersection: To eliminate sampling bias, we obtained the intersection of valid timestamps across all datasets and only retained timestamps where the TAHMO station and the GPPs had valid data.

This filtering process yielded 117 stations meeting the daily completeness criteria and 122 stations meeting the series continuity requirements. To maintain a consistent spatial footprint across all evaluations, the final observation network was restricted to the subset of 117 stations. Prior to evaluation, the datasets were aggregated into daily totals for the rolling windows analysis, and the sub-daily datasets (IMERG, ERA5 and TAHMO) into hourly intervals to evaluate the diurnal cycle.

2.3. Methodology

The sections below discuss the methodological framework that we applied in our analysis.

2.3.1. Variance Stabilisation and Transformed Space

Daily precipitation data exhibits extreme positive skewness and zero-inflation. To address this, we transformed all data into a variance-stabilised “transformed space” [21] using two methods:

• Box-Cox Transformation: Defined by Equation (1) below:

  $$y\left(\lambda \right)=\left\{\begin{array}{cc}{\displaystyle \frac{{(x+ϵ)}^{\lambda }-1}{\lambda }}\hfill & \mathrm{if}\lambda \ne 0\hfill \\ ln(x+ϵ)\hfill & \mathrm{if}\lambda =0\hfill \end{array}\right.$$

  (1)

  where:

  –

  $y\left(\lambda \right)$ represents the transformed precipitation value,

  –

  x is the raw precipitation intensity

  –

  $ϵ$ is an optimised offset derived from the data to accommodate zero-inflation

  –

  $\lambda$ is the shape parameter governing the severity of the transformation [28].

  These parameters are tuned independently; this effectively projects each GPP–TAHMO pair into a distinct, independent transformed space that accommodates the unique distributional characteristics (e.g., specific zero-inflation or tail behavior) of that product.
• Inverse Hyperbolic Sine (ArcSinh) Transformation: We applied the Inverse Hyperbolic Sine transformation as an alternative to the Box-Cox transformation to handle extreme values [22], as defined by Equation (2) below:

  $$y=\mathrm{asinh}\left(x\right)=ln\left(x+\sqrt{x^2+1}\right)$$

  (2)

  where:

  –

  y represents the transformed precipitation value

  –

  x is the raw precipitation intensity.

To maintain mathematical consistency for pairwise comparison, the optimal transformation parameters ($\lambda$ and offset $ϵ$) derived from a specific GPP were subsequently applied to the corresponding TAHMO station series. If the ground observations and the gridded predictions were transformed independently using their own optimal parameters, the two signals would be projected into entirely different dimensional spaces. This scale mismatch would render relative amplitude comparisons—such as the Kling-Gupta Bias and Variability ratios mathematically invalid.

Instead, our methodology projects TAHMO data into the specific variance-stabilised “transformed space” of each individual GPP. For example, in the Box-Cox space optimised for ERA5, physical rainfall thresholds (e.g., $T_{phys}=2.5\mathrm{mm}/\mathrm{day}$) are translated into an exact variance-stabilised equivalent ($T_{trans}$) using ERA5’s specific $\lambda$ and $ϵ$ parameters. Both the TAHMO and ERA5 signals are then evaluated against this shared $T_{trans}$ threshold.

This approach ensures that the evaluation metrics consider true phase dynamics while preserving the distinction between wet and dry states. Because GPPs parameterise dry states distinctly—such as ERA5 exhibiting a persistent trace drizzle bias [11], mapping TAHMO data into the GPP’s native transformed space forces the statistical evaluation to test whether the GPP correctly times its own version of a wet event against the physical occurrence of rain.

2.3.2. Spatial Interpolation, Variance Stabilisation, and Kernel Optimisation

To extract gridded data at specific station coordinates, we compared Nearest Neighbour interpolation against Bilinear Interpolation as a preliminary step.

This preliminary evaluation was conducted on the raw, unsmoothed daily time series ($\sigma =0$) prior to any variance stabilisation, ensuring that temporal smoothing did not blur the interpolation artifacts. For each station, we computed the Pearson correlation coefficient (r) between the ground gauge and the GPP signal. We then calculated the arithmetic mean of the station-level correlation differences ($\Delta r=r_{bilinear}-r_{nearest}$) across the observational network.

We also compared the Pearson correlation values of each GPP-TAHMO pair at various kernel widths for the Gaussian temporal smoothing (ranging from 0 to 21 days) to determine the optimal kernel width.

We then evaluated the correlation of the entire historical time series without applying temporal rolling windows. The Pearson correlation coefficient (r) was calculated independently at every valid TAHMO station to capture localised spatial dynamics. To quantify the overall efficacy of a transformation method (e.g., Box-Cox, ArcSinh), we aggregated these station-level correlations to construct a network-wide performance distribution. The most appropriate transformation for each GPP was then determined by evaluating the central tendency (median and arithmetic mean) and the interquartile spread of these station-level correlations. This ensured that the selected mathematical transformation structurally preserved or enhanced the baseline temporal signal across the majority of the spatial domain prior to conducting the localised phase-shift analysis.

2.3.3. Gaussian Temporal Smoothing

To analyze the consistency of large-scale weather patterns and dampen the effect of minor timing mismatches, we applied a 1D Gaussian smoothing kernel to all time series as shown by Equation (3) below:

$$S_{smooth}\left(t\right)={\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}\sum _{\tau }{S}_{raw}\left(\tau \right)exp\left(-{\displaystyle \frac{{(t-\tau )}^2}{2{\sigma }^2}}\right)$$

(3)

where:

• $S_{smooth}\left(t\right)$ is the temporally smoothed precipitation value at the target day t,
• $S_{raw}\left(\tau \right)$ is the raw, unshifted precipitation value at day $\tau$,
• $\sigma$ is the standard deviation of the Gaussian kernel, governing the temporal spread of the smoothing filter,
• $\tau$ represents the temporal indices within the localised rolling window surrounding t.

We used this to compare the Pearson correlation values of each GPP-TAHMO pair at various kernel widths for the Gaussian temporal smoothing (ranging from 0 to 21 days) to determine the optimal kernel width.

2.3.4. Rolling Window Cross-Correlation Analysis

To diagnose non-stationary temporal misalignments, we performed a rolling window analysis on the transformed data. The time series was partitioned into 30-day windows (total $N=362$ windows). Within each window, we:

• Applied a temporal shift $\Delta t$ to the GPP signal ranging from $[-3,+3]$ days.
• Calculated the similarity between the shifted GPP signal and the unshifted TAHMO signal.
• Identified the optimal shift that maximised the performance metric.

Windows exhibiting zero variance (for example all zeros) were flagged and excluded from the correlation analysis to prevent numerical instability.

2.4. Performance Metrics and Scoring

To evaluate the spatiotemporal coherence of the GPPs, we employed a scoring framework that prioritises phase dynamics and event detection capability over simple amplitude matching. This approach mitigates the “double penalty” problem inherent in point-to-grid comparisons, particularly within the transformed space.

2.4.1. Window-Level Optimisation and Event Detection

The fundamental unit of analysis was the 30-day time window ($C_k$). For each window k, station s, and dataset d, we evaluated the optimal temporal alignment between each GPP–TAHMO pair by maximising the cross-correlation function over a temporal lag window $\tau \in [-3,+3]$ days.

To evaluate the underlying phase of the atmospheric signal, the transformed GPP time series was first shifted by lag $\tau$, after which both the baseline transformed TAHMO signal and the shifted GPP signal were processed concurrently through a Gaussian smoothing filter. The optimal Pearson Correlation ($r_{opt}$) is calculated on these smoothed continuous signals using Equation (4) below:

$$r_{opt}(s,k)=\underset{\tau }{max}\left[{\displaystyle \frac{{\sum }_{t\in C_k}(O_{sm}\left(t\right)-{\overline{O}}_{sm})(M_{sm}(t+\tau )-{\overline{M}}_{sm})}{\sqrt{{\sum }_{t\in C_k}{(O_{sm}\left(t\right)-{\overline{O}}_{sm})}^2{\sum }_{t\in C_k}{(M_{sm}(t+\tau )-{\overline{M}}_{sm})}^2}}}\right]$$

(4)

where:

• $r_{opt}(s,k)$ is the optimal correlation achieved at station s during window k
• $\tau$ is the applied temporal lag, $C_k$ is the set of temporal indices within window k
• $O_{sm}\left(t\right)$ represents the transformed and smoothed TAHMO signal at time t
• $M_{sm}(t+\tau )$ represents the transformed and smoothed GPP signal at time t offset by the temporal lag $\tau$
• ${\overline{O}}_{sm}$ represents the arithmetic mean of $O_{sm}$ over the window
• ${\overline{M}}_{sm}$ represents the arithmetic mean of $M_{sm}(t+\tau )$ over the window.

Event Detection

Evaluating rain event detection on smoothed data artificially flattens sharp, high-intensity convective peaks, causing them to fall below detection thresholds and unfairly penalise datasets. To address this, we evaluated event detection using the Critical Success Index (CSI) [29] via a two-tiered approach: temporal point-to-point CSI and temporal neighbourhood (fuzzy) CSI.

A physical rainfall threshold of $T_{phys}=2.5\mathrm{mm}/\mathrm{day}$ was translated into the equivalent variance-stabilised threshold ($T_{trans}$) using the respective Box-Cox or ArcSinh parameters.

First, the point-to-point ($CS{I}_{point}$) was computed on the unsmoothed transformed data to evaluate exact temporal precision using Equation (5):

$$CS{I}_{point}(s,k)={\displaystyle \frac{H_{raw}}{H_{raw}+F_{raw}+Z_{raw}}}$$

(5)

where:

• $H_{raw}$ are hits
• $F_{raw}$ are false alarms
• $Z_{raw}$ are misses,

which was calculated by evaluating the unshifted raw signals directly against $T_{trans}$.

Second, we computed a fuzzy Critical Success Index ($CS{I}_{\mathrm{fuzzy}}$) by transitioning from strict point-to-point matching to a temporal neighbourhood approach. While this metric retains the exact formulation of the point-to-point CSI, the data processing pipeline proceeded in three sequential steps:

1.

Binarisation: Transformed precipitation signals were converted to binary events, where daily totals exceeding the physical threshold ($T_{\mathrm{trans}}$) equalled 1.0, and values below equalled 0.0.

2.

Gaussian Smoothing: A Gaussian filter ($\sigma =2$ days) was applied to the binary vectors. This transformed isolated, discrete rain days into continuous probability curves bounded within $[0,1]$.

3.

Threshold Evaluation: These smoothed curves were evaluated against a strict probability threshold to derive the resulting hits, misses, and false alarms.

For an isolated precipitation event recorded by TAHMO at day $t_0$, the continuous probability profile $P(\Delta t)$ at a temporal shift of $\Delta t=|t-t_0|$ days is defined by the Gaussian kernel function in Equation (6):

$$P(\Delta t)={\displaystyle \frac{1}{\sigma \sqrt{2\pi }}}exp\left(-{\displaystyle \frac{\Delta t^2}{2{\sigma }^2}}\right)$$

(6)

Using a kernel width of $\sigma =2$ days, setting the threshold at $P\ge 0.15$ creates an operational window to evaluate how closely the GPPs align with TAHMO:

• Exact Alignment (0-day shift): A direct daily match ($\Delta t=0$) yields a peak probability of $P\approx 0.199$, successfully registering as a Hit.
• Acceptable Shift ($\pm 1$-day shift): A 1-day shift ($\Delta t=1$) yields $P\approx 0.176$. Because this stays above the $0.15$ boundary, the dataset is rewarded for maintaining close temporal consistency.
• Unacceptable Shift ($\pm 2$-day shift): A 2-day shift ($\Delta t=2$) yields $P\approx 0.121$. Falling below the $0.15$ threshold, this divergence is penalised as a Miss.

By positioning the threshold between the 1-day and 2-day discrete intervals—which equates to a continuous temporal radius of $\Delta t=1.51$ days—the framework establishes a cutoff for daily observations. This selection is motivated by the operational requirements of the East African domain: it forgives minor 1-day physical timing shifts caused by convective phase lags, while rejecting displacements of 2 days or more that would be unsuitable for agricultural planning or early flood warnings.

Physical Scaling and Variability

We computed the Kling-Gupta Efficiency (KGE) components for each window on the smoothed signals to diagnose intensity and variability errors using ratios as shown in Equation (7):

$${\beta }_k={\displaystyle \frac{{\mu }_M}{{\mu }_O}}(\mathrm{Bias}\mathrm{Ratio}),{\alpha }_k={\displaystyle \frac{{\sigma }_M}{{\sigma }_O}}(\mathrm{Variability}\mathrm{Ratio})$$

(7)

where:

• ${\beta }_k$ is the volumetric bias ratio for the specific temporal window k
• ${\alpha }_k$ is the variability ratio for the specific temporal window k
• ${\mu }_M$ and ${\mu }_O$ represent the arithmetic means of the GPP and TAHMO precipitation signals, respectively, calculated over the duration of window k in the transformed space
• ${\sigma }_M$ and ${\sigma }_O$ denote the standard deviations of the GPP and TAHMO precipitation signals, respectively, within window k.

Because these statistics are computed on variance-stabilised data, their physical interpretation differs from standard metrics. ${\beta }_k$ quantifies the shift in the Typical Intensity State (e.g., indicating whether a GPP’s baseline behavior simulates physical dry days or parameterised drizzle), while ${\alpha }_k$ measures the signal’s ability to reproduce the full span of observed weather fluctuations.

2.4.2. Station-Level Scoring

To aggregate performance across all temporal windows for a single station, we defined a bounded Composite Score ($S_{comp}$). Aggregating the Kling-Gupta Bias ($\beta$) and Variability ($\alpha$) ratios presents two distinct mathematical challenges: ratio asymmetry and state mismatches.

First, simple arithmetic averaging of ratios introduces asymmetry errors. For example, a GPP that overestimates by double ($\beta =2.0$) in one window and subsequently underestimates by half ($\beta =0.5$) in the next yields an arithmetic mean of 1.25. This falsely implies a systematic 25% overestimation, despite the GPP being balanced overall.

Second, within the variance-stabilised transformed space, datasets may disagree on the weather state. For instance, a GPP might register a precipitation event (yielding a positive transformed value), while the gauge registers physical dry noise (yielding a negative transformed value). This “state mismatch” results in a negative ${\beta }_k$ ratio, which is physically uninterpretable for volumetric scaling and mathematically undefined for logarithmic transformation. Furthermore, the CSI scores already capture rain event detection, effectively penalising these state mismatches.

To address these issues, we first filtered the evaluation to isolate valid windows (V) where both datasets agreed on the overall distributional state (${\beta }_k>0$). The true aggregate physical bias was then calculated using the geometric mean of these valid ratios, preventing arithmetic asymmetry, as defined in Equation (8):

$${\mu }_{\beta }\left(s\right)=exp\left({\displaystyle \frac{1}{N_{valid}}}\sum _{k\in V}ln\left({\beta }_k\right)\right)$$

(8)

Because the final Composite Score requires evaluation metrics to operate on a normalised $[0,1]$ scale, these unbounded physical ratios (${\mu }_{\beta }$ and ${\mu }_{\alpha }$) must be scaled without introducing directional bias. We achieved this by applying a symmetric exponential decay function to the logarithmic space of the geometric means. As defined in Equation (9), this guarantees that proportional overestimations and underestimations decay toward zero at the exact same penalty rate (e.g., $\beta =2.0$ and $\beta =0.5$ both yield a symmetric score of $0.5$), bounding the final scores $S_{\beta }$ and $S_{\alpha }$ between $[0,1]$:

$$S_{\beta }\left(s\right)=exp\left(-\left|ln\left({\mu }_{\beta }\left(s\right)\right)\right|\right),S_{\alpha }\left(s\right)=exp\left(-\left|ln\left({\mu }_{\alpha }\left(s\right)\right)\right|\right)$$

(9)

Final Composite Formula

The final station-level Composite Score is a weighted linear combination designed to prioritise phase dynamics [30] over easily corrected volumetric amplitude errors as shown by Equation (10):

$$S_{comp}\left(s\right)=w_r{\overline{r}}_{opt}+w_{fuzzy}{\overline{CSI}}_{fuzzy}+w_{point}{\overline{CSI}}_{point}+w_{cons}{C}_{cons}+w_{\beta }{S}_{\beta }+w_{\alpha }{S}_{\alpha }+w_{imp}{S}_{imp}$$

(10)

where the weights reflect the prioritisation of phase alignment:

• Phase dynamics ($w_r=0.50$): Based on mean optimal correlation (${\overline{r}}_{opt}$).
• Fuzzy event detection ($w_{fuzzy}=0.15$): Based on the temporal neighbourhood event detection (${\overline{CSI}}_{flex}$).
• Point-to-point detection ($w_{point}=0.10$): Based on exact phase alignment (${\overline{CSI}}_{point}$).
• Stability ($w_{cons}=0.10$): Based on the inverse standard deviation of the optimal correlations across all windows ($C_{cons}$).
• Volume accuracy ($w_{\beta }=0.05,w_{\alpha }=0.05$): The normalised mean log bias ($S_{\beta }$) and log variability ($S_{\alpha }$) scores.
• Correlation Improvement ($w_{imp}=0.05$): Rewarding correlation improvements achieved through temporal shifting.

The final dataset ranking ($S_{global}$) was derived from the spatial arithmetic mean of the Composite Score across all valid ground stations.

2.4.3. Dataset-Level Ranking

The final dataset ranking is derived from the spatial arithmetic mean of the Composite Score across all M valid stations. Since the use of bilinear interpolation explicitly accounts for sub-grid spatial variability by reconstructing the continuous precipitation field at the exact station coordinates, no secondary distance-based weighting was required. We used Equation (11) to rank the datasets:

$$S_{global}\left(\mathcal{D}\right)={\displaystyle \frac{1}{M}}{\displaystyle \sum _{s=1}^M}{S}_{comp}(s,\mathcal{D})$$

(11)

where:

• $S_{global}\left(d\right)$ is the final global performance ranking for a specific GPP (dataset $\mathcal{D}$),
• M is the total number of valid spatial ground stations within the evaluated network,
• s is the index representing an individual ground station,
• $S_{comp}(s,\mathcal{D})$ is the final computed Composite Score for dataset $\mathcal{D}$ at station s.

2.5. Diurnal Cycle and Phase Analysis

We analysed the diurnal distribution of precipitation across multiple statistical dimensions: absolute volume, phase timing, event frequency, conditional intensity, and probabilistic extremes. To distinguish volume biases from timing mismatches, this analysis was conducted in both raw and transformed spaces. Because the distributional characteristics of precipitation vary with temporal resolution, dataset-specific Box-Cox parameters were independently re-optimised for the hourly data rather than reusing the parameters derived at the daily scale. It is worth noting that Arcsinh was not used here because both CHIRPS and TAMSAT do not have sub-daily resolution.

It is critical to note that the time-shifting and Gaussian smoothing techniques applied during the daily-scale analysis were intentionally omitted from the diurnal evaluation. To preserve peaks, we did not use smoothing on the time series signals with the exception of clustering the diurnal cycle regimes; this was done to aid in distinguishing the various regimes in the region by reducing hourly noise.

2.5.1. Spatial Aggregation via Clustering

To objectively identify localised diurnal cyle regimes, we employed a probabilistic clustering approach. To extract temporal features for the clustering algorithm and ensure temporal continuity across the midnight boundary, a periodic circular smoothing filter was applied to the raw hourly precipitation profiles of both TAHMO and the evaluated GPPs.

The smoothed precipitation intensity, ${\tilde{R}}_h$, for a given local hour $h\in [0,23]$, is calculated as a centered moving average with periodic boundary conditions, as defined in Equation (12):

$${\tilde{R}}_h={\displaystyle \frac{1}{W}}{\displaystyle \sum _{i=-r}^r}{R}_{(h+i)mod24}$$

(12)

where:

• ${\tilde{R}}_h$ is the smoothed precipitation intensity at hour h.
• R represents the raw, unsmoothed hourly mean precipitation profile.
• W is the temporal smoothing window size in hours (restricted to an odd integer).
• $r={\displaystyle \frac{W-1}{2}}$ is the radius of the smoothing window.
• The modulo operator ($mod24$) enforces the circular boundary condition, guaranteeing that hour 23 connects continuously to hour 0 without edge truncation.

Prior to executing the clustering algorithm, each station’s smoothed 24-hour profile was independently Min-Max normalised (scaled between 0 and 1). This row-wise normalisation isolates the temporal phase of the diurnal cycle from the absolute volumetric magnitude.

When preparing these profiles for clustering, the smoothing window must respect the physical lifespan of local convective events. If a filter is too wide, it destroys the true weather signal, artificially flattening peaks and shifting their timing [31]. Rain events in East Africa are largely anchored to the local terrain, driven by highly localised temperature gradients such as afternoon highland heating or nighttime lake-breeze convergence [25,26]. Satellite tracking demonstrates that individual East African convective cells grow rapidly, dropping their heaviest precipitation in sudden, short-lived bursts that frequently last less than an hour [32].

Because these localised events are so temporally concentrated, applying a broad smoothing filter introduces spectral leakage across the diurnal cycle. Empirical testing within our clustering framework demonstrated that expanding the moving average to 5 hours ($W=5$) artificially merged completely distinct regimes, blurring the early morning land-breeze showers at 06:00 Local Standard Time (LST) and the afternoon sea-breeze convection (14:00 LST) observed along the Indian Ocean coastline [26] into a single, non-physical regime.

Therefore, to balance physical amplitude preservation with noise reduction, the hyperparameter grid search was bounded to a smoothing window $W\in [1,3]$ hours. By forcing the statistical model to respect these physical and methodological constraints, the Bayesian Information Criterion (BIC) objectively determined that a 3-hour smoothing window ($W=3$, $r=1$) paired with 5 distinct climate zones ($K=5$) yielded the most robust spatial configuration.

To classify the stations into their respective temporal regimes, we constructed a Shared Phase Space. Rather than clustering the datasets independently, the smoothed, normalised 24-hour diurnal profiles from TAHMO, ERA5, and IMERG were vertically concatenated to form a unified $3N\times 24$ feature matrix (where $N=117$ stations). By structuring the data this way, the clustering algorithm evaluates the pure temporal shape of each profile, entirely blind to the specific dataset from which a row originates. A Gaussian Mixture Model (GMM) was subsequently fitted to this joint dataset. This methodological choice ensures that a specific mathematical cluster (e.g., Regime 1) represents the exact same physical convective timing across all datasets, thereby enabling robust, one-to-one spatial disagreement comparisons.

The GMM evaluates each temporal profile by modelling the overall probability density as a weighted sum of K multivariate Gaussian distributions, as shown by Equation (13):

$$p\left(\mathbf{x}\right)={\displaystyle \sum _{k=1}^K}{\pi }_k\mathcal{N}\left(\mathbf{x}\right|{\mu }_k,{\Sigma }_k)$$

(13)

where:

• $p\left(\mathbf{x}\right)$ is the probability density of observing the specific diurnal profile $\mathbf{x}$.
• $\mathbf{x}$ is the 24-dimensional input vector representing the smoothed and Min-Max normalised diurnal precipitation profile for a given station.
• K is the total number of temporal regimes (optimised via Bayesian Information Criterion to $K=5$).
• ${\pi }_k$ is the mixture weight (prior probability) of regime k, representing the overall proportion of the dataset belonging to that cluster, subject to the constraint ${\sum }_{k=1}^K{\pi }_k=1$.
• $\mathcal{N}\left(\mathbf{x}\right|{\mu }_k,{\Sigma }_k)$ is the multivariate Gaussian probability density function for regime k.
• ${\mu }_k$ is the mean vector (centroid) of regime k, representing the universal 24-hour temporal shape of that specific cluster across all datasets.
• ${\Sigma }_k$ is the covariance matrix for regime k, capturing the allowable variance and temporal spread around the centroid in the 24-dimensional space.

Finally, to quantify the spatial accuracy of the GPPs, we calculated the disagreement Rate. A GPP was flagged as a “disagreement” if its dominant GMM assigned cluster differed from the TAHMO cluster at a given station. Furthermore, we defined a “severe spatial disagreement” as an instance where a GPP disagreed with the regime registered by TAHMO, yet the GMM reported a classification confidence exceeding 75%. This isolates whether GPPs are failing gracefully in uncertain zones, or failing confidently within established core regimes.

2.5.2. Volumetric and Phase Analysis (Raw Space)

For each spatial scale, we first computed the arithmetic hourly mean rainfall (${\overline{R}}_h$) for each hour $h\in [0,23]$ across the entire study period. This profile was evaluated using two distinct metrics:

• Absolute Intensity and Mass Conservation: We compared the raw diurnal curves to assess the volumetric accuracy of the GPPs. To quantify the divergence in the total daily rainfall budget, we calculated a diurnal beta (${\beta }_{\mathrm{diurnal}}$) by integrating the area under the diurnal curve (AUC) using the trapezoidal rule as illustrated by Equation (14):

  $${\beta }_{\mathrm{diurnal}}={\displaystyle \frac{{\int }_0^{24}{\overline{R}}_{\mathrm{GPP}}\left(h\right)dh}{{\int }_0^{24}{\overline{R}}_{\mathrm{obs}}\left(h\right)dh}}$$

  (14)

  where:

  –

  ${\beta }_{\mathrm{diurnal}}$ represents the diurnal mass conservation ratio, functioning as an indicator of total volumetric accuracy,

  –

  ${\overline{R}}_{\mathrm{GPP}}\left(h\right)$ and ${\overline{R}}_{\mathrm{obs}}\left(h\right)$ denote the mean precipitation intensity at local hour h for the GPP and TAHMO, respectively

  –

  The integral over the 24-hour cycle ($dh$) calculates the total daily accumulated water volume.

  A ratio of ${\beta }_{\mathrm{diurnal}}>1$ indicates a net overestimation of the daily water budget by GPP, whereas ${\beta }_{\mathrm{diurnal}}<1$ signifies a volumetric underestimation.
• Normalised Phase (Timing Analysis): To isolate phase discrepancies from amplitude biases, we applied Min-Max normalisation to the diurnal cycle curves as shown by Equation (15) below:

  $${\widehat{R}}_h={\displaystyle \frac{{\overline{R}}_h-min\left(\overline{R}\right)}{max\left(\overline{R}\right)-min\left(\overline{R}\right)}}$$

  (15)

  where:

  –

  ${\widehat{R}}_h$ represents the normalised diurnal phase strength at hour h, bounded between $[0,1]$

  –

  ${\overline{R}}_h$ is the raw mean precipitation intensity at that specific hour

  –

  $min\left(\overline{R}\right)$ and $max\left(\overline{R}\right)$ denote the minimum baseline and maximum peak intensities observed across the entire 24-hour cycle, respectively

  –

  $max\left(\overline{R}\right)-min\left(\overline{R}\right)$, calculates the total diurnal amplitude. By scaling the hourly intensities against this amplitude, the normalisation isolates the timing of the event from its absolute volumetric magnitude.

2.5.3. Frequency and Conditional Intensity(Raw Space)

Because simply calculating the arithmetic mean blends high-frequency, low-intensity drizzle with rare, extreme convective events, we decomposed the raw diurnal cycle into distinct occurrence and intensity components. For a given hour h evaluated across N days, let $R_{h,i}$ denote the raw precipitation intensity on day i, and $I(·)$ denote the indicator function, which equals 1 if the condition is met (a wet hour) and 0 otherwise (a dry hour).

• Probability of Precipitation (PoP): The diurnal frequency of rainfall occurrence, defined as the percentage of days in the dataset where a specific hour experienced rainfall exceeding a defined threshold. This is calculated using Equation (16) below:

  $$Po{P}_h={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}I(R_{h,i}>T_{threshold})\times 100$$

  (16)

  where $Po{P}_h$ is the occurrence probability (expressed as a percentage) at hour h, and $T_{threshold}$ is the physical rainfall threshold used to define a valid precipitation event (e.g., $0.1\mathrm{mm}/\mathrm{hr}$ for trace detection). This metric allows us to diagnose whether GPPs correctly capture the timing and frequency of precipitation independent of the actual volume of rain that falls.
• Conditional Mean Intensity (${\mu }_{cond,h}$): The average rainfall intensity calculated over wet hours, isolating the physical strength of the rain event from its occurrence frequency. To calculate this, we used Equation (17) below:

  $${\mu }_{cond,h}={\displaystyle \frac{{\sum }_{i=1}^N{R}_{h,i}·I(R_{h,i}>T_{threshold})}{{\sum }_{i=1}^{N}I(R_{h,i}>T_{threshold})}}$$

  (17)

  where ${\mu }_{cond,h}$ represents the average precipitation rate at hour h, computed exclusively using the occurrences where the intensity successfully exceeded $T_{threshold}$. The denominator effectively counts the number of times that specific hour h was wet across the N days, ensuring the mean is not artificially lowered by dry instances.

2.5.4. Percentiles of the Diurnal Cycle

To evaluate the central tendency of the diurnal cycle without the skewing effect of extreme right-tail outliers inherent to the arithmetic mean, we computed the percentiles of the diurnal cycle. For each hour, we extracted specific percentiles ($Q_{50},Q_{90},Q_{95},Q_{99}$). The 50th percentile ($Q_{50}$) represents the median state, while the upper percentiles gives a picture of the extreme state, providing an outlier-resistant profile of diurnal convection.

2.5.5. Temporal Accuracy (Transformed Space)

While raw space analysis highlights volumetric contributions, it often obscures timing and location accuracy due to the spatial mismatch penalties; gauge point measurements inherently register higher raw extremes than smoothed gridded pixels.

To investigate the temporal accuracy of the precipitation signal, we repeated the mean diurnal analysis in the transformed space established above. In this transformed space, extreme outliers are heavily compressed. Consequently, the arithmetic mean functions as a proxy for the geometric mean in the raw space.

Therefore, rather than evaluating the “mean volume,” the hourly mean in the transformed space evaluates the typical state of the system. This approach effectively normalises the variance mismatch between gauges and satellites, allowing us to identify whether IMERG and ERA5 successfully capture the rhythm of hourly dryness and wetness in East Africa.

3. Results

3.1. Preliminary Analysis Results and Variance Stabilisation

These are the results of the preliminary analysis from Section 2.3.2. Evaluating the temporal and intensity dynamics of the GPPs, we addressed the inherent spatial mismatch between the TAHMO gauges and the areal gridded datasets. We compared the baseline correlations of the raw time series using Nearest Neighbour against Bilinear Interpolation at the daily temporal resolution.

As shown in Table 2, bilinear interpolation systematically improved the correlation across all datasets. By interpolating the grid data to the station coordinates, we smoothed the bounding-box artifacts and reduced the spatial representation error prior to statistical transformation.

The preliminary analysis on kernel width optimisation revealed that performance gains plateaued at a standard deviation of $\sigma =2.0$ days (Figure 2).

To isolate the timing of precipitation from absolute volumetric extremes, the daily signals were transformed. We first established the optimal Box-Cox transformation parameters by minimising the skewness of each dataset.

By optimising via skewness minimisation, Box-Cox “zooms in” on the physical threshold separating dry days from trace drizzle. This forces the correlation metric to reward the GPP for correctly distinguishing “zero” from “trace” rain. As shown in Table 3, the pairwise tuning accommodates the distinct behavior of each product: ERA5’s higher offset (0.08) explicitly acknowledges its “drizzle floor,” while CHIRPS’s lesser offset (0.0001) respects its “hard zeros.”

Both variance-stabilisation techniques vastly outperformed the original linear signals (Table 4). While ArcSinh achieved the highest overall average, this performance was driven by a degradation of Box-Cox when applied to CHIRPS and TAMSAT. For these specific satellite products, Box-Cox compression resulted in lower correlations than the raw time series, whereas ArcSinh preserved their variance. On the other hand, Box-Cox was highly effective for IMERG and ERA5. Consequently, ArcSinh was utilised for CHIRPS and TAMSAT, while optimised Box-Cox was applied to IMERG and ERA5.

Figure 3 presents the network-wide distribution of station-level Pearson correlations (r) for each dataset and transformation pairing, showing the degradation of Box-Cox by CHIRPS and TAMSAT.

3.2. Temporal Phase Alignment

To investigate the impact of shifting and smoothing, we compared a global uniform, static shift (shifting + smoothing without rolling windows) against the rolling-window shifting (At each window shift + smooth).

The 30-day rolling windows proved superior to applying a uniform, static time shift across the entire time series. TAMSAT and CHIRPS benefited the most from shift alignment. Rolling window shifting outperformed uniform shifting for all the GPPs.

• TAMSAT: Rolling window shifting outperformed uniform shifting in 98.0% of windows (Mean advantage: $+0.0537$).
• CHIRPS: Rolling window shifting outperformed in 97.0% of windows (Mean advantage: $+0.0529$).
• ERA5: Rolling window shifting outperformed in 96.5% of windows (Mean advantage: $+0.0260$).
• IMERG: Rolling window shifting outperformed in 90.2% of windows (Mean advantage: $+0.0187$).

From Figure 4 and Figure 5, smoothing increases the Pearson correlation between TAHMO-GPP pairs.

As shown in Table 5, the median shift for all products is exactly 0.0 days. The standard deviations approximate $\pm 1.7$ to $2.0$ days. This shows that phase mismatch is highly dependent on the localised weather system mismatches and station location.

As shown by the Figure 6, there is an absence of distinct systematic regional clustering. While ERA5 and CHIRPS exhibit a slight, domain-wide negative tendency, IMERG and TAMSAT exhibit a slight positive tendency.

3.3. Daily Scale Performance and Ranking

The datasets were ranked using the multi-component Composite Score derived from the rolling window analysis.

IMERG demonstrated superior capability, achieving the highest overall Composite Score (0.6059) as shown in Table 6, driven by better timing and volume accuracy, and rain event detection.

A definitive insight from this ranking analysis is the performance paradox of ERA5. Despite achieving the second-highest correlation, ERA5 registered the lowest Normalised Bias Score (0.6389) of all datasets. To interpret this, the threshold mechanics and exclusion logic of our methodology must be considered. In the optimised Box-Cox space of ERA5, the physical threshold defining a precipitation event (2.5 mm) translates to 0.9455 Box-Cox units. Our Bias Score is calculated conditionally: it evaluates intensity only when the GPP and the ground station both agree on the precipitation state.

If one dataset records a value above the dry threshold (Positive State) while the other falls below it (Negative State), the resulting ratio becomes negative. Because the logarithm of a negative number is undefined, these “state mismatches” are algorithmically excluded from the conditional Bias Score. Therefore, ERA5’s low Bias Score is not a simple volumetric error. It frequently triggers low-intensity precipitation during physically dry periods, creating state mismatches with TAHMO.

3.4. Diurnal Cycle Results

3.4.1. Diurnal Cycle Regime Clustering

Applying the GMM to the shared phase space isolated $K=5$ distinct, physically robust diurnal regimes across East Africa. Figure 7 shows the peaks of these diurnal regimes as well as their spatial mapping.

The spatial mapping of TAHMO reveals that these clusters map to known mechanisms:

• Regime 2 (Peak 06:00 LST) & Regime 4 (Peak 14:00 LST): These regimes geographically correspond to the well-documented land-water breeze circulations of East Africa. The morning precipitation peak (Regime 2) is clustered primarily over the open water and adjacent shores of the Indian Ocean and Lake Victoria. Literature largely attributes this specific morning regime to nocturnal land breezes converging over the water [26]. On the other hand, the early afternoon peak (Regime 4) is located slightly further inland, spatially aligning with the documented inland penetration of the daytime sea and lake breeze fronts [26].
• Regime 1 (Peak 16:00 LST): Representing classic afternoon convection, this regime is driven by daytime surface heating and subsequent atmospheric instability over continental terrain [26,33]. Within the observational network, this 16:00 LST footprint is predominantly captured across Western Kenya, scattered Eastern Tanzanian stations, East of Lake Kivu, the region around Kigali and South of Kigali, and Western Uganda.
• Regime 0 (Peak 19:00 LST): Clustered predominantly around areas of steep orography, this regime captures delayed evening precipitation. In Kenya, this spatial footprint dominates the North Rift and the Kenyan Highlands, with footprints southward around the Kenya-Tanzania border. In Uganda, it forms a distinct northeastern corridor extending from West Pian Upe Game Reserve towards Gulu. In Rwanda, it is highly localised to the regions north and south of Lake Kivu, specifically aligning with the Virunga Mountains and the steep terrain west of Nyungwe National Park. The spatial isolation of this 19:00 LST peak is consistent with observations of rain events initiating over mountain peaks during the afternoon and subsequently producing precipitation in the adjacent valleys and lower elevations during the evening [26].
• Regime 3 (Peak 01:00 LST): This nocturnal precipitation regime is heavily clustered in central Kenya, with scattered signatures in northeastern Kenya, along the coast, and specific corridors in northwestern Uganda. Notably, it is absent from the Tanzanian and Rwandan stations. The geographic footprint of this 01:00 LST peak strongly corresponds to the documented locations of mature Mesoscale Convective Systems (MCSs) that produce late-night precipitation maximums over inland plains [26].

While the GMM successfully isolated these physical regimes within the unified feature matrix, projecting the GPPs onto this space revealed substantial structural divergence. Across the 117 evaluated stations, the datasets achieved concurrent alignment (where both ERA5 and IMERG matched TAHMO) at only 50 stations, resulting in a network-wide divergence rate of 57.3%.

Relative to TAHMO, IMERG demonstrated higher spatial consistency, registering deviation at 35.0% (41 stations) compared to ERA5’s deviation rate of 50.4% (59 stations). Furthermore, significant inter-product divergence exists between the gridded datasets themselves; ERA5 and IMERG classified precipitation into entirely distinct diurnal regimes at 43 separate stations. A breakdown of the spatial mapping highlights the following dataset-specific structural tendencies:

• ERA5 Systemic Omissions: ERA5 exhibits parameterisation failures over complex terrain. It completely fails to identify Regime 0 across Uganda and Rwanda, defaulting almost entirely to the generic afternoon convection of Regime 1. Furthermore, ERA5 completely misses the Regime 3 nocturnal signatures in Uganda and classifies the majority of them in central Kenya as either Regime 0 or 1.
• IMERG Blurring: While IMERG captures the broader regional distributions better than ERA5, it struggles in certain regions. For example, IMERG frequently confuses Regime 0 and Regime 1 across the Kenya-Tanzania border (classifying only a single station near Arusha as Regime 0) as well as in the Kenyan highlands and North Rift. It also disagrees with TAHMO in stations around Lake Kivu in Rwanda and Murchison Falls in Uganda.

Out of ERA5’s 59 disagreeing stations, 57 were flagged as “severe spatial disagreements”—meaning the GPP chose an entirely convective regime but did so with a confidence exceeding 75%. Similarly, 38 of IMERG’s 41 disagreements were severe. This indicates that the GPPs are not failing gracefully in uncertain areas; rather, their underlying physics and parameterisation schemes are confidently simulating the wrong diurnal convective mechanisms across vast portions of the East African region.

3.4.2. Volumetric and Phase Analysis (Raw Space)

The domain-wide analysis of absolute intensity suggests a moderate volumetric underestimation by the GPPs, with IMERG and ERA5 yielding diurnal mass conservation ratios (${\beta }_{\mathrm{diurnal}}$) of 0.75 and 0.80, respectively, as shown in Figure 8. However, disaggregating the analysis into regional sub-domains reveals spatial heterogeneity that the domain-wide mean obscures.

The normalised phase analysis revealed that, on a domain-wide scale, TAHMO, ERA5 and IMERG all peak at 16:00 LST. However, spatial disaggregation demonstrates that these convective peaks, alongside the volumetric biases, vary drastically depending on the localised climatological regime (Table 7).

As shown in Table 7, the GPPs severely underestimate the total volume in Regime 0 (missing more than 50% of the water budget). The GPPs substantially miss the morning convection in Regime 0.

ERA5 overestimates the afternoon rainfall in Regimes 2,3 and 4 while underestimating the morning rainfall in those regimes as well as Regime 0.

In Regime 4, both GPPs overestimate the daily water budget (IMERG $\beta =1.19$, ERA5 $\beta =1.21$) as shown in Figure 9, likely struggling to correctly parameterise the intensity of the nocturnal land-breeze convergence [26].

More results focusing on seasonal aspects of this investigation can be found in Appendix A.

3.4.3. Occurrence vs. Intensity

Decomposing the arithmetic mean into Probability of Precipitation (PoP) and Conditional Intensity exposes the flaws of the GPPs, especially within ERA5.

At the threshold of $>0.1$ mm/hr, ERA5 exhibits a massive overestimation of rainfall frequency, reaching a PoP of 35% at 16:00 LST. In contrast, TAHMO hovers below 7% while IMERG remains below a 22% probability of occurrence at this threshold. However, when evaluating the Conditional Mean Intensity of these wet hours, TAHMO dominates. TAHMO registers sharp, high-magnitude spikes (e.g., 4.3 mm/hr and 5mm/hr at 08:00 LST and 12:00 LST respectively, and multiple nocturnal bursts $>2.9$ mm/hr), whereas ERA5’s conditional intensity rarely exceeds 1.1 mm/hr, as shown in Figure 10.

This inverse relationship—where ERA5 generates rainfall almost half the time but at highly suppressed magnitudes—confirms a “drizzle bias” inherent to its convective parameterisation scheme. It satisfies its daily volumetric water budget through persistent, light rain rather than resolving the discrete, violent mesoscale convective systems captured by the gauges.

An exhaustive breakdown of PoP and conditional intensities across certain physical thresholds and across the regimes is provided in Appendix A.

3.4.4. Percentiles of the Diurnal Cycle

The percentile analysis further validates the discrete nature of East African convection. As shown in Figure 11, across all regions, the median state ($Q_{50}$) for TAHMO and IMERG is completely dry (0.0 mm/hr), physically reflecting that rain does not occur on the majority of days at any given hour. ERA5, on the other hand, registers precipitation (trace rain) even at the 50th percentile, further exposing its over-triggering. The true convective signals for TAHMO and IMERG only emerge in the upper percentiles ($Q_{90}$ to $Q_{99}$), showing that regional precipitation is driven by rare, heavy-tail events rather than central-tendency means.

Percentiles of each regime are provided in Appendix A.

4. Discussion

4.1. Overcoming Spatial and Distributional Mismatches

The foundational step of our analysis highlighted a penalty of spatial support mismatch between areal and point-scale datasets as shown in the systematic correlation improvement achieved via bilinear interpolation (Table 2).

However, addressing spatial scale alone is insufficient due to the zero-inflated, heavy-tailed nature of precipitation. The transition into transformed spaces (Box-Cox and ArcSinh) proved critical (Table 4). By compressing extreme right-tail outliers, these transformations force the evaluation metrics to prioritise the event detection of rain vs. no-rain rather than being overwhelmingly skewed by a single extreme event.

Practically, this revealed that GPPs possess fundamentally different underlying distributions. The performance of Box-Cox when applied to CHIRPS and TAMSAT demonstrates that they produce “hard zeros” which are distorted by aggressive power transforms. On the other hand, ERA5 and IMERG require the specific parameterisation of Box-Cox to effectively stretch the boundary between dry days and their inherent trace-level noise. For hydrologists and machine learning practitioners, this indicates that applying a uniform statistical transformation across different precipitation products will inherently degrade the signals of some while enhancing others; transformation strategies must be dynamically tailored to the dataset’s specific climatological behavior.

4.2. The Physical Basis for Shifting Temporal Alignment

The comparison between static uniform shifting and 30-day rolling-window analysis proved that temporal displacement in GPPs is non-stationary. As demonstrated in Table 5, while the median shift across the datasets is 0.0 days, the standard deviations approximate $\pm 1.7$ to $2.0$ days. This indicates that phase mismatch is highly dependent on localised weather systems rather than a systematic global delay.

The findings from the hourly diurnal analysis provide the fundamental physical explanation for this daily-scale stochasticity. The diurnal phase analysis revealed that actual precipitation exhibits complex, multi-modal timing dictated by local topography as shown by the phase peaks and how the GPPs’ peaks align to TAHMO’s peaks in Table 7.

4.3. Diagnosing Parameterisation Flaws

By evaluating the datasets in an optimised transformed space, we were able to rank the GPPs based on their true physical accuracy (Table 6). While IMERG demonstrated superior correlation and volume accuracy (${\mu }_{\alpha }=1.0514$), the ranking analysis exposed a performance paradox within ERA5. Despite achieving the second-highest optimal correlation ($0.6306$), ERA5 registered the lowest Bias Score ($0.6389$) of all datasets. The diurnal cycle analysis definitively explains this paradox. At the hourly scale at a 0.1 mm/hr threshold, ERA5 exhibits a massive Probability of Precipitation (PoP of 35% at 16:00 LST), yet its Conditional Mean Intensity rarely exceeds 1.1 mm/hr. In contrast, TAHMO rains much less frequently but at highly concentrated, intense magnitudes (e.g., bursts $>2.5$ mm/hr) Figure 10.

In our daily-scale methodology, the conditional Bias Score utilises an exclusion logic that removes “State mismatches”—instances where one dataset records a dry state and the other records a wet state. ERA5’s low Bias Score exposes its “drizzle bias.” Because it frequently reports low-intensity precipitation during physically dry periods, forcing persistent state mismatches with TAHMO. ERA5 satisfies its total daily volumetric budget by drizzling lightly across the entire day.

4.4. Event Detection and Accuracy

IMERG demonstrated superior overall timing and volume accuracy. This conclusion strongly aligns with recent regional evaluations utilising the TAHMO network [34], who similarly identified IMERG as more reliable when compared to datasets like CHIRPS and ERA5. At the daily scale, IMERG and ERA5 achieved the highest rain event detection skill—answering the question, `Did it rain today?’—yielding point-to-point CSI scores of 0.3481 and 0.3068, respectively. On the other hand, CHIRPS and TAMSAT have lower detection skill with point-to-point CSI of 0.2914 and 0.2898, respectively, indicating a much higher rate of missed events and false alarms. The variance scores further highlight the performance of CHIRPS and TAMSAT, where lower values indicate a failure to capture the variability of rainfall. While IMERG and ERA5 show higher variance scores (0.8581 and 0.8145 respectively), CHIRPS and TAMSAT have lower scores of 0.6874 and 0.6958 (Table 6). Previous studies utilising the TAHMO network [20] have established that CHIRPS suffers from volumetric overestimations in East Africa using traditional raw-space metrics. Our variance-stabilised analysis builds upon this by revealing that CHIRPS does not simply overestimate rainfall amounts; it exhibits poor rain event detection and does not align them in time as well as IMERG or ERA5.

Furthermore, isolating the geometric mean of the bias (${\mu }_{\beta }$) and variability ratios (${\mu }_{\alpha }$) within the transformed space shows how this failure happens. As shown in Table 6, CHIRPS registers a geometric mean bias ratio of 1.8551 combined with a high geometric mean variability ratio of 1.6933 in the transformed space. Because the ArcSinh transformation actively compresses extreme right-tail outliers, this $85\%$ overestimation of the geometric mean suggests that the regional wet bias in CHIRPS is driven by systematic over-estimation of the typical state of day-to-day precipitation while producing spurious fluctuations between wet and dry states.

Operational Implications

Distinguishing between detection skill and volumetric accuracy carries profound implications for operational hydrology and climate downscaling in East Africa. The high event detection scores achieved by IMERG and ERA5 indicate that their internal parameterisations are highly optimised for occurrence-based applications. For tasks such as agricultural drought monitoring, crop-yield forecasting, and categorical anomaly detection, successfully capturing the duration of dry spells or the exact frequency of wet days is fundamentally more critical than estimating absolute rainfall intensity. In these contexts, both datasets provide highly robust utility.

Beyond basic event detection, impact models requiring strict mass conservation and dynamic phase accuracy—such as river routing, flash flood forecasting, and fine-scale atmospheric downscaling—demand rigorous volumetric precision. In these volume-dependent applications, relying on datasets with much lower detection skill such as CHIRPS and TAMSAT introduces substantial risk. IMERG maintains high detection skill while simultaneously preserving the heavy-tailed variance and the central tendency of the physical water budget, thus emerging as the optimal choice for these specific tasks.

5. Conclusion

This study demonstrates that evaluating GPPs in East Africa requires moving beyond traditional raw-space, static comparisons. Due to the heavy-tailed, zero-inflated nature of precipitation and the spatial support mismatch between point-scale gauges and areal grids, applying dataset-specific variance-stabilising transformations (Box-Cox and ArcSinh) is useful in finding spatial and timing patterns. Furthermore, our diurnal phase analysis revealed that temporal mismatches between GPPs and ground observations are not systematic global delays, but rather random artifacts driven by local topography and convective parameterisations.

By isolating timing accuracy from volume accuracy, our analysis reveals that each GPP should be utilised based on its strengths. IMERG achieved the highest overall timing and volume accuracy, establishing it as suitable for volume-dependent tasks like flash flood modeling. Meanwhile, despite exhibiting high baseline correlations, ERA5 is characterised by a drizzle bias. It satisfies its diurnal water budget through persistent, low-intensity drizzle while failing to capture the intense, discrete convective bursts that characterise East African hydroclimatology.

These findings carry implications for regional climate impact studies and the development of AI-driven weather emulators. First, utilising raw GPPs, especially ERA5 precipitation directly for local hydrological modelling without bias correction or post-processing is highly discouraged. Second, when training machine learning downscaling models, developers should consider how they evaluate the models. Standard evaluation metrics like RMSE rely on strict point-to-point comparisons, meaning they heavily penalise minor timing errors. If a model accurately simulates a rain event but places it two hours early, RMSE punishes it twice: once for a false alarm, and once for the missed event. These models should instead be rewarded for predicting the timing, volume and intensity of localised rain events.

Appendix A.1. Spatial Support Mismatch and the Normalisation Artifact

As illustrated in the domain-wide evaluation (Figure 8), a distinct “order flip” occurs between the raw intensity and the normalised phase profiles during the early morning hours (e.g., 01:00 to 04:00 LST). In the absolute intensity space, the TAHMO gauge network records higher physical precipitation volumes than the IMERG satellite product. However, when transitioning to the Min-Max normalised in the transformed space, this hierarchy inverts, and IMERG appears to exhibit a much stronger relative morning phase than TAHMO.

This phenomenon is not an error in the observational data, but rather a direct mathematical consequence of spatial support mismatch. The phase normalisation maps the diurnal cycle to a $[0,1]$ scale by dividing the hourly intensities by the total diurnal amplitude ($max\left(\overline{R}\right)-min\left(\overline{R}\right)$).

As demonstrated by Peleg et al. [35] , point-scale gauges inherently capture the extreme, highly localised cores of convective rain. When it rains over a TAHMO station, the gauge registers a massive, unsmoothed afternoon peak, resulting in an exceptionally large $max\left(\overline{R}\right)$ denominator. In contrast, areal gridded products like IMERG (representing a roughly $10\times 10$ km pixel on the ground) spatially average these sub-grid convective cores with the surrounding non-precipitating areas, artificially dampening their diurnal maximum.

Because TAHMO’s normalisation denominator is physically much larger, its non-peak nocturnal hours are compressed toward zero. IMERG, possessing a dampened afternoon peak and thus a smaller diurnal range, allocates a larger percentage of its total daily variance to the morning hours.

Therefore, the observed “flip” in Figure 8 does not outrightly indicate that IMERG measures more morning rain; rather, it highlights the spatial variance penalty.

Appendix A.2. Regime-Level Probabiliy of Precipitation and Conditional Mean Intensity

The figures below highlight the diurnal cycle considered at different precipitation thresholds for the various regimes.

Figure A1 below shows the probability of precipitation and conditional mean intensity for regime 0.

Figure A1. Regime 0 diurnal cycle PoP and mean intensity at various precipitation thresholds.

Figure A1. Regime 0 diurnal cycle PoP and mean intensity at various precipitation thresholds.

Figure A2 below shows the probability of precipitation and conditional mean intensity for regime 1.

Figure A2. Regime 1 diurnal cycle PoP and mean intensity at various precipitation thresholds.

Figure A2. Regime 1 diurnal cycle PoP and mean intensity at various precipitation thresholds.

Figure A3 below shows the probability of precipitation and conditional mean intensity for regime 2.

Figure A3. Regime 2 diurnal cycle PoP and mean intensity at various precipitation thresholds.

Figure A3. Regime 2 diurnal cycle PoP and mean intensity at various precipitation thresholds.

Figure A4 below shows the probability of precipitation and conditional mean intensity for regime 3.

Figure A4. Regime 3 diurnal cycle PoP and mean intensity at various precipitation thresholds.

Figure A4. Regime 3 diurnal cycle PoP and mean intensity at various precipitation thresholds.

Figure A5 below shows the probability of precipitation and conditional mean intensity for regime 4.

Figure A5. Regime 4 diurnal cycle PoP and mean intensity at various precipitation thresholds.

Figure A5. Regime 4 diurnal cycle PoP and mean intensity at various precipitation thresholds.

Appendix A.3. Seasonal Diurnal Cycle at Each Regime

The results shown in Figure A6 are those of Timing and Volume in Section 2.5.2 where we look at the mean diurnal cycle at seasonal levels. It is also influenced by the normalisation in Appendix A.1, in the transformed space.

Figure A6. Seasonal diurnal cycle phase alignment analysis in all the regimes.

Figure A6. Seasonal diurnal cycle phase alignment analysis in all the regimes.

The results indicate that the 19:00 LST in Regime 0’s mean diurnal cycle is largely as a result of the October-November-December short rains season.

It also shows that even within the regimes, there is a difference in the precipitation experienced in different seasons. For instance, in Regime 0, the absolute mean intensity during OND is on average greater than that of MAM, using TAHMO as our reference.

Appendix A.4. Diurnal Cycle Percentiles

This section shows the results of the percentile analysis at each regime.

Figure A7. Sub-regional percentiles of the diurnal cycle.

Figure A7. Sub-regional percentiles of the diurnal cycle.

The median state even at the regime level shows the ERA5 “drizzle bias” in all the regimes as shown by Figure A7.